我正在寻找一种快速有效的方法来计算一组数据的稳健、移动的规模估计。我正在处理通常为 3-400k 元素的一维数组。直到最近,我一直在处理模拟数据(没有灾难性的异常值),出色的 Bottleneck 包中的 move_std 函数对我很有用。但是,随着我过渡到嘈杂的数据,std 不再表现得足够好而有用。
过去,我使用一个非常简单的双权重中方差代码逐个元素来处理行为不良的分布问题:
def bwmv(data_array):
cent = np.median(data_array)
MAD = np.median(np.abs(data_array-cent))
u = (data_array-cent) / 9. / MAD
uu = u*u
I = np.asarray((uu <= 1.), dtype=int)
return np.sqrt(len(data_array) * np.sum((data_array-cent)**2 * (1.-uu)**4 * I)\
/(np.sum((1.-uu) * (1.-5*uu) * I)**2))
但是我现在正在使用的数组足够大,这太慢了。有谁知道提供这种估算器的软件包,或者对如何以快速有效的方式解决这个问题有任何建议?
最佳答案
我在类似情况下使用了一个简单的低通滤波器。
从概念上讲,您可以使用 fac = 0.99; filtered[k] = fac*filtered[k-1] + (1-fac)*data[k]
获得均值的移动估计,这在实现上非常有效(在 C 中)。一个比这个更花哨的 IIR 滤波器,butterworth 低通,很容易在 scipy 中设置:
b, a = scipy.signal.butter(2, 0.1)
filtered = scipy.signal.lfilter(b, a, data)
要获得“比例”的估计值,您可以从数据中减去这个“平均估计值”。这实际上将低通变成了高通滤波器。获取 abs() 并通过另一个低通滤波器运行它。
结果可能如下所示:
完整脚本:
from pylab import *
from scipy.signal import lfilter, butter
data = randn(1000)
data[300:] += 1.0
data[600:] *= 3.0
b, a = butter(2, 0.03)
mean_estimate = lfilter(b, a, data)
scale_estimate = lfilter(b, a, abs(data-mean_estimate))
plot(data, '.')
plot(mean_estimate)
plot(mean_estimate + scale_estimate, color='k')
plot(mean_estimate - scale_estimate, color='k')
show()
显然,butter() 参数需要根据您的问题进行调整。如果您将顺序设置为 1 而不是 2,您将获得我首先描述的简单过滤器。
免责声明:这是工程师对问题的看法。这种方法可能在任何统计或数学方面都不合理。另外,我不确定它是否真的解决了您的问题(如果没有,请更好地解释),但别担心,无论哪种方式,我都玩得很开心;-)
关于python - python数组的有效移动,稳健的尺度估计,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/16000637/