周期数据直方图统计

Question

对于(-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)

现在有几个注释:第一列给出了计算平均值/stdev。可以看出，这个平均值与scipy.stats.circmean非常一致(谢谢JoeKington指出这一点)。不幸的是，stdev不同。我稍后再看。第二列给出了完全错误的结果( numpy中的非周期均值/std显然在这里不起作用)。第三栏给出了一些我想从直方图数据中获得的东西(@JoeKington:我的原始数据不适合我的计算机内存。) @dmytro:谢谢您的输入:当然，垃圾箱大小会影响结果，但是在我的应用程序中，我没有太多的选择，即我不得不以某种方式减少数据。可以看出，平均值(第3列)是正确计算的，stdev需要进一步注意:)

Answer 1

看看scipy.stats.circmean和scipy.stats.circstd。

还是只有直方图计数，而不是“原始”数据？如果是这样的话，您可以将Von分布与直方图计数相匹配，并以这种方式近似平均值和stddev。

Answer 2

下面是如何得到近似的方法。

自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)