4

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

但是我现在使用的数组足够大，以至于速度非常慢。有谁知道提供这种估算器的软件包，或者对如何以快速有效的方式解决这个问题有任何建议？

4

1 回答 1

3

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

从概念上讲，您可以通过 获得均值的移动估计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() 参数需要针对您的问题进行调整。如果您将 order 设置为 1 而不是 2，您将得到我首先描述的简单过滤器。

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

于 2013-04-14T18:21:09.497 回答