我使用scipys gaussian_kde来得到一些双峰数据的概率密度。但是,由于我的数据是有Angular 的(它是以度为单位的方向)当值出现在极限附近时,我遇到了一个问题。下面的代码给出了两个kde的例子,当定义域为0-360时,由于无法处理数据的圆形性质,因此估计值偏低。pdf需要在单位圆上定义,但我不能在scipy.stats中找不到任何适合这种类型数据的东西(冯米塞斯分布在那里,但只适用于单峰数据)。以前有人遇到过这种情况吗?有什么东西(最好是基于python的)可以用来估计单位圆上的双峰pdf吗?
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')
4条答案
按热度按时间vpfxa7rd1#
Dave的答案是不正确的,因为
scipy
的vonmises
没有环绕[-pi, pi]
。你可以使用下面的代码,它基于同样的原理,基于numpy中描述的方程。
下面是一个例子
fzwojiic2#
以下是@kingjr给出的更准确答案的快速近似:
测试(使用tqdm进行进度条和计时,使用matplotlib验证结果):
结果:
(1945/ 135 =快14倍)
为了获得更快的速度,使用2的整数幂作为bin的数量。它的伸缩性也更好(即,它在许多bin和大量数据的情况下保持快速)。在我的PC上,它比原始答案(n_bins=1024)快118倍。
为什么会成功
两个信号(* 无 * 零填充)的FFT的乘积等于两个信号的circular (or cyclic) convolution。kernel density estimation基本上是与在每个数据点的位置处具有脉冲的信号卷积的核。
"为什么不准确"
由于我使用直方图来均匀地分布数据,因此我丢失了每个样本的确切位置,而只使用它所属的bin的中心。每个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的中心。jfgube3f3#
所以我有一个我认为合理的解决方案。基本上我使用冯米塞斯分布作为核密度估计的基函数。代码如下,以防对其他人有用。
6ie5vjzr4#
作为@kingjr回答的一个小而重要的补充:我预计许多开发人员会希望以极坐标投影显示此KDE,如