如何使用 Cholesky 分解或替代方法進行相關數據模擬

我使用 Cholesky 分解來模擬給定相關矩陣的相關隨機變量。問題是，結果永遠不會重現給出的相關結構。這是 Python 中的一個小示例來說明這種情況。

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))) 

這打印：

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

如您所見，事後估計的相關矩陣與之前的有很大不同。我的代碼中是否存在錯誤，或者是否有使用 Cholesky 分解的替代方法？

編輯

我請你原諒這個爛攤子。由於對我之前研究過的材料的一些誤解，我認為代碼和/或 Cholesky 分解的應用方式沒有錯誤。事實上，我確信該方法本身並不精確，並且在讓我發布這個問題之前，我一直對此表示滿意。感謝您指出我的誤解。我已經編輯了標題以更好地反映@Silverfish 提出的真實情況。

基於 Cholesky 分解的方法應該有效，它在此處進行了描述， 並在 Mark L. Stone 的答案中顯示，幾乎與該答案同時發布。

儘管如此，我有時會從多元正態分佈中生成繪圖 如下：

在哪裡是最後的抽籤，是從單變量標準正態分佈中抽取的，是一個包含目標矩陣的歸一化特徵向量的矩陣和是一個包含特徵值的對角矩陣以與列中的特徵向量相同的順序排列.

示例R（對不起，我沒有使用您在問題中使用的相同軟件）：

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

---

你可能也對 這篇文章 和這篇文章感興趣。

引用自：https://stats.stackexchange.com/questions/160054