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

我需要执行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倍。）