scipy gaussian_kde和循环数据

Question

我正在使用scipys gaussian_kde获取某些双峰数据的概率密度。但是，由于我的数据是有角度的(以度为单位的方向)，所以当值接近极限时会出现问题。下面的代码提供了两个示例kde，当域为0-360时，由于无法处理数据的循环性质，因此处于估计状态。 pdf需要在单位圆上定义，但我在scipy.stats中找不到适合此类数据的任何内容(存在冯·米斯分布，但仅适用于单峰数据)。外面有没有人遇到过这个？是否有任何可用于估算单位圆上的双峰pdf的(基于python的优先选择)？

import numpy as np
import scipy as sp
from pylab import plot,figure,subplot,show,hist
from scipy import stats



baz = np.array([-92.29061004, -85.42607874, -85.42607874, -70.01689348,
               -63.43494882, -63.43494882, -70.01689348, -70.01689348,
               -59.93141718, -63.43494882, -59.93141718, -63.43494882,
               -63.43494882, -63.43494882, -57.52880771, -53.61564818,
               -57.52880771, -63.43494882, -63.43494882, -92.29061004,
               -16.92751306, -99.09027692, -99.09027692, -16.92751306,
               -99.09027692, -16.92751306,  -9.86580694,  -8.74616226,
                -9.86580694,  -8.74616226,  -8.74616226,  -2.20259816,
                -2.20259816,  -2.20259816,  -9.86580694,  -2.20259816,
                -2.48955292,  -2.48955292,  -2.48955292,  -2.48955292,
                 4.96974073,   4.96974073,   4.96974073,   4.96974073,
                -2.48955292,  -2.48955292,  -2.48955292,  -2.48955292,
                -2.48955292,  -9.86580694,  -9.86580694,  -9.86580694,
               -16.92751306, -19.29004622, -19.29004622, -26.56505118,
               -19.29004622, -19.29004622, -19.29004622, -19.29004622])


xx = np.linspace(-180, 180, 181)
scipy_kde = stats.gaussian_kde(baz)              
print scipy_kde.integrate_box_1d(-180,180)

figure()
plot(xx, scipy_kde(xx), c='green')             

baz[baz<0] += 360             
xx = np.linspace(0, 360, 181)
scipy_kde = stats.gaussian_kde(baz)              
print scipy_kde.integrate_box_1d(-180,180)
plot(xx, scipy_kde(xx), c='red')             

Answer 1

这是对@kingjr的更准确答案的快速近似:

def vonmises_pdf(x, mu, kappa):
    return np.exp(kappa * np.cos(x - mu)) / (2. * np.pi * scipy.special.i0(kappa))


def vonmises_fft_kde(data, kappa, n_bins):
    bins = np.linspace(-np.pi, np.pi, n_bins + 1, endpoint=True)
    hist_n, bin_edges = np.histogram(data, bins=bins)
    bin_centers = np.mean([bin_edges[1:], bin_edges[:-1]], axis=0)
    kernel = vonmises_pdf(
        x=bin_centers,
        mu=0,
        kappa=kappa
    )
    kde = np.fft.fftshift(np.fft.irfft(np.fft.rfft(kernel) * np.fft.rfft(hist_n)))
    kde /= np.trapz(kde, x=bin_centers)
    return bin_centers, kde

测试(将tqdm用于进度条和计时，并使用matplotlib验证结果):

import numpy as np
from tqdm import tqdm
import scipy.stats
import matplotlib.pyplot as plt

n_runs = 1000
n_bins = 100
kappa = 10

for _ in tqdm(xrange(n_runs)):
    bins1, kde1 = vonmises_kde(
        data=np.r_[
            np.random.vonmises(-1, 5, 1000),
            np.random.vonmises(2, 10, 500),
            np.random.vonmises(3, 20, 100)
        ],
        kappa=kappa,
        n_bins=n_bins
    )


for _ in tqdm(xrange(n_runs)):
    bins2, kde2 = vonmises_fft_kde(
        data=np.r_[
            np.random.vonmises(-1, 5, 1000),
            np.random.vonmises(2, 10, 500),
            np.random.vonmises(3, 20, 100)
        ],
        kappa=kappa,
        n_bins=n_bins
    )

plt.figure()
plt.plot(bins1, kde1, label="kingjr's solution")
plt.plot(bins2, kde2, label="dolf's FFT solution")
plt.legend()
plt.show()

结果:

100%|██████████| 1000/1000 [00:07<00:00, 135.29it/s]
100%|██████████| 1000/1000 [00:00<00:00, 1945.14it/s]

(1945/135 =快14倍)



为了获得更高的速度，请使用2的整数次方作为箱数。它的伸缩性也更好(即在许多分档和大量数据的情况下保持快速运行)。在我的电脑上，它的速度是n_bins = 1024时原始答案的118倍。

为什么有效？

两个信号的FFT乘积(无零填充)等于两个信号的 circular (or cyclic) convolution。 kernel density estimation本质上是一个卷积有信号的内核，该信号在每个数据点的位置都有一个脉冲。

为什么不准确？

由于我使用直方图来均匀分布数据，因此我丢失了每个样本的确切位置，而只使用了它所属的bin的中心。每个仓中的样本数量用作该点的脉冲幅度。例如:暂时忽略归一化，如果您有一个0到1的bin，并且该bin中有两个样本(分别为0.1和0.2)，则 exact KDE将为 the kernel centred around 0.1 + the kernel centred around 0.2。近似值是2x内核在0.5的中心，该中心是bin的中心。