如何绘制具有不同 sigma 的二维高斯分布？答案

【问题标题】：How to plot a 2d gaussian with different sigma?如何绘制具有不同 sigma 的二维高斯分布？
【发布时间】：2021-02-03 20:10:30
【问题描述】：

我正在尝试制作和绘制具有两个不同标准偏差的二维高斯。他们给出了数学世界的方程式：http://mathworld.wolfram.com/GaussianFunction.html，但我似乎无法得到一个以零为中心的正确二维数组。

我得到了这个，但它不太好用。

x = np.array([np.arange(size)])
y = np.transpose(np.array([np.arange(size)]))

psf  = 1/(2*np.pi*sigma_x*sigma_y) * np.exp(-(x**2/(2*sigma_x**2) + y**2/(2*sigma_y**2)))

【问题讨论】：

您的 x 和 y 轴从 0 开始，因此您只能得到 2D 高斯的四分之一（落在正 xy 象限中。尝试从每个数组中减去 size/2。

标签： python plot gaussian

【解决方案1】：

对于@Coolcrab 来说，这个答案可能为时已晚，但我想把它留在这里以供将来参考。

您可以使用如下的多元高斯公式

更改均值元素会更改原点，而更改协方差元素会更改形状（从圆形变为椭圆）。

代码如下：

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D

# Our 2-dimensional distribution will be over variables X and Y
N = 40
X = np.linspace(-2, 2, N)
Y = np.linspace(-2, 2, N)
X, Y = np.meshgrid(X, Y)

# Mean vector and covariance matrix
mu = np.array([0., 0.])
Sigma = np.array([[ 1. , -0.5], [-0.5,  1.]])

# Pack X and Y into a single 3-dimensional array
pos = np.empty(X.shape + (2,))
pos[:, :, 0] = X
pos[:, :, 1] = Y

def multivariate_gaussian(pos, mu, Sigma):
    """Return the multivariate Gaussian distribution on array pos."""

    n = mu.shape[0]
    Sigma_det = np.linalg.det(Sigma)
    Sigma_inv = np.linalg.inv(Sigma)
    N = np.sqrt((2*np.pi)**n * Sigma_det)
    # This einsum call calculates (x-mu)T.Sigma-1.(x-mu) in a vectorized
    # way across all the input variables.
    fac = np.einsum('...k,kl,...l->...', pos-mu, Sigma_inv, pos-mu)

    return np.exp(-fac / 2) / N

# The distribution on the variables X, Y packed into pos.
Z = multivariate_gaussian(pos, mu, Sigma)

# plot using subplots
fig = plt.figure()
ax1 = fig.add_subplot(2,1,1,projection='3d')

ax1.plot_surface(X, Y, Z, rstride=3, cstride=3, linewidth=1, antialiased=True,
                cmap=cm.viridis)
ax1.view_init(55,-70)
ax1.set_xticks([])
ax1.set_yticks([])
ax1.set_zticks([])
ax1.set_xlabel(r'$x_1$')
ax1.set_ylabel(r'$x_2$')

ax2 = fig.add_subplot(2,1,2,projection='3d')
ax2.contourf(X, Y, Z, zdir='z', offset=0, cmap=cm.viridis)
ax2.view_init(90, 270)

ax2.grid(False)
ax2.set_xticks([])
ax2.set_yticks([])
ax2.set_zticks([])
ax2.set_xlabel(r'$x_1$')
ax2.set_ylabel(r'$x_2$')

plt.show()

【讨论】：

这是对该问题的最完整和最一般的答案。
有趣的代码与 scipy 博客上的代码相同 scipython.com/blog/…

【解决方案2】：

您的函数以零为中心，但您的坐标向量不是。试试：

size = 100
sigma_x = 6.
sigma_y = 2.

x = np.linspace(-10, 10, size)
y = np.linspace(-10, 10, size)

x, y = np.meshgrid(x, y)
z = (1/(2*np.pi*sigma_x*sigma_y) * np.exp(-(x**2/(2*sigma_x**2)
     + y**2/(2*sigma_y**2))))

plt.contourf(x, y, z, cmap='Blues')
plt.colorbar()
plt.show()

【讨论】：

【解决方案3】：

为此，您可以使用 scipy 包中的 multivariate_normal()，如下所示：

# Imports
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal
from mpl_toolkits.mplot3d import Axes3D

# Data
x = np.linspace(-10, 10, 500)
y = np.linspace(-10, 10, 500)
X, Y = np.meshgrid(x,y)

# Multivariate Normal
mu_x = np.mean(x)
sigma_x = np.std(x)
mu_y = np.mean(y)
sigma_y = np.std(y)
rv = multivariate_normal([mu_x, mu_y], [[sigma_x, 0], [0, sigma_y]])

# Probability Density
pos = np.empty(X.shape + (2,))
pos[:, :, 0] = X
pos[:, :, 1] = Y
pd = rv.pdf(pos)

# Plot
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot_surface(X, Y, pd, cmap='viridis', linewidth=0)
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Probability Density')
plt.title("Multivariate Normal Distribution")
plt.show()

【讨论】：

【解决方案4】：

您的高斯以 (0,0) 为中心，因此围绕该原点设置轴。例如，

In [40]: size = 200

In [41]: sigma_x,sigma_y = 50, 20

In [42]: x = np.array([np.arange(size)]) - size/2
In [43]: y = np.transpose(np.array([np.arange(size)])) - size /2
In [44]: psf  = 1/(2*np.pi*sigma_x*sigma_y) *
                np.exp(-(x**2/(2*sigma_x**2) + y**2/(2*sigma_y**2)))

In [45]: pylab.imshow(psf)
Out[45]: <matplotlib.image.AxesImage at 0x10bc07f10>

In [46]: pylab.show()

【讨论】：

【解决方案5】：

我觉得这确实不太友好。我将在这里编写代码并解释它的工作原理。

多元高斯方程如下：

在 2D 情况下，和是 2D 列向量，是 2x2 协方差矩阵并且 n=2。

所以在二维情况下，向量实际上是一个点(x,y)，我们要为其计算函数值，给定二维均值向量，即我们也可以写成(mX, mY)，协方差矩阵。

为了让实现更友好，我们计算的结果：

所以是column 向量(x - mX, y - mY)。所以的计算结果就是二维行向量：

(CI[0,0] * (x - mX) + CI[1,0] * (y - mY) , CI[0,1] * (x - mX) + CI[1, 1] * (y - mY))，其中CI是协方差矩阵的逆，在方程中显示为，它是一个2x2矩阵，如是。

然后，将当前结果（二维行向量）乘以列向量（内积），最终得到标量：

CI[0,0](x - mX)^2 + (CI[1,0] + CI[0,1])(y - mY)(x - mX) + CI[1, 1](y - mY)^2

与相比，使用 NumPy 实现此表达式将更容易，即使它们具有相同的值。

>>> m = np.array([[0.2],[0.6]])  # defining the mean of the Gaussian (mX = 0.2, mY=0.6)
>>> cov = np.array([[0.7, 0.4], [0.4, 0.25]])   # defining the covariance matrix
>>> cov_inv = np.linalg.inv(cov)  # inverse of covariance matrix
>>> cov_det = np.linalg.det(cov)  # determinant of covariance matrix
# Plotting
>>> x = np.linspace(-2, 2)
>>> y = np.linspace(-2, 2)
>>> X,Y = np.meshgrid(x,y)
>>> coe = 1.0 / ((2 * np.pi)**2 * cov_det)**0.5
>>> Z = coe * np.e ** (-0.5 * (cov_inv[0,0]*(X-m[0])**2 + (cov_inv[0,1] + cov_inv[1,0])*(X-m[0])*(Y-m[1]) + cov_inv[1,1]*(Y-m[1])**2))
>>> plt.contour(X,Y,Z)
<matplotlib.contour.QuadContourSet object at 0x00000252C55773C8>
>>> plt.show()

结果：

【讨论】：