二维高斯分布总和不等于一？

Question

我使用此处给出的等式在 Python 中构建了一个包裹的二元高斯分布： http://www.aos.wisc.edu/~dvimont/aos575/Handouts/bivariate_notes.pdf 但是，我不明白为什么我的分布尽管包含了一个归一化常数，但总和却不能为 1。

对于 U x U 晶格，

import numpy as np
from math import *

U = 60
m = np.arange(U)
i = m.reshape(U,1)
j = m.reshape(1,U)

sigma = 0.1
ii = np.minimum(i, U-i)
jj = np.minimum(j, U-j)
norm_constant = 1/(2*pi*sigma**2)
xmu = (ii-0)/sigma; ymu = (jj-0)/sigma
rhs = np.exp(-.5 * (xmu**2 + ymu**2))
ker = norm_constant * rhs

>> ker.sum() # area of each grid is 1 
15.915494309189533

我确信我的思考方式从根本上缺失，并怀疑需要进行某种额外的规范化，尽管我无法绕过它。

更新：

感谢其他人有见地的建议，我重写了我的代码以将 L1 规范化应用于内核。但是，似乎在通过 FFt 进行 2D 卷积的情况下，将范围保持为 [0, U] 仍然能够返回令人信服的结果：

U = 100
Ukern = np.copy(U)
#Ukern = 15

m = np.arange(U)
i = m.reshape(U,1)
j = m.reshape(1,U)

sigma = 2.
ii = np.minimum(i, Ukern-i)
jj = np.minimum(j, Ukern-j)
xmu = (ii-0)/sigma; ymu = (jj-0)/sigma
ker = np.exp(-.5 * (xmu**2 + ymu**2))
ker /= np.abs(ker).sum()

''' Point Density '''
ido = np.random.randint(U, size=(10,2)).astype(np.int)
og = np.zeros((U,U))
np.add.at(og, (ido[:,0], ido[:,1]), 1)

''' Convolution via FFT and inverse-FFT '''
v1 = np.fft.fft2(ker)
v2 = np.fft.fft2(og)
v0 = np.fft.ifft2(v2*v1)
dd = np.abs(v0)

plt.plot(ido[:,1], ido[:,0], 'ko', alpha=.3)
plt.imshow(dd, origin='origin')
plt.show()

另一方面，使用注释掉的行调整内核大小会给出这个错误的图：

Answer 1

注意：如下面的 cmets 所述，此解决方案仅在您尝试构建高斯卷积核（或高斯滤波器）用于图像处理时才有效。它不是一个适当归一化的高斯密度函数，但它是用于从图像中去除高斯噪声的形式。

---

您缺少 L1 标准化：

ker /= np.abs(ker).sum()

这将使您的内核表现得像一个实际的密度函数。由于您拥有的网格的值大小可能会有很大差异，因此需要上述标准化步骤。

事实上，你可以省略高斯归一化常数，只使用上面的 L1 范数。如果我没穿，你正在尝试创建一个高斯卷积，上面是应用于它的常用归一化技术。

正如@Praveen 所说，您的第二个错误是您需要从[-U//2, U//2] 中采样网格。你可以这样做：

i, j = np.mgrid[-U//2:U//2+1, -U//2:U//2+1]

最后，如果您要做的是构建一个高斯滤波器，内核的大小通常从 sigma 估计（以避免远离中心的零）为U//2 <= t * sigma，其中t 是截断参数通常设置t=3或t=4。

Answer 2

目前，ker 的（放大很多）等值线图如下所示：

如您所见，这看起来一点也不像高斯核。您的大多数函数从 0 到 1 消失。查看内核本身会发现所有值确实真的很快消失：

>>> ker[0:5, 0:5]
array([[  1.592e+001,   3.070e-021,   2.203e-086,   5.879e-195,   0.000e+000],
       [  3.070e-021,   5.921e-043,   4.248e-108,   1.134e-216,   0.000e+000],
       [  2.203e-086,   4.248e-108,   3.048e-173,   8.136e-282,   0.000e+000],
       [  5.879e-195,   1.134e-216,   8.136e-282,   0.000e+000,   0.000e+000],
       [  0.000e+000,   0.000e+000,   0.000e+000,   0.000e+000,   0.000e+000]])

你得到的 15.915 的总和值基本上就是 ker[0, 0]。这一切都告诉你，你没有正确地构建你的网格。

请记住，在计算机上创建内核时，您必须在适当的位置对其进行采样。采样太粗会导致你的总和不正确。

首先，如果您想要以mu=0 为中心的全密度，您必须将i 和j 从-U // 2 带到U // 2。但要解决您的分辨率问题，我建议将U 的点数取在-0.5 到0.5 之间。

import numpy as np
import matplotlib.pyplot as plt

U = 60
m = np.linspace(-0.5, 0.5, U)    # 60 points between -1 and 1
delta = m[1] - m[0]              # delta^2 is the area of each grid cell
(x, y) = np.meshgrid(m, m)       # Create the mesh

sigma = 0.1
norm_constant = 1 / (2 * np.pi * sigma**2)

rhs = np.exp(-.5 * (x**2 + y**2) / sigma**2)
ker = norm_constant * rhs
print(ker.sum() * delta**2)

plt.contour(x, y, ker)
plt.axis('equal')
plt.show()

这会产生接近 1.0 的总和，并且内核以mu=0 为中心，正如预期的那样。

在这种情况下，了解选择的范围（-0.5 到 0.5）取决于您的功能。例如，如果您现在采用sigma = 2，您会发现您的总和不会再次计算出来，因为现在您的采样太精细了。将范围设置为参数的函数（例如 -5 * sigma 到 5 * sigma）可能是最好的选择。