Scipy Circular Variance

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

Это circvar в соответствии с строкой документации

... использует определение круговой дисперсии, которое в пределе малых углов возвращает число, близкое к «линейной» дисперсии.

На самом деле, это квадрат circstd, из которого википедия сообщает

... значения от 0 до бесконечности.Это определение стандартного отклонения ... полезно, потому что для упакованного нормального распределения оно является оценкой стандартного отклонения базового нормального распределения.Поэтому это позволит стандартизировать круговое распределение, как в линейном случае, для малых значений стандартного отклонения.Это также относится к распределению фон Мизеса ...

также упоминается, что при небольшом разбросе два определения круговой дисперсии одинаковы до двух раз.

Answer 2

Менее полезный ответ будет, поскольку именно так его определяет 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

, что соответствует приведенным вами определениям.Это было исправлено в багфиксе в то же время, когда были добавлены тесты, но без какой-либо ссылки на то, откуда появились новые вычисления.

Некоторое обсуждение трескотного баг-трекера доступно по адресу https://github.com/scipy/scipy/pull/5747. Этопредполагает, что поведение является преднамеренным, и не будет исправлено.Есть еще одна реализация, доступная в astropy, http://docs.astropy.org/en/stable/api/astropy.stats.circvar.html,, которая отмечает

Определение, используемое здесь, отличается от определения в scipy.stats.circvar.Точно, Scipy Cirvar использует приближение, основанное на пределе малых углов, которое приближается к линейной дисперсии.

Итак, в итоге, по неизвестным причинам scipy использует приближение (которое кажется довольно плохимв некоторых случаях).Однако из-за обратной совместимости это не будет исправлено, поэтому вы можете использовать реализацию astropy.

Answer 3

Я разработал этот код, и он всегда дает мне разницу между 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

Answer 4

Наверное, не должно быть.Расчет для 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).Это может быть приближением, которого я никогда раньше не видел, но я не знаю - вероятно, я бы открыл для этого ошибку.