Scipy 循环方差

Question

根据我的理解，循环方差的范围在 0 到 1 之间。这在 wikipedia 和 here 中也得到了证实。但由于某些原因，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 的值

Answer 1

此代码还允许加权平均值。 它返回相应 Wikipedia 文章中定义的均值和方差。 计算以辐射为单位。

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.mean(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

Answer 2

我开发了这段代码，它总是给我 0-1 之间的差异。刚刚改编了我读到的here。

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

Answer 3

这个circvar根据文档字符串

... 使用循环方差的定义，在小范围内 角度返回接近“线性”方差的数字。

其实就是维基百科所说的circstd的平方

... 介于 0 和无穷大之间的值。本标准的定义 偏差 ... 很有用，因为对于包裹的正态分布，它 是基础正态的标准差的估计量 分配。因此，它将允许循环分布 在线性情况下标准化，对于标准的小值 偏差。这也适用于 von Mises 分布...

它还提到，对于小散布，循环方差的两个定义是相同的，最多可达两倍。

Answer 4

不太有用的答案是因为 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]，并且使用默认轴（在示例中均为 true），则可以将其简化为

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

这似乎符合您提供的定义。 在添加测试的同时在错误修复中对此进行了更改，但没有提及新计算的来源。

https://github.com/scipy/scipy/pull/5747 上提供了有关 scipy 错误跟踪器的一些讨论。 它表明该行为是故意的，不会被修复。 astropy 中还有另一种实现，http://docs.astropy.org/en/stable/api/astropy.stats.circvar.html，它指出

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

因此，总而言之，由于未知原因scipy 使用了近似值（在某些情况下似乎相当差）。但是，由于向后兼容，它不会被修复，因此您可能希望使用astropy 的实现。

Answer 5

这可能不应该。 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)。这可能是我以前从未见过的近似值，但我不知道 - 我可能会为此打开一个错误。