我需要代码来做二维核密度估计(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代码的其他带宽方法？
默认情况下,如果您对细节感到好奇,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,而是改变公式,并有一些其他的标量因子.其他一切都是一样的.