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")
上記の 2 つのサンプル セットがあります。1000 のランダム ポイントはベンチマーク用であり、コメント アウトされています。3 つのポイントはデバッグ用です。
後者の場合の結果のプロットは、この投稿の最後にあります。
ここに私の質問があります：
- ヒストグラムの .T と KDE1 の fftshift を避けることはできますか? なぜ必要なのかはわかりませんが、ガウス分布がないと間違った場所に表示されます。
- SciPy のスカラー帯域幅はどのように定義されていますか? 2 つの実装では、ガウス分布の幅が大きく異なります。
- 同様に、gaussian_kde にスカラー帯域幅を指定したにもかかわらず、SciPy 実装のガウス分布が放射状に対称でないのはなぜですか?
- FFT コード用に SciPy で利用可能な他の帯域幅メソッドを実装するにはどうすればよいですか?
(1000 ランダム ポイントの場合、FFT コードは SciPy コードよりも ~390 倍高速であることに注意してください。)
