在 Matlab 中生成多元正态分布随机数答案

【问题标题】：Generating multivariate normally distributed random numbers in Matlab在 Matlab 中生成多元正态分布随机数
【发布时间】：2015-03-22 08:08:39
【问题描述】：

本题是关于协方差矩阵在多维正态分布中的使用：

我想在 Matlab 中生成具有给定均值 mu 和协方差矩阵 Sigma 的多维随机数 x。假设Z 是一个标准的正态分布随机数（例如使用randn 生成），那么正确的代码是什么：

x = mu + chol(Sigma) * Z

或

x = mu + Sigma ^ 0.5 * Z

我不确定在多维正态分布的定义中使用协方差矩阵——分母中的行列式是平方根还是 Cholesky 因子...

【问题讨论】：

这是第一个选项，但我认为您需要包含一个复杂的转置。见here。或者，使用 Matlab 的 mvnrnd 函数，该函数接受 sigma 并在内部进行 Chlesky 分解
@LuisMendo, afaics ^0.5 或 sqrtm 应该也可以。 chol的优势是什么？我猜恢复 Sigma 的数值精度？
@A.Donda 对不起，我误读了Sigma^0.5。其实我不确定。我所知道的是第一个选项是正确的。也许第二个太

标签： matlab random multidimensional-array normal-distribution

【解决方案1】：

如果根据定义，您指的是多元正态分布的密度：

它既不包含 Cholesky 分解也不包含 Σ 的矩阵平方根，而是它的逆矩阵和行列式的标量平方根。

但是对于从这个分布中生成随机数，密度没有帮助。它甚至不是多元正态分布的最一般描述，因为密度公式仅对正定矩阵 Σ 有意义，而如果特征值为零，则分布也被定义——这仅意味着方向上的方差为 0各自的特征向量。

您的问题遵循从randn 产生的标准多元正态分布随机数Z 开始，然后应用线性变换的方法。假设mu 是一个p 维的行向量，我们需要一个nxp 维的随机矩阵（每行一个观察值，每列一个变量）：

Z = randn(n, p);
x = mu + Z * A;

我们需要一个矩阵A，使得x 的协方差为Sigma。由于Z 的协方差是单位矩阵，所以x 的协方差由A' * A 给出。 Cholesky decomposition 给出了一个解决方案，所以自然的选择是

A = chol(Sigma);

其中A 是一个上三角矩阵。

但是，我们也可以搜索 Hermitian 解，A' = A，然后 A' * A 变为 A^2，即矩阵平方。对此的解决方案由matrix square root 给出，其计算方法是将Sigma 的每个特征值替换为其平方根（或其负数）；一般来说，n 个正特征值有 2ⁿ 种可能的解。 Matlab 函数sqrtm 返回主矩阵平方根，这是唯一的非负定解。因此，

A = sqrtm(Sigma)

也可以。 A ^ 0.5 原则上应该这样做。

使用此代码进行模拟

p = 10;
n = 1000;

nr = 1000;
cp = nan(nr, 1);
sp = nan(nr, 1);
pp = nan(nr, 1);

for i = 1 : nr
    x = randn(n, p);
    Sigma = cov(x);

    cS = chol(Sigma);
    cp(i) = norm(cS' * cS - Sigma);

    sS = sqrtm(Sigma);
    sp(i) = norm(sS' * sS - Sigma);

    pS = Sigma ^ 0.5;
    pp(i) = norm(pS' * pS - Sigma);
end    

mean([cp sp pp])

chol 比其他两种方法更精确，并且分析表明它也快得多，对于 p = 10 和 p = 100。

然而，Cholesky 分解确实有一个缺点，它只为正定 Σ 定义，而矩阵平方根的要求仅仅是 Σ 是非负定的（sqrtm 会为奇异输入返回警告，但返回一个有效的结果）。

【讨论】：

干得好！我不知道sqrtm
很好的清晰全面的解释。非常感谢！
如果Sigma 有负值并且不是正定值怎么办？乔尔斯基会哭。有什么建议吗？
@makis，Sigma 可以有负条目（非对角线），只要它没有负特征值。如果是，则它不是一个有效的协方差矩阵，即不存在它是协方差的可能分布。你应该问自己的问题是为什么会有负特征值。