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

Question

我需要代码来做二维核密度估计(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定义标量带宽？高斯在两种实现中具有不同的宽度.
• 同样,为什么SciPy实现中的高斯不是径向对称的,即使我给了gaussian_kde一个标量带宽？
• 如何实现SciPy中可用于FFT代码的其他带宽方法？

(请注意,在1000个随机点的情况下,FFT代码比SciPy代码快〜390倍.)

Answer 1

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

默认情况下，使用Scott 规则gaussian_kde选择带宽。如果您对细节感到好奇， 请深入研究代码。下面的代码片段来自我不久前编写的代码，用于执行与您正在执行的操作类似的操作。（如果我没记错的话，该特定版本中有一个明显的错误，它确实不应该用于卷积，但带宽估计​​和归一化是正确的。）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

希望这有助于澄清事情。

至于你的其他问题：

我可以避免直方图的 .T 和 KDE1 的 fftshift 吗？我不确定为什么需要它们，但是如果没有它们，高斯函数就会出现在错误的地方。

我可能会误读您的代码，但我认为您只是进行了转置，因为您要从点坐标到索引坐标（即从<x, y>到<y, x>）。

同样，为什么 SciPy 实现中的高斯不是径向对称的，即使我给了 gaussian_kde 一个标量带宽？

这是因为 scipy 使用输入 x, y 点的完整协方差矩阵来确定高斯核。您的公式假设 x 和 y 不相关。gaussian_kde测试并使用结果中 x 和 y 之间的相关性。

如何为 FFT 代码实现 SciPy 中可用的其他带宽方法？

我会把这个问题留给你去解决。:) 不过，这并不太难。基本上，scotts_factor您可以更改公式并使用其他一些标量因子，而不是 。其他一切都一样。