在Python中生成3D高斯分布

Question

我想在Python中生成一个高斯分布，其中x和y代表位置，z代表某个量的大小。

这个分布的最大值是2e6，标准差sigma为0.025。

在MATLAB中，我可以这样做：

x1 = linspace(-1,1,30);
x2 = linspace(-1,1,30);

mu = [0,0];
Sigma = [.025,.025];

[X1,X2] = meshgrid(x1,x2);
F = mvnpdf([X1(:) X2(:)],mu,Sigma);
F = 314159.153*reshape(F,length(x2),length(x1));
surf(x1,x2,F);

在Python中，我目前有的代码是：

x = np.linspace(-1,1,30)
y = np.linspace(-1,1,30)

mu = (np.median(x),np.median(y))

sigma = (.025,.025)

有一个Numpy的函数叫做numpy.random.multivariate_normal，听说它可以做和MATLAB的mvnpdf一样的事情，但我对文档有点看不懂。特别是关于如何获取numpy.random.multivariate_normal所需的协方差矩阵。

Answer 1

我正在使用一个叫做 scikit-guess 的工具包，它里面有一些快速估算非线性拟合的功能。这个工具包里有一个函数 skg.ngauss.model（也可以用 skg.ngauss_fit.model 或 skg.ngauss.ngauss_fit.model 来访问），它正好能满足你的需求。最棒的是，它不是一个概率密度函数（PDF），所以你可以直接设置振幅：

import numpy as np
import skg.ngauss

a = 2e6
mu = 0, 0
sigma = 0.025, 0.025

x = y = np.linspace(-1, 1, 31)

cov = np.diag(sigma)**2
X = np.meshgrid(x, y)

data = skg.ngauss.model(X, a, mu, cov, axis=0)

你需要告诉它 axis=0，因为它会自动帮你把数组堆叠起来。如果你不想传这个参数，可以这样写：

X = np.stack(np.meshgrid(x, y), axis=-1)

你可以把结果画出来：

from matplotlib import pyplot as plt
plt.imshow(data)
plt.show()

这个分布看起来不太有趣，因为它的分散度太小，导致你在一个像素的距离内就得到了大约 2e-5 的值。你可能需要扩大采样范围，才能得到更有意义的分辨率。

注意：在写这段话的时候，拟合函数（ngauss_fit）还有一些bug，但模型已经成功测试过，只是还没有在这个工具包里使用。

免责声明：如果上面没有说明，我是 scikit-guess 的作者。

Answer 2

从scipy 0.14版本开始，你可以使用 scipy.stats.multivariate_normal.pdf() 这个功能。

import numpy as np
from scipy.stats import multivariate_normal

x, y = np.mgrid[-1.0:1.0:30j, -1.0:1.0:30j]
# Need an (N, 2) array of (x, y) pairs.
xy = np.column_stack([x.flat, y.flat])

mu = np.array([0.0, 0.0])

sigma = np.array([.025, .025])
covariance = np.diag(sigma**2)

z = multivariate_normal.pdf(xy, mean=mu, cov=covariance)

# Reshape back to a (30, 30) grid.
z = z.reshape(x.shape)