Cara menggunakan dekomposisi Cholesky, atau alternatif, untuk simulasi data berkorelasi

19

Saya menggunakan dekomposisi Cholesky untuk mensimulasikan variabel acak berkorelasi diberi matriks korelasi. Masalahnya, hasilnya tidak pernah mereproduksi struktur korelasi seperti yang diberikan. Berikut ini adalah contoh kecil dalam Python untuk menggambarkan situasi.

import numpy as np    

n_obs = 10000
means = [1, 2, 3]
sds = [1, 2, 3] # standard deviations 

# generating random independent variables 
observations = np.vstack([np.random.normal(loc=mean, scale=sd, size=n_obs)
                   for mean, sd in zip(means, sds)])  # observations, a row per variable

cor_matrix = np.array([[1.0, 0.6, 0.9],
                       [0.6, 1.0, 0.5],
                       [0.9, 0.5, 1.0]])

L = np.linalg.cholesky(cor_matrix)

print(np.corrcoef(L.dot(observations))) 

Ini mencetak:

[[ 1.          0.34450587  0.57515737]
 [ 0.34450587  1.          0.1488504 ]
 [ 0.57515737  0.1488504   1.        ]]

Seperti yang Anda lihat, matriks korelasi taksiran post-hoc berbeda secara drastis dari yang sebelumnya. Apakah ada bug dalam kode saya, atau ada beberapa alternatif untuk menggunakan dekomposisi Cholesky?

Edit

Maafkan saya atas kekacauan ini. Saya tidak berpikir ada kesalahan dalam kode dan / atau cara dekomposisi Cholesky diterapkan karena beberapa kesalahpahaman tentang materi yang saya pelajari sebelumnya. Sebenarnya saya yakin bahwa metode itu sendiri tidak dimaksudkan untuk menjadi tepat dan saya baik-baik saja dengan itu sampai situasi yang membuat saya memposting pertanyaan ini. Terima kasih telah menunjukkan kesalahpahaman yang saya miliki. Saya telah mengedit judul untuk lebih mencerminkan situasi sebenarnya seperti yang diusulkan oleh @Silverfish.

Eli Korvigo
sumber
1
Cholesky berfungsi dengan baik, dan ini benar-benar jenis pertanyaan "bisakah Anda menemukan bug dalam kode saya". Judul dan isi pertanyaan, seperti yang aslinya ditulis, pada dasarnya adalah "Cholesky tidak berfungsi, apa alternatifnya"? Itu akan sangat membingungkan bagi pengguna yang mencari situs ini. Haruskah pertanyaan ini diedit untuk mencerminkan hal ini? (Kelemahannya adalah bahwa jawaban javlacalle akan kurang relevan. Sisi baiknya adalah teks pertanyaan akan mencerminkan apa yang sebenarnya akan ditemukan oleh pencari di halaman.)
Silverfish
@Antoni Parellada Ya, saya pikir Anda telah menerjemahkan kode MATLAB saya untuk (a) cara yang benar untuk melakukannya ke Python numpy, lengkap dengan penyesuaian untuk np.linalg.cholesky menjadi segitiga bawah vs kol MATLAB menjadi segitiga atas. Saya sudah menerjemahkan kode OP yang salah ke MATLAB dan menggandakan hasil yang salah.
Mark L. Stone

Jawaban:

11

Pendekatan berdasarkan dekomposisi Cholesky harus bekerja, dijelaskan di sini dan ditunjukkan dalam jawaban oleh Mark L. Stone yang diposting hampir bersamaan dengan jawaban ini.

N(μ,Σ)

Y=QX+μ,denganQ=Λ1/2Φ,

YXΦΣΛΣΦ

Contoh dalam R(maaf saya tidak menggunakan perangkat lunak yang sama dengan yang Anda gunakan dalam pertanyaan):

n <- 10000
corM <- rbind(c(1.0, 0.6, 0.9), c(0.6, 1.0, 0.5), c(0.9, 0.5, 1.0))
set.seed(123)
SigmaEV <- eigen(corM)
eps <- rnorm(n * ncol(SigmaEV$vectors))
Meps <- matrix(eps, ncol = n, byrow = TRUE)    
Meps <- SigmaEV$vectors %*% diag(sqrt(SigmaEV$values)) %*% Meps
Meps <- t(Meps)
# target correlation matrix
corM
#      [,1] [,2] [,3]
# [1,]  1.0  0.6  0.9
# [2,]  0.6  1.0  0.5
# [3,]  0.9  0.5  1.0
# correlation matrix for simulated data
cor(Meps)
#           [,1]      [,2]      [,3]
# [1,] 1.0000000 0.6002078 0.8994329
# [2,] 0.6002078 1.0000000 0.5006346
# [3,] 0.8994329 0.5006346 1.0000000

Anda mungkin juga tertarik dengan posting ini dan posting ini .

javlacalle
sumber
Untuk membuat matriks korelasi yang direproduksi tepat, seseorang harus menghilangkan korelasi palsu dalam data acak dari generator acak sebelum menerapkannya pada prosedur pembuatan data. Misalnya, periksa korelasi data acak Anda dalam eps untuk melihat korelasi palsu itu terlebih dahulu.
Gottfried Helms
17

Orang-orang kemungkinan akan menemukan kesalahan Anda lebih cepat jika Anda menjelaskan apa yang Anda lakukan dengan kata-kata dan aljabar daripada kode (atau setidaknya menulisnya menggunakan pseudocode).

Anda tampaknya melakukan hal yang setara dengan ini (meskipun mungkin diubah):

  1. n×kZ

  2. σsayaμsaya

  3. Y=L.X

L.

Yang harus Anda lakukan adalah ini:

  1. n×kZ

  2. X=L.Z

  3. σsayaμsaya

Ada banyak penjelasan tentang algoritma ini di situs. misalnya

Bagaimana cara menghasilkan angka acak berkorelasi (diberikan berarti, varians dan tingkat korelasi)?

Bisakah saya menggunakan metode Cholesky untuk menghasilkan variabel acak berkorelasi dengan mean yang diberikan?

Yang ini membahasnya langsung dalam hal matriks kovarians yang diinginkan, dan juga memberikan algoritma untuk mendapatkan sampel kovarians yang diinginkan :

Menghasilkan data dengan matriks kovarians sampel yang diberikan

Glen_b -Reinstate Monica
sumber
11

Tidak ada yang salah dengan faktorisasi Cholesky. Ada kesalahan dalam kode Anda. Lihat edit di bawah.

Berikut adalah kode dan hasil MATLAB, pertama untuk n_obs = 10.000 seperti yang Anda miliki, kemudian untuk n_obs = 1e8. Untuk kesederhanaan, karena itu tidak mempengaruhi hasil, saya tidak peduli dengan cara, yaitu, saya membuatnya nol. Perhatikan bahwa chol MATLAB menghasilkan faktor Cholesky segitiga atas R dari matriks M sehingga R '* R = M. numpy.linalg.cholesky menghasilkan faktor Cholesky segitiga yang lebih rendah, sehingga diperlukan penyesuaian vs kode saya; tapi saya percaya kode Anda baik-baik saja dalam hal itu.

   >> correlation_matrix = [1.0, 0.6, 0.9; 0.6, 1.0, 0.5;0.9, 0.5, 1.0];
   >> SD = diag([1 2 3]);
   >> covariance_matrix = SD*correlation_matrix*SD
   covariance_matrix =
      1.000000000000000   1.200000000000000   2.700000000000000
      1.200000000000000   4.000000000000000   3.000000000000000
      2.700000000000000   3.000000000000000   9.000000000000000
   >> n_obs = 10000;
   >> Random_sample = randn(n_obs,3)*chol(covariance_matrix);
   >> disp(corr(Random_sample))
      1.000000000000000   0.599105015695768   0.898395949647890
      0.599105015695768   1.000000000000000   0.495147514173305
      0.898395949647890   0.495147514173305   1.000000000000000
   >> n_obs = 1e8;
   >> Random_sample = randn(n_obs,3)*chol(covariance_matrix);
   >> disp(corr(Random_sample))
      1.000000000000000   0.600101477583914   0.899986072541418
      0.600101477583914   1.000000000000000   0.500112824962378
      0.899986072541418   0.500112824962378   1.000000000000000

Sunting: Saya menemukan kesalahan Anda. Anda salah menerapkan standar deviasi. Ini setara dengan apa yang Anda lakukan, yang salah.

   >> n_obs = 10000;
   >> Random_sample = randn(n_obs,3)*SD*chol(correlation_matrix);
   >> disp(corr(Random_sample))
      1.000000000000000   0.336292731308138   0.562331469857830
      0.336292731308138   1.000000000000000   0.131270077244625
      0.562331469857830   0.131270077244625   1.000000000000000
   >> n_obs=1e8;
   >> Random_sample = randn(n_obs,3)*SD*chol(correlation_matrix);
   >> disp(corr(Random_sample))
      1.000000000000000   0.351254525742470   0.568291702131030
      0.351254525742470   1.000000000000000   0.140443281045496
      0.568291702131030   0.140443281045496   1.000000000000000
Mark L. Stone
sumber
6

CV bukan tentang kode, tetapi saya tertarik untuk melihat bagaimana ini akan terlihat setelah semua jawaban yang baik, dan khususnya kontribusi @Mark L. Stone. Jawaban aktual untuk pertanyaan diberikan pada posnya (mohon kredit posnya jika ragu). Saya sedang memindahkan info tambahan ini di sini untuk memudahkan pengambilan posting ini di masa mendatang. Tanpa mengecilkan salah satu jawaban luar biasa lainnya, setelah jawaban Markus, ini menyelesaikan masalah dengan memperbaiki pos di OP.

Sumber

DI PYTHON:

import numpy as np

no_obs = 1000             # Number of observations per column
means = [1, 2, 3]         # Mean values of each column
no_cols = 3               # Number of columns

sds = [1, 2, 3]           # SD of each column
sd = np.diag(sds)         # SD in a diagonal matrix for later operations

observations = np.random.normal(0, 1, (no_cols, no_obs)) # Rd draws N(0,1) in [3 x 1,000]

cor_matrix = np.array([[1.0, 0.6, 0.9],
                       [0.6, 1.0, 0.5],
                       [0.9, 0.5, 1.0]])          # The correlation matrix [3 x 3]

cov_matrix = np.dot(sd, np.dot(cor_matrix, sd))   # The covariance matrix

Chol = np.linalg.cholesky(cov_matrix)             # Cholesky decomposition

array([[ 1.        ,  0.        ,  0.        ],
       [ 1.2       ,  1.6       ,  0.        ],
       [ 2.7       , -0.15      ,  1.29903811]])

sam_eq_mean = Chol .dot(observations)             # Generating random MVN (0, cov_matrix)

s = sam_eq_mean.transpose() + means               # Adding the means column wise
samples = s.transpose()                           # Transposing back

print(np.corrcoef(samples))                       # Checking correlation consistency.

[[ 1.          0.59167434  0.90182308]
 [ 0.59167434  1.          0.49279316]
 [ 0.90182308  0.49279316  1.        ]]

DALAM [R]:

no_obs = 1000             # Number of observations per column
means = 1:3               # Mean values of each column
no_cols = 3               # Number of columns

sds = 1:3                 # SD of each column
sd = diag(sds)         # SD in a diagonal matrix for later operations

observations = matrix(rnorm(no_cols * no_obs), nrow = no_cols) # Rd draws N(0,1)

cor_matrix = matrix(c(1.0, 0.6, 0.9,
                      0.6, 1.0, 0.5,
                      0.9, 0.5, 1.0), byrow = T, nrow = 3)     # cor matrix [3 x 3]

cov_matrix = sd %*% cor_matrix %*% sd                          # The covariance matrix

Chol = chol(cov_matrix)                                        # Cholesky decomposition

     [,1] [,2]      [,3]
[1,]    1  1.2  2.700000
[2,]    0  1.6 -0.150000
[3,]    0  0.0  1.299038

sam_eq_mean = t(observations) %*% Chol          # Generating random MVN (0, cov_matrix)

samples = t(sam_eq_mean) + means

cor(t(samples))

          [,1]      [,2]      [,3]
[1,] 1.0000000 0.6071067 0.8857339
[2,] 0.6071067 1.0000000 0.4655579
[3,] 0.8857339 0.4655579 1.0000000

colMeans(t(samples))
[1] 1.035056 2.099352 3.065797
apply(t(samples), 2, sd)
[1] 0.9543873 1.9788250 2.8903964
Antoni Parellada
sumber
1

Seperti yang telah ditunjukkan orang lain: karya cholesky. Di sini sepotong kode yang sangat pendek dan sangat dekat dengan pseudocode: sebuah codepiece di MatMate:

Co = {{1.0, 0.6, 0.9},  _
      {0.6, 1.0, 0.5},  _
      {0.9, 0.5, 1.0}}           // make correlation matrix


chol = cholesky(co)              // do cholesky-decomposition           
data = chol * unkorrzl(randomn(3,100,0,1))  
                                 // dot-multiply cholesky with random-
                                 // vectors with mean=0, sdev=1  
                                 //(refined by a "decorrelation" 
                                 //to remove spurious/random correlations)   


chk = data *' /100               // check the correlation of the data
list chk

1.0000  0.6000  0.9000
0.6000  1.0000  0.5000
0.9000  0.5000  1.0000
Gottfried Helms
sumber