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

Dave 的回答不正确，因为scipy'svonmises没有环绕[-pi, pi]。

相反，您可以使用以下基于相同原理的代码。它基于numpy 中描述的方程。

def vonmises_kde(data, kappa, n_bins=100):
    from scipy.special import i0
    bins = np.linspace(-np.pi, np.pi, n_bins)
    x = np.linspace(-np.pi, np.pi, n_bins)
    # integrate vonmises kernels
    kde = np.exp(kappa*np.cos(x[:, None]-data[None, :])).sum(1)/(2*np.pi*i0(kappa))
    kde /= np.trapz(kde, x=bins)
    return bins, kde

这是一个例子

import matplotlib.pyplot as plt
import numpy as np
from numpy.random import vonmises

# generate complex circular distribution
data = np.r_[vonmises(-1, 5, 1000), vonmises(2, 10, 500), vonmises(3, 20, 100)]

# plot data histogram
fig, axes = plt.subplots(2, 1)
axes[0].hist(data, 100)

# plot kernel density estimates
x, kde = vonmises_kde(data, 20)
axes[1].plot(x, kde)

直方图和核密度图

Answer 2

这是一个快速近似这是@kingjr 更准确答案的

\n\n

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

\n\n

测试（使用 tqdm 进行进度条和计时，使用 matplotlib 验证结果）：

\n\n

import numpy as np\nfrom tqdm import tqdm\nimport scipy.stats\nimport matplotlib.pyplot as plt\n\nn_runs = 1000\nn_bins = 100\nkappa = 10\n\nfor _ in tqdm(xrange(n_runs)):\n    bins1, kde1 = vonmises_kde(\n        data=np.r_[\n            np.random.vonmises(-1, 5, 1000),\n            np.random.vonmises(2, 10, 500),\n            np.random.vonmises(3, 20, 100)\n        ],\n        kappa=kappa,\n        n_bins=n_bins\n    )\n\n\nfor _ in tqdm(xrange(n_runs)):\n    bins2, kde2 = vonmises_fft_kde(\n        data=np.r_[\n            np.random.vonmises(-1, 5, 1000),\n            np.random.vonmises(2, 10, 500),\n            np.random.vonmises(3, 20, 100)\n        ],\n        kappa=kappa,\n        n_bins=n_bins\n    )\n\nplt.figure()\nplt.plot(bins1, kde1, label="kingjr\'s solution")\nplt.plot(bins2, kde2, label="dolf\'s FFT solution")\nplt.legend()\nplt.show()\n

\n\n

结果：

\n\n

100%|\xe2\x96\x88\xe2\x96\x88\xe2\x96\x88\xe2\x96\x88\xe2\x96\x88\xe2\x96\x88\xe2\x96\x88\xe2\x96\x88\xe2\x96\x88\xe2\x96\x88| 1000/1000 [00:07<00:00, 135.29it/s]\n100%|\xe2\x96\x88\xe2\x96\x88\xe2\x96\x88\xe2\x96\x88\xe2\x96\x88\xe2\x96\x88\xe2\x96\x88\xe2\x96\x88\xe2\x96\x88\xe2\x96\x88| 1000/1000 [00:00<00:00, 1945.14it/s]\n

\n\n

（1945 / 135 = 快 14 倍）

\n\n

\n\n

为了获得更快的速度，请使用 2 的整数次方作为 bin 的数量。它还具有更好的扩展性（即，它可以在处理许多数据箱和大量数据时保持快速）。在我的 PC 上，n_bins=1024 时，它比原始答案快 118 倍。

\n\n

为什么它有效？

\n\n

两个信号（无补零）的 FFT 乘积等于两个信号的循环（或循环）卷积。核密度估计基本上是与在每个数据点的位置处具有脉冲的信号进行卷积的核。

\n\n

为什么不准确？

\n\n

由于我使用直方图均匀地间隔数据，因此我丢失了每个样本的确切位置，并且仅使用它所属的箱的中心。每个箱中的样本数量用作该点的脉冲幅度。例如：暂时忽略归一化，如果您有一个从 0 到 1 的 bin，并且该 bin 中有两个样本，分别为 0.1 和 0.2，则 KDEexact将为the kernel centred around 0.1+ the kernel centred around 0.2。近似值为 2x“内核以 0.5 为中心，即 bin 的中心。

\n