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

概述我需要代码来做二维核密度估计(KDE),我发现SciPy实现太慢了.所以,我已经编写了一个基于FFT的实现,但有些事情让我很困惑. (FFT实现还强制执行周期性边界条件,这就是我想要的.)该实现基于从样本创建简单的直方图,然后使用高斯进行卷积.这是执行此 *** 作的代码,并将其与SciPy结果进行比较.from numpy import * from scipy.

我需要代码来做二维核密度估计(KDE),我发现SciPy实现太慢了.所以,我已经编写了一个基于FFT的实现,但有些事情让我很困惑. (FFT实现还强制执行周期性边界条件,这就是我想要的.)

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

from numpy import *from scipy.stats import *from numpy.fft import *from matplotlib.pyplot import *from time import clockion()#ParaMETERSN   = 512   #number of histogram bins; want 2^n for maximum FFT speed?nSamp   = 1000  #number of samples if using the ranom variableh   = 0.1   #wIDth of gaussianwh  = 1.0   #wIDth and height of square domain#VARIABLES FROM ParaMETERSrv  = uniform(loc=-wh,scale=2*wh)   #random variable that can generate samplesxyBnds  = linspace(-1.0,1.0,N+1)  #boundarIEs of histogram binsxy  = (xyBnds[1:] + xyBnds[:-1])/2      #centers of histogram binsxx,yy = meshgrID(xy,xy)#define SAMPLES,TWO OPTIONS#samples = rv.rvs(size=(nSamp,2))samples = array([[0.5,0.5],[0.2,0.2]])#DEFinitioNS FOR FFT IMPLEMENTATIONker = exp(-(xx**2 + yy**2)/2/h**2)/h/sqrt(2*pi) #Gaussian kernelfKer = fft2(ker) #DFT of kernel#FFT IMPLEMENTATIONstime = 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))/Netime = 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,axis=1)#NON-FFT IMPLEMTATION FROM SCIPYstime = 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 RESulTSfig = figure()ax = fig.add_subplot(111,aspect='equal')c = contour(xy,xy,KDE1.real)clabel(c)Title("FFT Implementation Results")#PRINT SCIPY IMPLEMENTATION RESulTSfig = figure()ax = fig.add_subplot(111,KDE2)clabel(c)Title("SciPy Implementation Results")

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

后一种情况的结果图是在这篇文章的最后.

这是我的问题：

>我可以避免直方图的.T和KDE1的fftshift吗？我不确定他们为什么需要它们,但是如果没有它们,高斯人就会出现在错误的地方.
>如何为SciPy定义标量带宽？高斯在两种实现中具有不同的宽度.
>沿着同样的路线,为什么SciPy实现中的高斯不是径向对称的,即使我给了gaussian_kde一个标量带宽？
>我如何实现SciPy中可用于FFT代码的其他带宽方法？

(请注意,在1000个随机点的情况下,FFT代码比SciPy代码快〜390倍.)最佳答案正如您已经注意到的那样,您所看到的差异是由带宽和缩放因素造成的.

默认情况下,如果您对细节感到好奇,gaussian_kde会使用Scott’s rule. Dig into the code选择带宽.下面的代码片段是从我写的quite awhile ago to do something similar到你正在做的事情. (如果我没记错的话,那个特定版本有一个明显的错误,它确实不应该使用scipy.signal进行卷积,但带宽估计和规范化是正确的.)

# Calculate the covariance matrix (in pixel coords)cov = np.cov(xyi)# Scaling factor for banDWIDthscotts_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 kernelinv_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**2norm_factor = np.linalg.det(norm_factor)norm_factor = n * dx * dy * np.sqrt(norm_factor)# normalize the resultgrID /= norm_factor

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

至于你的其他问题：

Can I avoID the .T for the histogram and the fftshift for KDE1? I’m
not sure why they’re needed,but the gaussians show up in the wrong
places without them.

我可能会误读你的代码,但我认为你只是有转置因为你从点坐标到索引坐标(即从< x,y>到< y,x>).