运行或滑动中位数、均值和标准差

Question

我正在尝试计算大型数组的运行中位数、均值和标准差。我知道如何计算运行均值如下：

def running_mean(x, N):
    cumsum = np.cumsum(np.insert(x, 0, 0))
    return (cumsum[N:] - cumsum[:-N]) / float(N)

这非常有效。但是我不太明白为什么(cumsum[N:] - cumsum[:-N]) / float(N)可以给出平均值（我是从别人那里借来的）。

我试图添加另一个 return 句子来计算中位数，但它没有达到我想要的效果。

return (cumsum[N:] - cumsum[:-N]) / float(N), np.median(cumsum[N:] - cumsum[:-N])

有没有人给我一些提示来解决这个问题？非常感谢你。

张华年

Answer 1

那个 cumsum 技巧是特定于查找 sum 或 average 值的，不要认为您可以简单地扩展它以获得 median 和 std 值。在 1D 数组上的 sliding/running window 中执行通用 ufunc 操作的一种方法是创建一系列基于 1D 滑动 windows 的索引堆叠作为二维数组，然后沿堆叠轴应用 ufunc。要获取这些索引，您可以使用 broadcasting.

因此，为了执行运行意思，它看起来像这样 -

idx = np.arange(N) + np.arange(len(x)-N+1)[:,None]
out = np.mean(x[idx],axis=1)

对于运行median和std，只需分别替换np.mean with np.median and np.std即可。

Answer 2

为了估计给定样本集的均值和标准差，存在增量算法 (std, mean)，可帮助您保持较低的计算负荷并进行在线估计。中位数的计算应用排序。您可以近似计算中位数。设 x(t) 是给定时间 t 的数据，m(t) 是时间 t 的中值，m(t-1) 是 e 之前的中值，一个小数，例如e = 0.001 比

m(t) = m(t-1) + e, if m(t-1) < x(t)

m(t) = m(t-1) - e, if m(t-1) > x(t)

m(t) = m(t), else

如果您对中位数 m(0) 有一个很好的初始猜测，那么这很有效。应该根据您的值范围和期望的样本数量来选择 e。例如。如果 x = [-4 2 7.5 2], e = 0.05 会很好，如果 x = [1000 , 3153, -586, -29], e = 10.

Answer 3

让我介绍一个包装器以开始移动"anything":

import numpy as np

def runningFoo(operation):
    """ Make function that applies central running window operation
    """
    assert hasattr(np, operation), f"numpy has no method '{operation}'"

    method = getattr(np, operation)
    assert callable(method), f"numpy.{operation} is not callable"

    def closure(X, windowSize):
        assert windowSize % 2 == 1, "window size must be odd"
        assert windowSize <= len(X), "sequence must be longer than window"

        # setup index matrix
        half = windowSize // 2
        row = np.arange(windowSize) - half
        col = np.arange(len(X))
        index = row + col[:, None]

        # reflect boundaries
        row, col = np.triu_indices(half)
        upper = (row, half - 1 - col)
        index[upper] = np.abs(index[upper]) % len(X)
        lower = (len(X) - 1 - row, windowSize - 1 - upper[1])
        index[lower] = (len(X) - 2 - index[lower]) % len(X)

        return method(X[index], axis=1)

    return closure

例如，如果您希望运行表示您可以调用 runningFoo("mean")。实际上，您可以在 NumPy 中调用任何适当的方法。例如，runningFoo("max") 将是形态学膨胀操作，而 runningFoo("min") 将是形态学腐蚀：

runningStd = runningFoo("std")
runningStd(np.arange(10), windowSize=3)

确保 window 尺寸为奇数。另外，请注意边界点被反映出来。

运行或滑动中位数、均值和标准差

Running or sliding median, mean and standard deviation

python

numpy

mean

median

运行 或滑动中位数、均值和标准差

Running or sliding median, mean and standard deviation

python

numpy

mean

median

运行或滑动中位数、均值和标准差