在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所需的协方差矩阵方面。

- Jonny

---

我认为你错了，认为numpy.random.multivariate_normal()做的是同样的事情，因为它并没有给出分布的概率密度函数，它只是从协方差矩阵和期望值mu定义的分布中绘制随机数。 - Nras

我明白你的意思。那么，你有什么建议来实现它吗？ - Jonny

我看到你的_xy_分布是可分离的，也就是说，它是一个_x_高斯分布和一个_y_高斯分布的乘积。也许这可以帮助你使用Python。 - Luis Mendo

Answer 1

从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)

- Robert Kern

Answer 2

我正在开发一个名为scikit-guess的scikit，其中包含一些用于非线性拟合的快速估计例程。它有一个名为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）仍存在错误，但该模型已经成功测试，只是没有在scikit中进行。

免责声明：如果以上内容不明显，我是scikit-guess的作者。

- Mad Physicist