Bagaimana saya bisa menghasilkan data dengan matriks korelasi yang ditentukan sebelumnya?

19

Saya mencoba untuk menghasilkan urutan acak berkorelasi dengan mean = 0 , varians = 1 , koefisien korelasi = 0.8 . Dalam kode di bawah ini, saya menggunakan s1& s2sebagai standar deviasi, dan m1& m2sebagai sarana.

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

Ini memberi saya benar corrcoef()0,8 antara xdan y. Pertanyaan saya adalah bagaimana saya bisa menghasilkan seri artinya jika saya ingin zitu juga berkorelasi dengan y(dengan korelasi yang sama r=0.8 ), tetapi tidak dengan x. Apakah ada formula khusus yang perlu saya ketahui? Saya menemukan satu tetapi tidak bisa memahaminya.

anisa
sumber

Jawaban:

21

Tampaknya Anda bertanya bagaimana cara menghasilkan data dengan matriks korelasi tertentu.

Sebuah fakta yang berguna adalah bahwa jika Anda memiliki vektor acak dengan kovarian matriks Σ , maka random vektor A x memiliki rata-rata A E ( x ) dan kovariansi matriks Ω = A Σ A T . Jadi, jika Anda mulai dengan data yang berarti nol, mengalikan dengan A tidak akan mengubahnya, sehingga persyaratan pertama Anda mudah dipenuhi. xΣAxAE(x)Ω=AΣATA

Katakanlah Anda mulai dengan (mean nol) Data berkorelasi (yaitu matriks kovarians diagonal) - karena kita sedang berbicara tentang matriks korelasi, mari kita hanya mengambil . Anda dapat mengubah ini menjadi data dengan matriks kovarians yang diberikan dengan memilih A menjadi akar kuadrat cholesky dari Ω - maka A x akan memiliki matriks kovarians yang diinginkan Ω .Σ=IAΩAxΩ

Dalam contoh Anda, Anda tampaknya menginginkan sesuatu seperti ini:

Ω=(1.80.81.80.81)

Sayangnya matriks itu tidak pasti positif, jadi tidak bisa menjadi matriks kovarians - Anda dapat memeriksanya dengan melihat bahwa determinannya negatif. Mungkin sebaliknya

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

sudah cukup. Saya tidak yakin bagaimana menghitung root kuadrat cholesky di matlab (yang tampaknya adalah apa yang Anda gunakan) tetapi RAnda dapat menggunakan chol()fungsinya.

Dalam contoh ini, selama dua s tercantum di atas kelipatan matriks yang tepat (masing-masing) akanΩ

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

The Rkode yang digunakan untuk sampai pada ini:

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
sumber
1
Fungsi MATLAB juga disebut chol. Perhatikan bahwa ini bisa sangat tidak stabil secara numerik jika hampir tunggal. Dalam hal itu, menggunakan akar kuadrat simetris yang diperoleh, misalnya, melalui SVD, seringkali merupakan pilihan yang lebih baik dalam hal stabilitas numerik. :)Ω
kardinal
1
Tentu saja itu benar @ kardinal - banyak hal yang secara teoritis dibenarkan menjadi buruk ketika Anda mencoba melakukan hal-hal secara numerik dengan matriks yang hampir tunggal. Saya (dengan mudah) membayangkan situasi di mana matriks korelasi target tidak ada di ranah di mana ini menjadi masalah. Ada baiknya Anda menunjukkan ini - terima kasih (dan terima kasih atas hasil edit untuk jawaban saya yang lain)
Makro
1
Alasan utama saya berpikir tentang ini adalah karena mata Anda yang tajam dalam mengakui bahwa saran pertama OP bahkan tidak pasti positif. Dan, semoga hasil edit untuk pertanyaan lain tidak terlalu bersemangat; Saya suka kedua jawaban ini.
kardinal
7

Jika Anda menggunakan R, Anda juga dapat menggunakan fungsi mvrnorm dari paket MASS, dengan asumsi Anda ingin variabel yang terdistribusi normal. Implementasinya mirip dengan deskripsi Makro di atas, tetapi menggunakan vektor eigen dari matriks korelasi alih-alih dekomposisi cholesky dan penskalaan dengan dekomposisi nilai singular (jika opsi empiris disetel ke true).

Jika XΣγλΣ

X=γλXT

ΣX

Perhatikan bahwa matriks korelasi harus pasti positif, tetapi mengonversinya dengan fungsi nearPD dari paket Matrix di R akan bermanfaat.

zzk
sumber
1

Solusi alternatif tanpa faktorisasi cholesky adalah sebagai berikut. Biarkan kovarians matriks yang diinginkan dan anggaplah Anda memiliki data x dengan Σ x = I . Misalkan ΣΣyxΣx=IΣyΛV

Anda dapat menulis Σy=VΛVT=(VΛ)(ΛTVT)=AAT

y=Ax

Mario Sansone
sumber