Jak mogę generować dane z uprzednio określoną macierzą korelacji?

19

Próbuję wygenerować skorelowaną losową sekwencję ze średnią = , wariancja = , współczynnik korelacji = . W poniższym kodzie używam & jako standardowych odchyleń i & jako środków.010.8s1s2m1m2

p = 0.8 
u = randn(1, n)
v = randn(1, n)
x = s1 * u + m1
y = s2 * (p * u + sqrt(1 - p^2) * v) + m2

To daje mi poprawkę corrcoef()0,8 pomiędzy xi y. Moje pytanie brzmi: w jaki sposób mogę wygenerować szereg oznacza, że ​​jeśli chcę z, to jest również skorelowane y(z tą samą korelacją ), ale nie z . Czy jest jakiś szczególny wzór, który muszę znać? Znalazłem jeden, ale nie mogłem go zrozumieć.r=0.8x

anisa
źródło

Odpowiedzi:

21

Wygląda na to, że pytasz, jak wygenerować dane z określoną macierzą korelacji.

Przydatnym faktem jest to, że jeśli masz losowy wektor z macierzą kowariancji , to losowy wektor ma średnią i macierz kowariancji . Tak więc, jeśli zaczniesz od danych oznaczających zero, pomnożenie przez tego nie zmieni, więc twoje pierwsze wymaganie jest łatwo spełnione. xΣAxAE(x)Ω=AΣATA

Powiedzmy zacząć (średnia zero) danych nieskorelowanych (tj macierz kowariancji jest przekątna) - skoro mówimy o macierzy korelacji, po prostu wziąć . Możesz przekształcić to w dane za pomocą danej macierzy kowariancji, wybierając jako cholesky pierwiastek kwadratowy z - wtedy miałby pożądaną macierz kowariancji .Σ=IAΩAxΩ

W twoim przykładzie wydaje się, że chcesz czegoś takiego:

Ω=(1.80.81.80.81)

Niestety, macierz ta nie jest jednoznacznie określona, ​​więc nie może być macierzą kowariancji - można to sprawdzić, widząc, że wyznacznik jest ujemny. Być może zamiast tego

Ω=(1.8.3.81.8.3.81)    or   Ω=(12/302/312/302/31)

wystarczyłoby. Nie jestem pewien, jak obliczyć cholesky pierwiastek kwadratowy w Matlabie (który wydaje się być tym, którego używasz), ale Rmożesz użyć tej chol()funkcji.

W tym przykładzie dla dwóch wymienionych wyżej odpowiednimi mnożnikami macierzy byłyby odpowiednioΩ

A=(100.8.60.3.933.1972)    or   A=(1002/3.745300.8944.4472)

Do tego Rcelu użyto kodu:

x = matrix(0,3,3)
x[1,]=c(1,.8,.3)
x[2,]=c(.8,1,.8)
x[3,]=c(.3,.8,1)
t(chol(x))

     [,1]      [,2]      [,3]
[1,]  1.0 0.0000000 0.0000000
[2,]  0.8 0.6000000 0.0000000
[3,]  0.3 0.9333333 0.1972027

x[1,]=c(1,2/3,0)
x[2,]=c(2/3,1,2/3)
x[3,]=c(0,2/3,1)
t(chol(x))

      [,1]      [,2]      [,3]
[1,] 1.0000000 0.0000000 0.0000000
[2,] 0.6666667 0.7453560 0.0000000
[3,] 0.0000000 0.8944272 0.4472136
Makro
źródło
1
Funkcja MATLAB jest również nazywany chol. Zauważ, że może to być dość niestabilna liczbowo, jeśli jest prawie pojedyncza. W takim przypadku zastosowanie uzyskanego symetrycznego pierwiastka kwadratowego, np. Za pomocą SVD, jest często lepszym wyborem pod względem stabilności numerycznej. :)Ω
kardynał
1
Oczywiście to prawda @ cardinal - wiele teoretycznie usprawiedliwionych rzeczy psuje się, gdy próbujesz robić rzeczy numerycznie z prawie pojedynczymi macierzami. (Wygodnie) wyobrażałem sobie sytuację, w której docelowa matryca korelacji nie znajdowała się w dziedzinie, w której był to problem. Dobrze, że zwróciłeś na to uwagę - dzięki (i dzięki za edycję mojej drugiej odpowiedzi)
Macro
1
Głównym powodem, dla którego o tym myślałem, było twoje bystre oko w rozpoznaniu, że pierwsza sugestia PO nie była nawet pozytywnie określona. I miejmy nadzieję, że redakcja drugiego pytania nie była nadgorliwa; Lubię obie te odpowiedzi.
kardynał
7

Jeśli używasz R, możesz również użyć funkcji mvrnorm z pakietu MASS, zakładając, że chcesz normalnie rozproszonych zmiennych. Implementacja jest podobna do powyższego opisu Makra, ale wykorzystuje wektory własne macierzy korelacji zamiast chłodnego rozkładu i skalowania z rozkładem pojedynczej wartości (jeśli opcja empiryczna jest ustawiona na true).

XΣγλΣ

X=γλXT.

W przypadku, gdy X”oznacza zwykle rozmieszczone macierz macierzy korelacji i środków kolumn jest taka sama jak .ΣX

Zauważ, że macierz korelacji musi być definitywnie dodatnia, ale przydatne będzie przekonwertowanie jej za pomocą funkcji nearPD z pakietu Matrix w R.

zzk
źródło
1

ΣyxΣx=IΣyΛV

Σy=V.ΛV.T.=(V.Λ)(ΛT.V.T.)=ZAZAT.

y=ZAx

Mario Sansone
źródło