在 python 中有效生成 "shifted" 高斯核

Question

我有一个（非常多的）数据点，每个数据点都包含一个 x 和 y 坐标以及一个 sigma 不确定性（sigma 在 x 和 y 方向上相同；所有三个变量都是浮点数）。对于每个数据点，我想在标准网格上生成一个二维数组，实际值在该位置的概率。

例如，如果 x=5.0、y=5.0、sigma=1.0，在 (0,0)->(9,9) 网格上，我希望生成：

   [ 0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ],
   [ 0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ],
   [ 0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ],
   [ 0.  ,  0.  ,  0.  ,  0.  ,  0.01,  0.02,  0.01,  0.  ,  0.  ,  0.  ],
   [ 0.  ,  0.  ,  0.  ,  0.01,  0.06,  0.1 ,  0.06,  0.01,  0.  ,  0.  ],
   [ 0.  ,  0.  ,  0.  ,  0.02,  0.1 ,  0.16,  0.1 ,  0.02,  0.  ,  0.  ],
   [ 0.  ,  0.  ,  0.  ,  0.01,  0.06,  0.1 ,  0.06,  0.01,  0.  ,  0.  ],
   [ 0.  ,  0.  ,  0.  ,  0.  ,  0.01,  0.02,  0.01,  0.  ,  0.  ,  0.  ],
   [ 0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ],
   [ 0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ]]

上面是通过创建一个带有零的 numpy 数组生成的，[5,5] = 1，然后应用 ndimage.filters.gaussian_filter 和 1 的西格玛。我感觉我可以通过分布在附近的整数值上来处理非整数 x 和 y 并得到一个很好的近似值。

然而，以这种方式获得我的结果数组感觉有点矫枉过正，因为 scipy 必须考虑所有值，而不仅仅是位置 [5, 5] 中的 1，即使它们都是 0。对于 64x64 网格只需要 300us，但是，我仍然想知道是否没有更有效的方法来获得具有任意 x、y 和 sigma 的高斯内核的 X*Y numpy 数组。

Answer 1

一种相当快速的方法是注意高斯是可分离的，因此您可以计算 x 和 y 的一维高斯，然后取外积：

import numpy as np
import matplotlib.pyplot as plt

x0, y0, sigma = 5.5, 4.2, 1.4

x, y = np.arange(9), np.arange(9)

gx = np.exp(-(x-x0)**2/(2*sigma**2))
gy = np.exp(-(y-y0)**2/(2*sigma**2))
g = np.outer(gx, gy)
g /= np.sum(g)  # normalize, if you want that

plt.imshow(g, interpolation="nearest", origin="lower")
plt.show()

Answer 2

@tom10 的外积答案可能是这个特定案例的最佳答案。如果您想从二维（或更多）维度的任意函数中创建内核，您可能需要查看 np.indices 或 np.meshgrid.

例如：

def gaussian(x, mu=0, sigma=1):
    n = np.prod(sigma)*np.sqrt(2*np.pi)**len(x)
    return np.exp(-0.5*(((x-mu)/sigma)**2).sum(0))/n

gaussian(np.indices((10,10)), mu=5, sigma=1)

产生：

array([[ 0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   ],
       [ 0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   ],
       [ 0.   , 0.   , 0.   , 0.   , 0.001, 0.002, 0.001, 0.   , 0.   , 0.   ],
       [ 0.   , 0.   , 0.   , 0.003, 0.013, 0.022, 0.013, 0.003, 0.   , 0.   ],
       [ 0.   , 0.   , 0.001, 0.013, 0.059, 0.097, 0.059, 0.013, 0.001, 0.   ],
       [ 0.   , 0.   , 0.002, 0.022, 0.097, 0.159, 0.097, 0.022, 0.002, 0.   ],
       [ 0.   , 0.   , 0.001, 0.013, 0.059, 0.097, 0.059, 0.013, 0.001, 0.   ],
       [ 0.   , 0.   , 0.   , 0.003, 0.013, 0.022, 0.013, 0.003, 0.   , 0.   ],
       [ 0.   , 0.   , 0.   , 0.   , 0.001, 0.002, 0.001, 0.   , 0.   , 0.   ],
       [ 0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   ]])

为了更加灵活，您可以使用 np.meshgrid 来控制域的规模和范围：

kern = gaussian(np.meshgrid(np.linspace(-10, 5), np.linspace(-2, 2)))

为此，kern.shape 将是 (50, 50)，因为 50 是 np.linspace 的默认长度，而 meshgrid 定义 x 和 y 轴由传递给它的数组。这样做的等效方法是 np.mgrid[-10:5:50j, -2:2:50j]

在 python 中有效生成 "shifted" 高斯核

efficiently generate "shifted" gaussian kernel in python

python

numpy

scipy