Wygeneruj zmienną losową ze zdefiniowaną korelacją z istniejącą zmienną (zmiennymi)

71

Dla badań symulacyjnych mam do generowania zmiennych losowych, które wykazują prefined (populacji) korelację do istniejącej zmiennej .Y

I spojrzał w Ropakowaniach copula, a CDVinektóre mogą powodować przypadkowe wielowymiarowych rozkładów danej struktury zależności. Nie można jednak naprawić jednej z powstałych zmiennych do istniejącej zmiennej.

Wszelkie pomysły i linki do istniejących funkcji są mile widziane!


Wniosek: pojawiły się dwie ważne odpowiedzi, z różnymi rozwiązaniami:

  1. R Scenariusz według caracal, co wylicza zmienną losową z dokładnym (próbki) korelacji ustalonej zmiennej
  2. R Funkcja znalazłem się, co wylicza zmienną losową o określonej populacji korelacji do predefiniowanej zmiennej

[Dodanie @ttnphns: Zezwoliłem na rozszerzenie tytułu pytania z pojedynczej stałej zmiennej na dowolną liczbę stałych zmiennych; tj. jak wygenerować zmienną mającą predefiniowane korelacje z pewnymi stałymi, istniejącymi zmiennymi]

Felix S.
źródło
2
Zobacz to powiązane pytanie stats.stackexchange.com/questions/13382/..., które bezpośrednio odnosi się do twojego pytania (przynajmniej od strony teoretycznej).
Makro
Poniższe Q jest również silnie powiązane i będzie interesujące: Jak wygenerować skorelowane liczby losowe (podane średnie wariancje i stopień korelacji) .
gung

Odpowiedzi:

56

Oto kolejny: dla wektorów ze średnią 0 ich korelacja jest równa cosinus ich kąta. Zatem jednym sposobem na znalezienie wektora o dokładnie pożądanej korelacji , odpowiadającej kątowi :r θxrθ

  1. uzyskaj ustalony wektor i losowy wektorx 2x1x2
  2. wyśrodkuj oba wektory (średnia 0), dając wektory , ˙ x 2x˙1x˙2
  3. uczyń prostopadłą do (rzut na podprzestrzeń ortogonalną), dając ˙ x 1 ˙ x 2x˙2x˙1x˙2
  4. przeskaluj i do długości 1, dając i ˙ x 2 ˉ x 1 ˉ x 2x˙1x˙2x¯1x¯2
  5. ˉ x 1θ ˉ x 1rx1x¯2+(1/tan(θ))x¯1 to wektor, którego kąt do to , a którego korelacja z jest więc . Jest to również korelacja z ponieważ transformacje liniowe pozostawiają korelację bez zmian.x¯1θx¯1rx1

Oto kod:

n     <- 20                    # length of vector
rho   <- 0.6                   # desired correlation = cos(angle)
theta <- acos(rho)             # corresponding angle
x1    <- rnorm(n, 1, 1)        # fixed given data
x2    <- rnorm(n, 2, 0.5)      # new random data
X     <- cbind(x1, x2)         # matrix
Xctr  <- scale(X, center=TRUE, scale=FALSE)   # centered columns (mean 0)

Id   <- diag(n)                               # identity matrix
Q    <- qr.Q(qr(Xctr[ , 1, drop=FALSE]))      # QR-decomposition, just matrix Q
P    <- tcrossprod(Q)          # = Q Q'       # projection onto space defined by x1
x2o  <- (Id-P) %*% Xctr[ , 2]                 # x2ctr made orthogonal to x1ctr
Xc2  <- cbind(Xctr[ , 1], x2o)                # bind to matrix
Y    <- Xc2 %*% diag(1/sqrt(colSums(Xc2^2)))  # scale columns to length 1

x <- Y[ , 2] + (1 / tan(theta)) * Y[ , 1]     # final new vector
cor(x1, x)                                    # check correlation = rho

wprowadź opis zdjęcia tutaj

Na rzut prostopadły , to używany Q R -decomposition poprawić stabilność liczbową, gdyż wtedy po prostu P = Q, Q ' .PQRP=QQ

karakal
źródło
Próbowałem przepisać kod na składnię SPSS. Potykam się o twój rozkład QR, który zwraca kolumnę 20x1. W SPSS mam ortonormalizację Gram-Schmidta (która jest również rozkładem QR), ale nie mogę odtworzyć wynikowej kolumny Q. Czy możesz mi przeżuć swoją akcję QR. Lub wskaż kilka obejść, aby uzyskać projekcję. Dzięki.
ttnphns
@ caracal, P <- X %*% solve(t(X) %*% X) %*% t(X)nie produkuje r = 0,6, więc to nie jest obejście. Wciąż jestem zdezorientowany. (Z przyjemnością naśladuję twój wyraz Q <- qr.Q(qr(Xctr[ , 1, drop=FALSE]))w SPSS, ale nie wiem jak.)
ttnphns
@ttnphns Przepraszam za zamieszanie, mój komentarz dotyczył ogólnej sprawy. Zastosowanie do sytuacji w przykładzie: Uzyskanie macierzy projekcji poprzez rozkład QR jest tylko dla stabilności numerycznej. Można uzyskać macierz projekcji jako , jeśli podprzestrzeń jest łączonych przez kolumny macierzy . W R możesz tutaj pisać, ponieważ podprzestrzeń jest rozłożona na pierwszą kolumnę . Macierz rzutu na dopełniacz ortogonalny to wtedy IP. P=X(XX)1XXXctr[ , 1] %*% solve(t(Xctr[ , 1]) %*% Xctr[ , 1]) %*% t(Xctr[ , 1])Xctr
caracal
4
Czy ktoś mógłby wyjaśnić, jak wykonać coś podobnego dla więcej niż dwóch próbek? Powiedzmy, jeśli chciałem 3 próbki skorelowane parami przez rho, jak mogę przekształcić to rozwiązanie, aby to osiągnąć?
Andre Terra,
w przypadku limitu rho=1uznałem za użyteczne zrobić coś takiego: if (isTRUE(all.equal(rho, 1))) rho <- 1-10*.Machine$double.epsw przeciwnym razie dostawałem NaNs
PatrickT
19

Opiszę najogólniejsze możliwe rozwiązanie. Rozwiązanie problemu w tej ogólności pozwala nam osiągnąć niezwykle kompaktową implementację oprogramowania: wystarczy tylko dwie krótkie linie Rkodu.

Wybierz wektor o tej samej długości co , zgodnie z dowolnym rozkładem. Niech być pozostałości z regresji metodą najmniejszych kwadratów z z : ten wyodrębnia elementu z . Przez ponowne dodanie odpowiedniego wielokrotność do możemy wytworzenia wektora posiadającego dowolną korelacji z . Rozwiązaniem jest dowolna dowolna stała addytywna i dodatnia stała mnożąca - którą możesz dowolnie wybraćY Y X Y Y X Y Y ρ YXYYXYYXYYρY

XY;ρ=ρSD(Y)Y+1ρ2SD(Y)Y.

(„ nazwa ” oznacza wszelkie obliczenia proporcjonalne do odchylenia standardowego.)SD


Oto działający Rkod. Jeśli nie podasz , kod pobierze swoje wartości ze standardowego rozkładu normalnego na wielu odmianach.X

complement <- function(y, rho, x) {
  if (missing(x)) x <- rnorm(length(y)) # Optional: supply a default if `x` is not given
  y.perp <- residuals(lm(x ~ y))
  rho * sd(y.perp) * y + y.perp * sd(y) * sqrt(1 - rho^2)
}

W celu zilustrowania, że generowane losowo z elementów, a wytwarzane o różnych określonych korelacji z tym . Wszystkie zostały utworzone przy użyciu tego samego wektora początkowego . Oto ich wykresy rozrzutu. „Wykresy rugowe” u dołu każdego panelu pokazują wspólny wektor50 X Y ; ρ Y X = ( 1 , 2 , , 50 ) YY50XY;ρYX=(1,2,,50)Y

Postać

Istnieje niezwykłe podobieństwo między fabułami, czyż nie :-).


Jeśli chcesz eksperymentować, oto kod, który wygenerował te dane i rysunek. (Nie zawracałem sobie głowy skorzystaniem ze swobody, aby przesuwać i skalować wyniki, które są łatwymi operacjami).

y <- rnorm(50, sd=10)
x <- 1:50 # Optional
rho <- seq(0, 1, length.out=6) * rep(c(-1,1), 3)
X <- data.frame(z=as.vector(sapply(rho, function(rho) complement(y, rho, x))),
                rho=ordered(rep(signif(rho, 2), each=length(y))),
                y=rep(y, length(rho)))

library(ggplot2)
ggplot(X, aes(y,z, group=rho)) + 
  geom_smooth(method="lm", color="Black") + 
  geom_rug(sides="b") + 
  geom_point(aes(fill=rho), alpha=1/2, shape=21) +
  facet_wrap(~ rho, scales="free")

BTW, ta metoda z łatwością uogólnia na więcej niż jedno : jeśli jest to matematycznie możliwe, znajdzie po określeniu korelacji z całością zestaw . Wystarczy użyć zwykłych najmniejszych kwadratów, aby wyjąć efekty wszystkich z i utworzyć odpowiednią liniową kombinację i reszt. (Pomaga to zrobić w kategoriach podwójnej podstawy dla , która jest uzyskiwana przez obliczenie pseudo-odwrotności. Poniższy kod używa SVD dla osiągnięcia tego.)X Y 1 , Y 2 , , Y k ; ρ 1 , ρ 2 , , ρ k Y i Y i X Y i Y YYXY1,Y2,,Yk;ρ1,ρ2,,ρkYiYiXYiYY

Oto szkic algorytmu, w Rktórym są podane jako kolumny macierzy :Yiy

y <- scale(y)             # Makes computations simpler
e <- residuals(lm(x ~ y)) # Take out the columns of matrix `y`
y.dual <- with(svd(y), (n-1)*u %*% diag(ifelse(d > 0, 1/d, 0)) %*% t(v))
sigma2 <- c((1 - rho %*% cov(y.dual) %*% rho) / var(e))
return(y.dual %*% rho + sqrt(sigma2)*e)

Poniżej znajduje się pełniejsza implementacja dla tych, którzy chcieliby eksperymentować.

complement <- function(y, rho, x) {
  #
  # Process the arguments.
  #
  if(!is.matrix(y)) y <- matrix(y, ncol=1)
  if (missing(x)) x <- rnorm(n)
  d <- ncol(y)
  n <- nrow(y)
  y <- scale(y) # Makes computations simpler
  #
  # Remove the effects of `y` on `x`.
  #
  e <- residuals(lm(x ~ y))
  #
  # Calculate the coefficient `sigma` of `e` so that the correlation of
  # `y` with the linear combination y.dual %*% rho + sigma*e is the desired
  # vector.
  #
  y.dual <- with(svd(y), (n-1)*u %*% diag(ifelse(d > 0, 1/d, 0)) %*% t(v))
  sigma2 <- c((1 - rho %*% cov(y.dual) %*% rho) / var(e))
  #
  # Return this linear combination.
  #
  if (sigma2 >= 0) {
    sigma <- sqrt(sigma2) 
    z <- y.dual %*% rho + sigma*e
  } else {
    warning("Correlations are impossible.")
    z <- rep(0, n)
  }
  return(z)
}
#
# Set up the problem.
#
d <- 3           # Number of given variables
n <- 50          # Dimension of all vectors
x <- 1:n         # Optionally: specify `x` or draw from any distribution
y <- matrix(rnorm(d*n), ncol=d) # Create `d` original variables in any way
rho <- c(0.5, -0.5, 0)          # Specify the correlations
#
# Verify the results.
#
z <- complement(y, rho, x)
cbind('Actual correlations' = cor(cbind(z, y))[1,-1],
      'Target correlations' = rho)
#
# Display them.
#
colnames(y) <- paste0("y.", 1:d)
colnames(z) <- "z"
pairs(cbind(z, y))
Whuber
źródło
YBTW, this method readily generalizes to more... Just use ordinary least squares... and form a suitable linear combination
1
@ttnphns Zrobiłem to.
whuber
1
Dziękuję bardzo! Rozumiem i zakodowałem dziś twoje podejście w SPSS. Naprawdę świetna propozycja. Nigdy nie myślałem, że pojęcie podwójnej podstawy ma zastosowanie do rozwiązania zadania.
ttnphns
Czy można zastosować podobne podejście, aby uzyskać jednolicie rozłożony wektor? To znaczy, mam istniejący wektor xi chcę wygenerować nowy wektor yskorelowany z, xale także chcę, aby ywektor był równomiernie rozłożony.
Skumin
@Skumin Zastanów się nad użyciem do tego kopuły, aby kontrolować relacje między dwoma wektorami.
whuber
6

Oto inne podejście obliczeniowe (rozwiązanie zostało zaadaptowane z postu na forum Enrico Schumanna). Według Wolfganga (patrz komentarze) jest to obliczeniowo identyczne z rozwiązaniem zaproponowanym przez ttnphns.

W przeciwieństwie do rozwiązania karakala nie wytwarza próbki o dokładnej korelacji , ale dwa wektory, których korelacja populacji jest równa .ρρρ

Poniższa funkcja może obliczyć dwuwymiarowy rozkład próbek pobranych z populacji o danym . Oblicza dwie zmienne losowe lub pobiera jedną istniejącą zmienną (przekazaną jako parametr ) i tworzy drugą zmienną o pożądanej korelacji:ρx

# returns a data frame of two variables which correlate with a population correlation of rho
# If desired, one of both variables can be fixed to an existing variable by specifying x
getBiCop <- function(n, rho, mar.fun=rnorm, x = NULL, ...) {
     if (!is.null(x)) {X1 <- x} else {X1 <- mar.fun(n, ...)}
     if (!is.null(x) & length(x) != n) warning("Variable x does not have the same length as n!")

     C <- matrix(rho, nrow = 2, ncol = 2)
     diag(C) <- 1

     C <- chol(C)

     X2 <- mar.fun(n)
     X <- cbind(X1,X2)

     # induce correlation (does not change X1)
     df <- X %*% C

     ## if desired: check results
     #all.equal(X1,X[,1])
     #cor(X)

     return(df)
}

Funkcja może również wykorzystywać niestandardowe rozkłady brzeżne poprzez dostosowanie parametru mar.fun. Należy jednak pamiętać, że ustalenie jednej zmiennej tylko wydaje się działać z zmiennej o rozkładzie normalnym x! (co może odnosić się do komentarza Makra).

Należy również zauważyć, że „mały współczynnik korygujący” z pierwotnego postu został usunięty, ponieważ wydaje się, że przesądza powstałe korelacje, przynajmniej w przypadku rozkładów Gaussa i korelacji Pearsona (patrz także komentarze).

Felix S.
źródło
ρ
1
Łatwo jest wykazać, że z wyjątkiem „małej korekty do rho” (której cel w tym kontekście umyka mi), jest to dokładnie to samo, co sugerowane wcześniej przez ttnphns. Metoda opiera się po prostu na rozkładzie Choleskiego macierzy korelacji w celu uzyskania pożądanej macierzy transformacji. Patrz na przykład: en.wikipedia.org/wiki/… . I tak, to da ci tylko dwa wektory, których korelacja populacji jest równa rho.
Wolfgang,
„Mała korekta do rho” była w oryginalnym poście i została opisana tutaj . Właściwie to tak naprawdę nie rozumiem; ale badanie 50000 symulowanych korelacji z rho = .3 pokazuje, że bez „małej korekty” powstaje średnia z r .299, podczas gdy z korektą średnia z .312 (co jest wartością skorygowanego rho) wynosi wytworzony. Dlatego usunąłem tę część z funkcji.
Felix S
Wiem, że to jest stare, ale chcę również zauważyć, że ta metoda nie będzie działać w przypadku dodatnich macierzy korelacji. Np. Korelacja -1.
zzk
1
Dzięki; Zauważyłem, że jeśli x1 nie jest znormalizowana średnia = 0, sd = 1, a nie wolę przeskalować go, trzeba zmodyfikować linię: X2 <- mar.fun(n)do X2 <- mar.fun(n,mean(x),sd(x))zdobycia pożądanego korelację między x1 i x2
Dave M
6

XYXrXrY=rX+EE0sd=1r2XYrXYXρ=r

rEXEXYX1,X2,X3,...

XrYYrY


Zaktualizuj 11 listopada 2017 r. Dzisiaj spotkałem ten stary wątek i postanowiłem rozszerzyć moją odpowiedź, pokazując algorytm iteracyjnego dopasowania, o którym mówiłem na początku.

Y X

Disclamer: To iteracyjne rozwiązanie, które znalazłem gorsze od doskonałego, oparte na znalezieniu podwójnej podstawy i zaproponowane przez @whuber w tym wątku dzisiaj. @ rozwiązanie Whubera nie jest iteracyjne i, co ważniejsze, wydaje mi się, że wpływa na wartości wejściowej zmiennej „świnia” nieco mniej niż algorytm „mój” (byłoby to atutem, gdyby zadaniem było „poprawić” istniejąca zmienna i nie generować losowych zmiennych od zera). Nadal publikuję moje z ciekawości i dlatego, że to działa (patrz także przypis).

X1,X2,...,XmYYr1,r2,...,rmX

YXYY

  1. rdf=n1Sj=rjdfjX

  2. dfYXdf

  3. YXrb=(XX)1S

  4. YY^=Xb

  5. E=YY^

  6. SSS=dfSSY^

  7. EXjCj=i=1nEiXij

  8. EC0i

    Ei[corrected]=Eij=1mCjXijnj=1mXij2

    (mianownik nie zmienia się w iteracjach, oblicz go wcześniej)

    E0 EC

    Ei[corrected]=Eij=1mCjXij3i=1nXij2j=1mXij2

    1

  9. SSEEi[corrected]=EiSSS/SSE

    mrSSSn

  10. CErYY[corrected]=Y^+E

  11. Y

  12. Yr

YrY


1YX

ttnphns
źródło
1
Dzięki za odpowiedź. To było rozwiązanie empiryczne / iteracyjne, o którym również myślałem. Do moich symulacji potrzebuję jednak bardziej analitycznego rozwiązania bez kosztownej procedury dopasowania. Na szczęście właśnie znalazłem rozwiązanie, które wkrótce opublikuję ...
Felix S
Działa to w przypadku generowania dwuwymiarowych normalnych, ale nie działa w przypadku arbitralnej dystrybucji (lub dowolnej innej „dodatkowej” dystrybucji)
Makro
1
Nie rozumiem, dlaczego proponujesz iterację, kiedy możesz bezpośrednio wytworzyć cały stożek rozwiązań. Czy to podejście ma jakiś szczególny cel?
whuber
1
Y
1
@ whuber, twój komentarz jest tym, na co czekałem; właściwie moja odpowiedź (o heteroscedastyczności, do której linkuję) była dla ciebie wyzwaniem: być może jest to zaproszenie do opublikowania twojego rozwiązania - tak dokładnego i błyskotliwego, jak zwykle.
ttnphns,
4

Miałem ochotę trochę programować, więc wziąłem usuniętą odpowiedź @ Adama i postanowiłem napisać fajną implementację w języku R. Skupiam się na używaniu stylu zorientowanego funkcjonalnie (tj. Pętli stylu lapply). Ogólna idea polega na pobraniu dwóch wektorów, losowym permutacji jednego z wektorów, dopóki nie zostanie osiągnięta pewna korelacja między nimi. To podejście jest bardzo brutalne, ale łatwe do wdrożenia.

Najpierw tworzymy funkcję, która losowo permutuje wektor wejściowy:

randomly_permute = function(vec) vec[sample.int(length(vec))]
randomly_permute(1:100)
  [1]  71  34   8  98   3  86  28  37   5  47  88  35  43 100  68  58  67  82
 [19]  13   9  61  10  94  29  81  63  14  48  76   6  78  91  74  69  18  12
 [37]   1  97  49  66  44  40  65  59  31  54  90  36  41  93  24  11  77  85
 [55]  32  79  84  15  89  45  53  22  17  16  92  55  83  42  96  72  21  95
 [73]  33  20  87  60  38   7   4  52  27   2  80  99  26  70  50  75  57  19
 [91]  73  62  23  25  64  51  30  46  56  39

... i utwórz przykładowe dane

vec1 = runif(100)
vec2 = runif(100)

... napisz funkcję, która permutuje wektor wejściowy i koreluje go z wektorem referencyjnym:

permute_and_correlate = function(vec, reference_vec) {
    perm_vec = randomly_permute(vec)
    cor_value = cor(perm_vec, reference_vec)
    return(list(vec = perm_vec, cor = cor_value))
  }
permute_and_correlate(vec2, vec1)
$vec
  [1] 0.79072381 0.23440845 0.35554970 0.95114398 0.77785348 0.74418811
  [7] 0.47871491 0.55981826 0.08801319 0.35698405 0.52140366 0.73996913
 [13] 0.67369873 0.85240338 0.57461506 0.14830718 0.40796732 0.67532970
 [19] 0.71901990 0.52031017 0.41357545 0.91780357 0.82437619 0.89799621
 [25] 0.07077250 0.12056045 0.46456652 0.21050067 0.30868672 0.55623242
 [31] 0.84776853 0.57217746 0.08626022 0.71740151 0.87959539 0.82931652
 [37] 0.93903143 0.74439384 0.25931398 0.99006038 0.08939812 0.69356590
 [43] 0.29254936 0.02674156 0.77182339 0.30047034 0.91790830 0.45862163
 [49] 0.27077191 0.74445997 0.34622648 0.58727094 0.92285322 0.83244284
 [55] 0.61397396 0.40616274 0.32203732 0.84003379 0.81109473 0.50573325
 [61] 0.86719899 0.45393971 0.19701975 0.63877904 0.11796154 0.26986325
 [67] 0.01581969 0.52571331 0.27087693 0.33821824 0.52590383 0.11261002
 [73] 0.89840404 0.82685046 0.83349287 0.46724807 0.15345334 0.60854785
 [79] 0.78854984 0.95770015 0.89193212 0.18885955 0.34303707 0.87332019
 [85] 0.08890968 0.22376395 0.02641979 0.43377516 0.58667068 0.22736077
 [91] 0.75948043 0.49734797 0.25235660 0.40125309 0.72147500 0.92423638
 [97] 0.27980561 0.71627101 0.07729027 0.05244047

$cor
[1] 0.1037542

... i iteruj tysiąc razy:

n_iterations = lapply(1:1000, function(x) permute_and_correlate(vec2, vec1))

Zauważ, że reguły określania zakresu R zapewniają vec1i vec2znajdują się w środowisku globalnym, poza anonimową funkcją używaną powyżej. Zatem permutacje są względne w stosunku do oryginalnych zestawów danych testowych, które wygenerowaliśmy.

Następnie znajdujemy maksymalną korelację:

cor_values = sapply(n_iterations, '[[', 'cor')
n_iterations[[which.max(cor_values)]]
$vec
  [1] 0.89799621 0.67532970 0.46456652 0.75948043 0.30868672 0.83244284
  [7] 0.86719899 0.55623242 0.63877904 0.73996913 0.71901990 0.85240338
 [13] 0.81109473 0.52571331 0.82931652 0.60854785 0.19701975 0.26986325
 [19] 0.58667068 0.52140366 0.40796732 0.22736077 0.74445997 0.40125309
 [25] 0.89193212 0.52031017 0.92285322 0.91790830 0.91780357 0.49734797
 [31] 0.07729027 0.11796154 0.69356590 0.95770015 0.74418811 0.43377516
 [37] 0.55981826 0.93903143 0.30047034 0.84776853 0.32203732 0.25235660
 [43] 0.79072381 0.58727094 0.99006038 0.01581969 0.41357545 0.52590383
 [49] 0.27980561 0.50573325 0.92423638 0.11261002 0.89840404 0.15345334
 [55] 0.61397396 0.27077191 0.12056045 0.45862163 0.18885955 0.77785348
 [61] 0.23440845 0.05244047 0.25931398 0.57217746 0.35554970 0.34622648
 [67] 0.21050067 0.08890968 0.84003379 0.95114398 0.83349287 0.82437619
 [73] 0.46724807 0.02641979 0.71740151 0.74439384 0.14830718 0.82685046
 [79] 0.33821824 0.71627101 0.77182339 0.72147500 0.08801319 0.08626022
 [85] 0.87332019 0.34303707 0.45393971 0.47871491 0.29254936 0.08939812
 [91] 0.35698405 0.67369873 0.27087693 0.78854984 0.87959539 0.22376395
 [97] 0.02674156 0.07077250 0.57461506 0.40616274

$cor
[1] 0.3166681

... lub znajdź wartość najbliższą korelacji 0,2:

n_iterations[[which.min(abs(cor_values - 0.2))]]
$vec
  [1] 0.02641979 0.49734797 0.32203732 0.95770015 0.82931652 0.52571331
  [7] 0.25931398 0.30047034 0.55981826 0.08801319 0.29254936 0.23440845
 [13] 0.12056045 0.89799621 0.57461506 0.99006038 0.27077191 0.08626022
 [19] 0.14830718 0.45393971 0.22376395 0.89840404 0.08890968 0.15345334
 [25] 0.87332019 0.92285322 0.50573325 0.40796732 0.91780357 0.57217746
 [31] 0.52590383 0.84003379 0.52031017 0.67532970 0.83244284 0.95114398
 [37] 0.81109473 0.35554970 0.92423638 0.83349287 0.34622648 0.18885955
 [43] 0.61397396 0.89193212 0.74445997 0.46724807 0.72147500 0.33821824
 [49] 0.71740151 0.75948043 0.52140366 0.69356590 0.41357545 0.21050067
 [55] 0.87959539 0.11796154 0.73996913 0.30868672 0.47871491 0.63877904
 [61] 0.22736077 0.40125309 0.02674156 0.26986325 0.43377516 0.07077250
 [67] 0.79072381 0.08939812 0.86719899 0.55623242 0.60854785 0.71627101
 [73] 0.40616274 0.35698405 0.67369873 0.82437619 0.27980561 0.77182339
 [79] 0.19701975 0.82685046 0.74418811 0.58667068 0.93903143 0.74439384
 [85] 0.46456652 0.85240338 0.34303707 0.45862163 0.91790830 0.84776853
 [91] 0.78854984 0.05244047 0.58727094 0.77785348 0.01581969 0.27087693
 [97] 0.07729027 0.71901990 0.25235660 0.11261002

$cor
[1] 0.2000199

Aby uzyskać wyższą korelację, musisz zwiększyć liczbę iteracji.

Paul Hiemstra
źródło
2

Y1Y2,,YnR

Rozwiązanie:

  1. CCT=R
  2. X2,,XnY1
  3. Y1
  4. Y=CXYiY1

Kod Python:

import numpy as np
import math
from scipy.linalg import toeplitz, cholesky
from statsmodels.stats.moment_helpers import cov2corr

# create the large correlation matrix R
p = 4
h = 2/p
v = np.linspace(1,-1+h,p)
R = cov2corr(toeplitz(v))

# create the first variable
T = 1000;
y = np.random.randn(T)

# generate p-1 correlated randoms
X = np.random.randn(T,p)
X[:,0] = y
C = cholesky(R)
Y = np.matmul(X,C)

# check that Y didn't change
print(np.max(np.abs(Y[:,0]-y)))

# check the correlation matrix
print(R)
print(np.corrcoef(np.transpose(Y)))

Wyjście testowe:

0.0
[[ 1.   0.5  0.  -0.5]
 [ 0.5  1.   0.5  0. ]
 [ 0.   0.5  1.   0.5]
 [-0.5  0.   0.5  1. ]]
[[ 1.          0.50261766  0.02553882 -0.46259665]
 [ 0.50261766  1.          0.51162821  0.05748082]
 [ 0.02553882  0.51162821  1.          0.51403266]
 [-0.46259665  0.05748082  0.51403266  1.        ]]
Aksakal
źródło
Y1
@ whuber to była literówka
Aksakal,
0

Wygeneruj zmienne normalne z podaną macierzą kowariancji SAMPLING

covsam <- function(nobs,covm, seed=1237) {; 
          library (expm);
          # nons=number of observations, covm = given covariance matrix ; 
          nvar <- ncol(covm); 
          tot <- nvar*nobs;
          dat <- matrix(rnorm(tot), ncol=nvar); 
          covmat <- cov(dat); 
          a2 <- sqrtm(solve(covmat)); 
          m2 <- sqrtm(covm);
          dat2 <- dat %*% a2 %*% m2 ; 
          rc <- cov(dat2);};
          cm <- matrix(c(1,0.5,0.1,0.5,1,0.5,0.1,0.5,1),ncol=3);
          cm; 
          res <- covsam(10,cm)  ;
          res;

Wygeneruj zmienne normalne z podaną macierzą kowariancji LUDNOŚCI

covpop <- function(nobs,covm, seed=1237) {; 
          library (expm); 
          # nons=number of observations, covm = given covariance matrix;
          nvar <- ncol(covm); 
          tot <- nvar*nobs;  
          dat <- matrix(rnorm(tot), ncol=nvar); 
          m2 <- sqrtm(covm);
          dat2 <- dat %*% m2;  
          rc <- cov(dat2); }; 
          cm <- matrix(c(1,0.5,0.1,0.5,1,0.5,0.1,0.5,1),ncol=3);
          cm; 
          res <- covpop(10,cm); 
          res
użytkownik3635627
źródło
2
Musisz nauczyć się formatować kod w odpowiedzi! Istnieje specjalna opcja oznaczania tekstu jako fragmentów kodu, użyj go!
kjetil b halvorsen
-6

Po prostu utwórz losowy wektor i sortuj, aż uzyskasz pożądane r.

Adam
źródło
W jakich sytuacjach byłoby to lepsze niż powyższe rozwiązania?
Andy W
Sytuacja, w której użytkownik chce prostej odpowiedzi. Przeczytałem podobne pytanie na forum r, i jego odpowiedź, która została udzielona.
Adam
3
r
3
Jeśli ta odpowiedź została podana na forum r-help, podejrzewam, że była albo (a) ironiczna (tj. Przeznaczona jako żart), albo (b) oferowana przez kogoś, kto nie jest zbyt wyrafinowany statystycznie. Mówiąc bardziej zwięźle, jest to słaba odpowiedź na pytanie. -1
gung