在 python 中实现 2D、基于 FFT 的核密度估计器，并将其与 SciPy 实现进行比较

Posted 2023-03-12

技术标签:

【中文标题】在 python 中实现 2D、基于 FFT 的核密度估计器，并将其与 SciPy 实现进行比较【英文标题】：Implementing a 2D, FFT-based Kernel Density Estimator in python, and comparing it to the SciPy implimentation 【发布时间】：2013-09-26 02:47:25 【问题描述】：

我需要代码来执行 2D 核密度估计 (KDE)，但我发现 SciPy 实现太慢了。所以，我写了一个基于 FFT 的实现，但有几件事让我感到困惑。（FFT 实现还强制执行周期性边界条件，这正是我想要的。）

该实现基于从样本创建一个简单的直方图，然后将其与高斯卷积。这是执行此操作并将其与 SciPy 结果进行比较的代码。

from numpy import *
from scipy.stats import *
from numpy.fft import *
from matplotlib.pyplot import *
from time import clock

ion()

#PARAMETERS
N   = 512   #number of histogram bins; want 2^n for maximum FFT speed?
nSamp   = 1000  #number of samples if using the ranom variable
h   = 0.1   #width of gaussian
wh  = 1.0   #width and height of square domain

#VARIABLES FROM PARAMETERS
rv  = uniform(loc=-wh,scale=2*wh)   #random variable that can generate samples
xyBnds  = linspace(-1.0, 1.0, N+1)  #boundaries of histogram bins
xy  = (xyBnds[1:] + xyBnds[:-1])/2      #centers of histogram bins
xx, yy = meshgrid(xy,xy)

#DEFINE SAMPLES, TWO OPTIONS
#samples = rv.rvs(size=(nSamp,2))
samples = array([[0.5,0.5],[0.2,0.5],[0.2,0.2]])

#DEFINITIONS FOR FFT IMPLEMENTATION
ker = exp(-(xx**2 + yy**2)/2/h**2)/h/sqrt(2*pi) #Gaussian kernel
fKer = fft2(ker) #DFT of kernel

#FFT IMPLEMENTATION
stime = clock()
#generate normalized histogram. Note sure why .T is needed:
hst = histogram2d(samples[:,0], samples[:,1], bins=xyBnds)[0].T / (xy[-1] - xy[0])**2
#convolve histogram with kernel. Not sure why fftshift is neeed:
KDE1 = fftshift(ifft2(fft2(hst)*fKer))/N
etime = clock()
print "FFT method time:", etime - stime

#DEFINITIONS FOR NON-FFT IMPLEMTATION FROM SCIPY
#points to sample the KDE at, in a form gaussian_kde likes:
grid_coords = append(xx.reshape(-1,1),yy.reshape(-1,1),axis=1)

#NON-FFT IMPLEMTATION FROM SCIPY
stime = clock()
KDEfn = gaussian_kde(samples.T, bw_method=h)
KDE2 = KDEfn(grid_coords.T).reshape((N,N))
etime = clock()
print "SciPy time:", etime - stime

#PLOT FFT IMPLEMENTATION RESULTS
fig = figure()
ax = fig.add_subplot(111, aspect='equal')
c = contour(xy, xy, KDE1.real)
clabel(c)
title("FFT Implementation Results")

#PRINT SCIPY IMPLEMENTATION RESULTS
fig = figure()
ax = fig.add_subplot(111, aspect='equal')
c = contour(xy, xy, KDE2)
clabel(c)
title("SciPy Implementation Results")

上面有两组样本。 1000个随机点用于基准测试并被注释掉；三点是调试用的。

后一种情况的结果图在本文末尾。

这是我的问题：

我可以避免直方图的 .T 和 KDE1 的 fftshift 吗？我不确定为什么需要它们，但没有它们，高斯人会出现在错误的地方。如何为 SciPy 定义标量带宽？高斯在两种实现中的宽度有很大不同。同理，为什么即使我给 gaussian_kde 一个标量带宽，SciPy 实现中的高斯函数也不是径向对称的？如何为 FFT 代码实现 SciPy 中可用的其他带宽方法？

（请注意，在 1000 个随机点的情况下，FFT 代码比 SciPy 代码快约 390 倍。）

【问题讨论】：

对于问题 2 和 3，您最好查看 scipy 代码。它使用数据协方差矩阵和 Scott 或 Silverman 带宽规则。我相信上述实现显示了关于核密度估计的非常规选择，例如检查Tarn Duong's explanation；要点是您不应该对直方图进行卷积，而是对原始数据（即位于数据点的增量函数的组合）进行卷积。再一次，一如既往——我可能错了；见jakevdp's FFT-based-computation 【参考方案1】：

正如您已经注意到的那样，您看到的差异是由于带宽和缩放因素造成的。

默认情况下，gaussian_kde 使用Scott's rule. 选择带宽挖掘into the code，如果您对细节感到好奇。下面的代码 sn-ps 来自我写的 quite awhile ago to do something similar 到你正在做的事情。（如果我没记错的话，那个特定版本有一个明显的错误，它确实不应该使用scipy.signal 进行卷积，但带宽估计和归一化是正确的。）

# Calculate the covariance matrix (in pixel coords)
cov = np.cov(xyi)

# Scaling factor for bandwidth
scotts_factor = np.power(n, -1.0 / 6) # For 2D

#---- Make the gaussian kernel -------------------------------------------

# First, determine how big the gridded kernel needs to be (2 stdev radius) 
# (do we need to convolve with a 5x5 array or a 100x100 array?)
std_devs = np.diag(np.sqrt(cov))
kern_nx, kern_ny = np.round(scotts_factor * 2 * np.pi * std_devs)

# Determine the bandwidth to use for the gaussian kernel
inv_cov = np.linalg.inv(cov * scotts_factor**2)

卷积之后，网格被归一化：

# Normalization factor to divide result by so that units are in the same
# units as scipy.stats.kde.gaussian_kde's output.  (Sums to 1 over infinity)
norm_factor = 2 * np.pi * cov * scotts_factor**2
norm_factor = np.linalg.det(norm_factor)
norm_factor = n * dx * dy * np.sqrt(norm_factor)

# Normalize the result
grid /= norm_factor

希望这有助于澄清一些事情。

至于你的其他问题：