python - python数组的高效移动、稳健尺度估计

Question

我正在寻找一种快速有效的方法来计算一组数据的稳健的移动尺度估计。我正在处理通常包含 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))

然而，我现在使用的数组足够大，速度太慢了。有谁知道提供此类估算器的软件包，或者对如何以快速有效的方式处理此问题有任何建议？

Answer 1

我在类似情况下使用了一个简单的低通滤波器。

从概念上讲，您可以使用 fac = 0.99; 获得均值的移动估计值； filtered[k] = fac*filtered[k-1] + (1-fac)*data[k]，实现起来非常高效(在 C 中)。一个比这个稍微更花哨的 IIR 滤波器，巴特沃斯低通，很容易在 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，您将获得我首先描述的简单过滤器。

免责声明:这是工程师对问题的看法。这种方法在任何统计或数学方面可能都不可靠。另外，我不确定它是否真的解决了你的问题(如果没有，请更好地解释)，但别担心，无论哪种方式，我都玩得很开心 ;-)