高效移动,健壮的python数组估计

我正在寻找一种快速有效的方法来计算一组数据的稳健的移动尺度估计.我正在使用通常3-400k元素的1d阵列.直到最近,我一直在使用模拟数据(没有灾难性的异常值),而优秀的Bottleneck软件包中的move_std函数对我很有帮助.但是,当我转换为噪声数据时,std不再表现得非常有用.

在过去,我使用了一个非常简单的biweight mid-variance代码元素来处理行为不良的问题：

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; filter [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,那么您将获得我首先描述的简单过滤器.

免责声明：这是工程师对此问题的看法.这种方法可能在任何统计或数学方面都不合理.另外,我不确定它是否真的解决了你的问题(如果没有,请更好地解释),但不要担心,无论如何我都有一些乐趣;-)