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

Posted 2023-02-16

【中文标题】运行或滑动中位数、平均值和标准差

【英文标题】：Running or sliding median, mean and standard deviation

【发布时间】：2016-02-08 17:00:56

【问题描述】：

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

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 (cumsum[N:] - cumsum[:-N]) / float(N), np.median(cumsum[N:] - cumsum[:-N])

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

张华年

【问题讨论】：

您的目标是了解正在发生的事情，还是只想使用这些功能？在后一种情况下，您可以在pandas 中找到许多预定义的 - 请参阅pandas.pydata.org/pandas-docs/stable/…

median 是一个奇怪的统计数据。它只是中间值，或 2 个中间值的平均值。这需要对值进行排序，这与对它们求和完全不同。

【参考方案1】：

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

因此，对于执行运行均值，它看起来像这样 -

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

对于运行median 和std，只需将np.mean 分别替换为np.median 和np.std。

【讨论】：

我试过你的方法，尽管我不太明白你是如何进行滑动/运行的。但它给出了一些错误。我把这两句话放在我的running_mean的定义下，然后返回，但它说“TypeError：只有一个元素的整数数组可以转换为索引”。你能给我解释一下吗？非常感谢。

@HuanianZhang 让我问你 - x 是一维 NumPy 数组，N 是标量吗？

可能我的 x 只是一个列表，我需要转移到一维 numpy 数组。非常感谢。

这两行代码非常适合小型 np.array。但它总是为一个巨大的数组返回一些错误。我正在处理超过 100 万个元素的数组。我正在尝试修改您的代码以避免此类错误，但我没有成功。你能给我一些想法吗？非常感谢mcuh~

非常感谢。我已经修好了。看起来我有一些错字。对此感到抱歉。

【参考方案2】：

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

m(t) = m(t-1) + e，如果 m(t-1)

m(t) = m(t-1) - e，如果 m(t-1) > x(t)

m(t) = m(t)，否则

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

【参考方案3】：

让我介绍一个包装器来移动“任何东西”：

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)

确保窗口大小是奇数。另外，请注意边界点会被反射。