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

Question

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

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

x = mu + chol(Sigma) * Z

或

x = mu + Sigma ^ 0.5 * Z

?

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

Answer 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函数sqrtmreturns求主矩阵平方根，也就是nonnegative-definite的唯一解。因此，

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 分解确实有缺点，它仅针对 positive-definite Σ 定义，而矩阵平方根的要求仅仅是 Σ 是 nonnegative-definite（sqrtm returns 单个输入的警告，但 returns 有效结果）。

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

Generating multivariate normally distributed random numbers in Matlab

matlab

random

multidimensional-array

normal-distribution