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

Question

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

matlabrandommultidimensional-arraynormal-distribution

3

这个问题涉及到多维正态分布中协方差矩阵的使用：

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

x = mu + chol(Sigma) * Z

或者

x = mu + Sigma ^ 0.5 * Z

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

- Futurist

这是第一个选项，但我认为你需要包括一个复杂的转置。请参见此处。或者，使用Matlab的mvnrnd函数，它接受sigma并在内部执行Chlesky分解。 - Luis Mendo

@LuisMendo，据我看来 ^0.5 或 sqrtm 也可以工作。那么 chol 的优点是什么？我猜是恢复 Sigma 的数值精度更高？ - A. Donda

1

@A.Donda 抱歉，我误读了 Sigma^0.5。实际上我不确定。我知道的是第一个选项是正确的。也许第二个选项也是。 - Luis Mendo

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- A. Donda · Accepted Answer

如果你指的是多元正态分布的密度：density，那么它包含了Σ的逆矩阵和行列式的标量平方根，但不包括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是一个上三角矩阵。

然而，我们也可以寻找一个共轭解A' = A，那么A' * A就成为了矩阵平方A^2。这个问题的解决方法是通过计算矩阵平方根来替换Sigma的每个特征值为其平方根（或其相反数）；一般情况下，对于n个正特征值，可能有2ⁿ个解。Matlab函数sqrtm返回主矩阵平方根，即唯一的非负定解。因此，

A = sqrtm(Sigma)

这也适用于IT技术相关内容。 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 对奇异的输入会返回警告，但会返回有效结果）。