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

Question

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

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

x = mu + chol(Sigma) * Z

或者

x = mu + Sigma ^ 0.5 * Z

?

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

Answer 1

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

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

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

你的问题如下，从标准的多元正态分布的随机数启动方式Z，通过产生的randn，然后应用线性变换。假设这mu是一个一p维行向量，我们想要一个nxp维随机矩阵（每行一个观察，每列一个变量）：

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

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

A = chol(Sigma);

其中A是上三角矩阵。

但是，我们也可以搜索 Hermitian 解，A' = A，然后A' * A变成A^2，矩阵平方。对此的解决方案由矩阵平方根给出，该矩阵通过将 的每个特征值替换Sigma为其平方根（或其负数）来计算；一般有2个？n 个正特征值的可能解。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对于单一输入返回警告，但返回有效结果）。