如何从标准正态值生成多元正态分布？

Question

我需要仅使用随机值生成器而不使用 scipy 或 numpy 生成器来生成多元正态分布。

我需要生成以下内容

这是我的尝试

V = np.array([
    [1, 2], 
    [2, 5]])
B = np.linalg.cholesky(V)
A = np.array([1,2])

# norm() return one number from standard normal distribution 
n1 = np.array([norm() for _ in range(40)])
n2 = np.array([norm() for _ in range(40)])
np.array([n1,n2]).T.dot(B) + A

这里，我使用 Cholesky 分解为in this post

但是，我认为这是不正确的。

Answer 1

您的代码几乎是正确的，但如果您应用numpy's cov 函数，您可以检查您的数字是否不具有所需的协方差属性：

res = np.array([n1,n2]).T.dot(B) + A
np.cov(res.T).round()
# returns ~ 
# array([[5., 2.],
#        [2., 1.]])

请注意，与所需值相比，元素 1,1 和 2,2 是交换的。

要利用 numpy 的 CPU 向量化矩阵乘法，您可以使用 numpy's dot 函数。您将 N 个二维输入向量 Z 正确排列为 Nx2 维向量 np.array([n1,n2]).T。但是正如您在Cholesky decomposition and variance 问题中指出的那样，Z 值必须从左侧乘以B，您还想将其合并到dot 函数的广播规则中，问题就在这里。代码np.array([n1,n2]).T.dot(B) 将（数组）Z 从右侧相乘，而不是从左侧相乘。要计算 B 的左积，您需要使用 dot(B.T)

这个例子也表明协方差矩阵有正确的形式

import random
import numpy as np

random.seed(0)
N=10000

V = np.array([
    [1, 2],
    [2, 5]])
B = np.linalg.cholesky(V)
A = np.array([1, 2])

# norm() return one number from standard normal distribution
n1 = np.array([random.gauss(0, 1) for _ in range(10000)])
n2 = np.array([random.gauss(0, 1) for _ in range(10000)])
res = np.array([n1, n2]).T.dot(B.T) + A

np.cov(res.T).round()
# returns ~  array([[1., 2.],
#                   [2., 5.]])

如图。在随机点下方绘制，以及协方差矩阵的特征向量及其特征值的平方根长度，如Wikipedia。