机器算法验证 - 如何使用 Cholesky 分解或替代方法进行相关数据模拟 - 吾爱随笔录

如何使用 Cholesky 分解或替代方法进行相关数据模拟

机器算法验证相关性随机生成胆汁分解

2022-01-26 21:42:52

我使用 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 提出的真实情况。

4个回答

如果您用文字和代数而不是代码（或至少使用伪代码编写）来解释您所做的事情，人们可能会更快地发现您的错误。

您似乎正在做与此等效的操作（尽管可能已转置）：

生成一个 $n\times k$ 标准法线矩阵， $Z$
将列乘以 $\sigma_i$ 并添加 $\mu_i$ 获得非标准法线
计算 $Y=LX$ 得到相关的法线。

在哪里 $L$ 是相关矩阵的左 Cholesky 因子。

你应该做的是：

生成一个 $n\times k$ 标准法线矩阵， $Z$
计算 $X=LZ$ 得到相关的法线。
将列乘以 $\sigma_i$ 并添加 $\mu_i$ 获得非标准法线

现场有很多关于这个算法的解释。例如

如何生成相关随机数（给定均值、方差和相关程度）？

我可以使用 Cholesky 方法生成具有给定均值的相关随机变量吗？

这个直接根据期望的协方差矩阵讨论它，并且还给出了获得期望的样本协方差的算法：

使用给定的样本协方差矩阵生成数据

基于 Cholesky 分解的方法应该有效，它在此处进行了描述，并在 Mark L. Stone 的答案中显示，几乎与该答案同时发布。

尽管如此，我有时会从多元正态分布中生成绘图 $N(\vec\mu, \Sigma)$ 如下：

Y = Q X + \vec{μ}, with Q = Λ^{1 / 2} Φ,

$Y = Q X + \vec\mu \,, \quad \hbox{with}\quad Q=\Lambda^{1/2}\Phi \,,$

在哪里 $Y$ 是最后的抽签， $X$ 是从单变量标准正态分布中抽取的， $\Phi$ 是一个包含目标矩阵的归一化特征向量的矩阵 $\Sigma$ 和 $\Lambda$ 是一个包含特征值的对角矩阵 $\Sigma$ 以与列中的特征向量相同的顺序排列 $\Phi$ .

示例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

你可能也对这篇文章和这篇文章感兴趣。

Cholesky 分解没有任何问题。您的代码中有错误。请参阅下面的编辑。

这是 MATLAB 代码和结果，首先是 n_obs = 10000，然后是 n_obs = 1e8。为简单起见，因为它不影响结果，所以我不关心手段，即我将它们设为零。请注意，MATLAB 的 chol 产生矩阵 M 的上三角 Cholesky 因子 R，使得 R' * R = M。numpy.linalg.cholesky 产生下三角 Cholesky 因子，因此需要对我的代码进行调整；但我相信你的代码在这方面很好。

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

编辑：我发现你的错误。您错误地应用了标准偏差。这相当于你所做的，这是错误的。

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

CV 与代码无关，但我很想知道在所有好的答案之后，这将如何看待，特别是 @Mark L. Stone 的贡献。该问题的实际答案在他的帖子中提供（如有疑问，请记下他的帖子）。我将这里附加的信息移动到这里，以方便将来检索这篇文章。在马克的回答之后，在不淡化任何其他优秀答案的情况下，通过更正 OP 中的帖子来解决问题。

来源

在 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.        ]]

在 [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

其它你可能感兴趣的问题

上一篇贝叶斯统计的温和方法下一篇为什么逻辑回归在高维度上特别容易过拟合？