статистика для гистограммы периодических данных

Для ряда значений угла в диапазоне (-pi, pi) я строю гистограмму. Есть ли эффективный способ расчета среднего и модального (вероятностного) значения? Рассмотрим следующие примеры:

import numpy as N, cmath
deg = N.pi/180.
d = N.array([-175., 170, 175, 179, -179])*deg
i = N.sum(N.exp(1j*d))
ave = cmath.phase(i)
i /= float(d.size)
stdev = -2. * N.log(N.sqrt(i.real**2 + i.imag**2))

print ave/deg, stdev/deg

Теперь построим гистограмму:

counts, bins = N.histogram(data, N.linspace(-N.pi, N.pi, 360))

Можно ли рассчитать среднее значение, режим с подсчетами и ячейками? Для непериодических данных расчет среднего значения прост:

ave = sum(counts*bins[:-1])

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

cmax = bins[N.argmax(counts)]
mode = N.mean(N.take(bins, N.nonzero(counts == cmax)[0]))

Я понятия не имею, как рассчитать стандартное отклонение от таких данных. Одним из очевидных решений всех моих проблем (по крайней мере, описанных выше) является преобразование данных гистограммы в ряд данных, а затем использование их в вычислениях. Однако это не элегантно и неэффективно.

Любые подсказки будут очень признательны.

Это частичное решение, которое я написал.

import numpy as N, cmath
import scipy.stats as ST

d = [-175, 170.2, 175.57, 179, -179, 170.2, 175.57, 170.2]
deg = N.pi/180.
data = N.array(d)*deg

i = N.sum(N.exp(1j*data))
ave = cmath.phase(i)  # correct and exact mean for periodic data
wrong_ave = N.mean(d)

i /= float(data.size)
stdev = -2. * N.log(N.sqrt(i.real**2 + i.imag**2))
wrong_stdev = N.std(d)

bins = N.linspace(-N.pi, N.pi, 360)
counts, bins = N.histogram(data, bins, normed=False)
# consider it weighted vector addition
nz = N.nonzero(counts)[0]
weight = counts[nz]
i = N.sum(weight * N.exp(1j*bins[nz])/len(nz))
pave = cmath.phase(i)  # correct and approximated mean for periodic data
i /= sum(weight)/float(len(nz))
pstdev = -2. * N.log(N.sqrt(i.real**2 + i.imag**2))
print
print 'scipy: %12.3f (mean) %12.3f (stdev)' % (ST.circmean(data)/deg, \
                                               ST.circstd(data)/deg)

При запуске дает следующие результаты:

 mean:      175.840       85.843      175.360
stdev:        0.472      151.785        0.430

scipy:      175.840 (mean)        3.673 (stdev)

Теперь несколько комментариев: в первом столбце указано среднее/стандартное отклонение. Как видно, среднее значение хорошо согласуется с scipy.stats.circmean (спасибо JoeKington за указание на это). К сожалению, stdev отличается. Я посмотрю на это позже. Второй столбец дает совершенно неверные результаты (непериодическое среднее/стандартное значение из numpy, очевидно, здесь не работает). В 3-м столбце указано то, что я хотел получить из данных гистограммы (@JoeKington: мои необработанные данные не помещаются в памяти моего компьютера.., @dmytro: спасибо за ваш вклад: конечно, размер бина повлияет на результат, но в моем приложения у меня нет большого выбора, т.е. я должен как-то сократить данные). Как видно, среднее значение (3-й столбец) рассчитано правильно, stdev требует дальнейшего внимания :)


python numpy statistics complex-numbers
person krzym    schedule 22.04.2012    source источник
comment
Если я правильно понял, вы хотите рассчитать среднее значение данных, режим, стандартное отклонение и т. д. из данных гистограммы? Если да, то мне это кажется невозможным, потому что вы теряете много информации, взяв гистограмму данных. Все, что вы можете получить, это приближение, которое становится хуже с более широкими бинами. Или это то, что вы ищете?   -  person dmytro    schedule 22.04.2012
comment
Взгляните на дистрибутив фон Мизеса: en.wikipedia.org/wiki/Von_Mises_distribution . Если вам нужна книга, «Статистический анализ круговых данных» Фишера является стандартным учебником и обычно стоит довольно разумно.   -  person Joe Kington    schedule 22.04.2012

Ответы (2)


arrow_upward
5
arrow_downward

Посмотрите на scipy.stats.circmean и scipy.stats.circstd.

Или у вас есть только гистограмма, а не «сырые» данные? Если это так, вы можете подогнать распределение фон Мизеса к значениям гистограммы и аппроксимировать среднее значение и стандартное отклонение в сюда.

person Joe Kington    schedule 22.04.2012
comment
что, если данные далеки от нормального распределения? - person dmytro; 22.04.2012
comment
@JoeKington: спасибо, что указали на scipy.stats.{circmean,circstd}. Среднее значение, которое я вычислил, точно такое же, как среднее значение. Я изучу код circstd, чтобы выяснить, почему мои результаты отличаются. Я также благодарен за то, что обратил мое внимание на распространение фон Мизеса. Совет по подгонке тоже отличный. На самом деле, прежде чем я придумаю частичное решение (см. редактирование), я самостоятельно придумал аналогичную идею, и она прекрасно работает. - person krzym; 23.04.2012
comment
@dmytro: вы правы, нормальное распределение не является общим решением, в моем случае я установил p[0]*sin(a) exp(-0,5 (a/p[0])**2) с хорошим результатом . Таким образом, подгонка функции к данным гистограммы может быть решением в некоторых случаях. - person krzym; 23.04.2012
comment
@krzym - Для cricstd (и circmean) обязательно правильно укажите kwargs high и low. По умолчанию предполагается, что ваши данные находятся в диапазоне от 0 до 2pi. Если ваши данные указаны в градусах или имеют другой диапазон (например, от -pi до pi), вам необходимо указать другие low и high. - person Joe Kington; 23.04.2012

arrow_upward
1
arrow_downward

Вот как получить приближение.

С Var(x) = <x^2> - <x>^2 имеем:

meanX = N.sum(counts * bins[:-1]) / N.sum(counts)
meanX2 = N.sum(counts * bins[:-1]**2) / N.sum(counts)
std = N.sqrt(meanX2 - meanX**2)
person dmytro    schedule 22.04.2012
comment
Они не применяются к циклическим данным, чего бы это ни стоило. Среднее значение — это не просто среднее :) (например, 359 градусов и 0 градусов отличаются только на 1 градус) - person Joe Kington; 22.04.2012
comment
@JoeKington, достаточно честно. Автор, однако, упомянул непериодические данные и, кажется, совершенно согласен со своим sum(counts * bins[:-1]), поэтому я предположил, что вопрос больше касается оценки моментов по гистограмме. - person dmytro; 22.04.2012
comment
@dmytro: Что я напутал в своем первоначальном вопросе, так это способ расчета среднего значения для непериодических данных (моя исходная гистограмма нормализована, и поэтому я пренебрег делением на сумму отсчетов). На самом деле в моем коде мне нужны оба случая: т.е. мне нужно обрабатывать периодические и непериодические данные, поэтому ваше решение для расчета stdev очень ценится. - person krzym; 23.04.2012