Scipy循环方差

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

请问为什么从函数circvar返回的值有可能大于1？

- Khalil Al Hooti

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使用2*log(1/R)而不是1-R来转换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的实现。

- user2699

Answer 2

根据文档字符串，circvar使用圆形方差的定义，在小角度极限下返回接近于“线性”方差的数字。

......使用圆形方差的定义，在小角度极限下返回接近于“线性”方差的数字。

实际上，它是circstd的平方，维基百科表示：

......取值范围在0到无穷大之间。这个标准偏差的定义......很有用，因为对于一个包裹正态分布，它是潜在正态分布的标准偏差的估计量。因此，它将允许像线性情况一样标准化圆形分布，对于标准偏差的小值。这也适用于von Mises分布......

它还提到，对于小的扩散，两种圆形方差的定义相同，最多相差一个因子二。

- Paul Panzer

Answer 3

也许不应该这样。计算 circstd 看起来很正常：

return ((high - low)/2.0/pi) * sqrt(-2*log(R))

计算circvar的方式似乎不对：

return ((high - low)/2.0/pi)**2 * 2 * log(1/R)

我不知道为什么它要将圆形方差计算为2*ln(1/R)。这可能是我从未见过的一种近似方法，但我不确定——我可能会为此打开一个错误报告。

- CJR

---

这个 var 只是 std 的平方。我并不是说这是否适用于此处。文档字符串说：“这使用了一个圆形方差的定义，在小角度极限下返回一个接近于‘线性’方差的数字。” 这正是维基百科关于 std 的说法 - 所以我认为这看起来没问题。 - Paul Panzer

https://github.com/scipy/scipy/pull/140 是拉取请求，其中方差计算从 1 - R 更改为 2 * log (1/r)，所以那个人可能知道。 - CJR

@CJ59，这个人是numpy和scipy软件包的创始人！ - Khalil Al Hooti

Answer 4

我开发了这段代码，它总是给我一个0-1之间的差异。只是根据我在这里读到的内容进行了适应。

def variance_angle(deg):
    """
    deg: angles in degrees 
    """
    deg = np.deg2rad(deg)
    deg = deg[~np.isnan(deg)]

    S = np.array(deg)
    C = np.array(deg)

    length = C.size

    S = np.sum(np.sin(S))
    C = np.sum(np.cos(C))
    R = np.sqrt(S**2 + C**2)
    R_avg = R/length
    V = 1- R_avg

    return V

- Khalil Al Hooti

Answer 5

这段代码还允许进行加权平均。它返回与适当的维基百科文章中定义的平均值和方差。计算是以弧度为单位的。

def circular_mean(angles, weights=None):
    
    # https://en.wikipedia.org/wiki/Circular_mean
    
    if weights is None:
        weights = np.ones(len(angles))
        
    vectors = [ [w*np.cos(a), w*np.sin(a)]  for a,w in zip(angles,weights) ]
    
    vector = np.sum(vectors, axis=0) / np.sum(weights)
    
    x,y = vector
    
    angle_mean = np.arctan2(y,x)
    angle_variance = 1. - np.linalg.norm(vector)  # x*2+y*2 = hypot(x,y)
    
    return angle_mean, angle_variance

确保权重的总和为正数（不为零），并且所有的权重非负。它们不必被归一化，因为这将在函数中通过除以权重总和来完成。

或者您也可以使用从角度导出的单位向量进行加权平均数的np.average。

此外，请记住存在不同的环形方差公约。

您可以可视化这些向量：

plt.scatter(vectors[:,0], vectors[:,1])
plt.scatter([0],[0],color="black")
plt.grid()    
plt.scatter([x],[y])
plt.show()

- pas-calc