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

前端之家收集整理的这篇文章主要介绍了在python中实现基于FFT的基于FFT的核密度估计器,并将其与SciPy实现进行比较前端之家小编觉得挺不错的,现在分享给大家,也给大家做个参考。

我需要代码来做二维核密度估计(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.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,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,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 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

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

至于你的其他问题:

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>).

Along the same lines,why are the gaussians in the SciPy
implementation not radially symmetric even though I gave gaussian_kde
a scalar bandwidth?

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

How could I implement the other bandwidth methods available in SciPy
for the FFT code?

我会留下那个让你弄明白的. :)但这并不难.基本上,你不是改变scotts_factor,而是改变公式,并有一些其他的标量因子.其他一切都是一样的.

原文链接:https://www.f2er.com/python/439622.html

猜你在找的Python相关文章