Scipy Circular Variaance

Question

根据我的理解,循环方差的范围在0到1之间.这也在维基百科 以及这里得到了证实.但由于某些原因,圆形方差函数scipy.stats给出的值大于1.

import numpy as np
from scipy.stats import circmean, circvar

a = np.random.randint(0, high=360, size=10)

print(a)
print(circmean(a, 0, 360))
print(circvar(np.deg2rad(a)))
[143 116 152 172 349 152 182 306 345  81]
135.34974541954665
2.2576538466653857

有人可以告诉我为什么我从函数中获得高于1的值 circvar

Answer 1

不太有用的答案是,因为scipy定义它的方式,所以你最好让开发人员得到明确的答案.真.来自文档的例子是

from scipy.stats import circvar
circvar([0, 2*np.pi/3, 5*np.pi/3])
2.19722457734

所以你不能说这种行为是未被发现的.但为什么这样做呢？

您的第二个链接定义了一组n个角度a_1,... a_n的圆形方差

V = 1 - \hat {R_1}

哪里

\ hat {R_1} = R_1/n R_1 =\sqrt {C ^ 2 + S ^ 2}

和

C =\sum_ {i = 1} ^ n cos(a_i)S =\sum_ {i = 1} ^ n sin(a_i)

scipy库通过查找圆形方差

ang = (samples - low)*2.*pi / (high - low)
S = sin(ang).mean(axis=axis)
C = cos(ang).mean(axis=axis)
R = hypot(S, C)
return ((high - low)/2.0/pi)**2 * 2 * log(1/R)

理解这有点棘手.如果我们假设样本是零均值,则范围是[0,2*pi],并且正在使用默认轴(在示例中全部为真),它可以简化为

S = mean(sin(samples))
C = mean(cos(samples))
R = hypot(S, C)
V = 2 * log(1/R)

因此scipy使用的定义将R转换为2*log(1/R),而不是1-R.这看起来很奇怪.纵观历史,https://github.com/scipy/scipy/blame/v1.1.0/scipy/stats/morestats.py#L2696-L2733,一度使用

ang = (samples - low)*2*pi / (high-low)
res = stats.mean(exp(1j*ang))
V = 1-abs(res)
return ((high-low)/2.0/pi)**2 * V

这似乎与您提供的定义一致.在添加测试的同时在错误修正中进行了更改,但没有任何关于新计算来自何处的参考.

有关scipy bug跟踪器的一些讨论,请访问https://github.com/scipy/scipy/pull/5747.它表明这种行为是有意的,不会被修复.有一个在astropy提供另一种实现方式,http://docs.astropy.org/en/stable/api/astropy.stats.circvar.html,其中指出

这里使用的定义与scipy.stats.circvar中的定义不同.准确地说,Scipy circvar使用基于接近线性方差的小角度极限的近似值.

因此,总之,由于未知原因scipy使用近似值(在某些情况下似乎相当差).但是,由于向后兼容性,它将无法修复,因此您可能希望使用astropy的实现.