5

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

使用 FFT 三点调试 KDE。 使用 SciPy 三点调试 KDE

4

1 回答 1

4

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

默认情况下，使用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 吗？我不确定为什么需要它们，但是没有它们，高斯人会出现在错误的地方。

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

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

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

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

我把那个留给你去弄清楚。:) 不过，这并不难。基本上，scotts_factor您将更改公式并使用其他一些标量因子，而不是 。其他一切都是一样的。

于 2013-09-20T16:52:49.393 回答