使用Python numpy.random.randn绘制多元高斯分布样本

Question

我正在研究高斯混合模型，并遇到了这段代码，该代码从2个双变量高斯分布中抽取了许多样本。我不理解的是代码中使用的技术：

import numpy as np

# Number of samples per component
n_samples = 500

# Generate random sample, two components
np.random.seed(0)
C = np.array([[0., -0.1], [1.7, .4]])
X = np.r_[np.dot(np.random.randn(n_samples, 2), C),
          .7 * np.random.randn(n_samples, 2) + np.array([-6, 3])]

（原始链接：http : //scikit-learn.org/stable/auto_examples/mixture/plot_gmm_selection.html#sphx-glr-auto-examples-mixture-plot-gmm-selection-py）

根据此Wikipedia链接，我们可以通过Cholesky分解协方差矩阵来生成多元高斯样本，然后将其乘以由标准正态分布的成分组成的向量。

我的问题是代码中的C变量不是下三角矩阵，那么它在多元高斯随机生成中如何有意义？

Answer 1

X是两个双变量正态分布的混合。一半的样品与计算np.dot(np.random.randn(n_samples, 2), C)，其中C = np.array([[0., -0.1], [1.7, .4]])。此分布等效于其协方差为的分布C.T.dot(C)。也就是说，您可以使用来从相同的分布生成样本np.random.multivariate_normal([0, 0], C.T.dot(C), n_samples)。

请参阅我之前写的这些注释：“相关随机样本”。（在这些注释中，3x3矩阵C在右侧乘以一个形状样本(3, num_samples)。换句话说，这些注释使用此处使用的公式的转置，因此协方差矩阵为C.dot(C.T)。） C不必是较低的三角形。但是通常情况下，您会得到一个协方差矩阵，并且想找到C。如果使用Cholesky分解来查找C，则通过构造将其变为较低的三角形。

这个ipython会话演示了这两种方法从相同的分布生成样本：

In [60]: C = np.array([[0., -0.1], [1.7, .4]])

In [61]: X1 = np.dot(np.random.randn(n_samples, 2), C)

In [62]: X2 = np.random.multivariate_normal([0, 0], C.T.dot(C), n_samples)

In [63]: plot(X1[:,0], X1[:,1], 'g*', alpha=0.2)
Out[63]: [<matplotlib.lines.Line2D at 0x113c17550>]

In [64]: plot(X2[:,0], X2[:,1], 'ko', alpha=0.2, ms=4)
Out[64]: [<matplotlib.lines.Line2D at 0x113c3ba58>]