正如@hpaulj建议的那样，可以从标准多元正态分布生成样本，然后使用einsum和/或broadcasting来转换样本。缩放是通过将标准样本点乘以协方差矩阵的平方根来完成的。在下面，我使用scipy.linalg.sqrtm来计算矩阵平方根，使用numpy.einsum来执行矩阵乘法。

import numpy as np
from scipy.linalg import sqrtm
import matplotlib.pyplot as plt

# Sequence of means
means = np.array([[-15., 0.], [15., 0.], [0., 0.]])
# Sequence of covariance matrices.  Must be the same length as means.
covs = np.array([[[ 3, -1],
                  [-1,  2]],
                 [[ 1,  2],
                  [ 2,  5]],
                 [[ 1,  0],
                  [ 0,  1]]])
# Number of samples to generate for each (mean, cov) pair.
nsamples = 4000

# Compute the matrix square root of each covariance matrix.
sqrtcovs = np.array([sqrtm(c) for c in covs])

# Generate samples from the standard multivariate normal distribution.
dim = len(means[0])
u = np.random.multivariate_normal(np.zeros(dim), np.eye(dim),
                                  size=(len(means), nsamples,))
# u has shape (len(means), nsamples, dim)

# Transform u.
v = np.einsum('ijk,ikl->ijl', u, sqrtcovs)
m = np.expand_dims(means, 1)
t = v + m

# t also has shape (len(means), nsamples, dim).
# t[i] holds the nsamples sampled from the distribution with mean means[i]
# and covariance cov[i].

plt.subplot(2, 1, 1)
plt.plot(t[...,0].ravel(), t[...,1].ravel(), '.', alpha=0.02)
plt.axis('equal')
plt.xlim(-25, 25)
plt.ylim(-8, 8)
plt.grid()

# Make another plot, where we generate the samples by passing the given
# means and covs to np.random.multivariate_normal.  This plot should look
# the same as the first plot.
plt.subplot(2, 1, 2)
p0 = np.random.multivariate_normal(means[0], covs[0], size=nsamples)
p1 = np.random.multivariate_normal(means[1], covs[1], size=nsamples)
p2 = np.random.multivariate_normal(means[2], covs[2], size=nsamples)

plt.plot(p0[:,0], p0[:,1], 'b.', alpha=0.02)
plt.plot(p1[:,0], p1[:,1], 'g.', alpha=0.02)
plt.plot(p2[:,0], p2[:,1], 'r.', alpha=0.02)
plt.axis('equal')
plt.xlim(-25, 25)
plt.ylim(-8, 8)
plt.grid()

与循环遍历means和covs数组并为每对数组调用一次multivariate_normal相比，此方法可能不会更快（mean，cov）。这种方法最给予好处的情况是当你有 * 许多 * 不同的均值和协方差，并且每对生成少量样本时。即使这样，它也可能不会更快，因为脚本使用Python循环遍历covs数组，为每个协方差矩阵调用sqrtm。如果性能很关键，请使用实际数据进行测试。

由于我在任何地方都没有找到答案，我只需要为每个mean, std对计算一次pdf(X)。
我直接从公式中矢量化（所以它只适用于pdf（但其他函数也可以类似地编写）：

lpi = (2*np.pi)**3
def vectorized_normal_pdf(X, means, stds):
    ndev = (X - means)/stds
    exp = (ndev[:,None,:] @ (X - means)[:,:,None]).squeeze()
    return np.exp(-0.5*exp)/np.sqrt(lpi*stds.prod(axis=1))

其中，
所有X、means、stds的形状均为[N, d]（N个多变量数据点，每个数据点具有d个值）
并且输出是[N]
我已经验证了它给出了正确的答案（在1e-14的小误差范围内，我不知道为什么它不相等，也许他们在除法中添加了一个小ε，从而使用了一些数值稳定性问题）
并且快得多（在大小仅为10^4的情况下，我获得了**~4300倍**的加速比）：

X = np.random.rand(10000, 3)
means = np.random.rand(10000, 3)
stds = np.random.rand(10000, 3)

>>> %timeit norm_pdf(X, means, stds)
684 µs ± 12.7 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

>>> %timeit [multivariate_normal(means[i], stds[i]).pdf(X[i]) for i in range(10000)]
2.94 s ± 207 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

>>> ( (res1 - res2) < 1e-14 ).all()
True

这对于像高斯混合模型这样应用于图像的应用来说是至关重要的，因为我们需要为每个像素处理多个高斯，所以对于非常小/低分辨率的240 * 320图像，它是76800个高斯。

注意，这不处理协方差矩阵，但通常您可以只使用标准差而不是整个矩阵