如何用NumPy计算移动平均值？

Question

似乎没有函数可以简单地计算numpy/scipy上的移动平均值，从而得到convoluted solutions。

我的问题有两个：

• 用numpy（正确地）实现移动平均线最简单的方法是什么？
• 既然这看起来不简单而且容易出错，那么在这种情况下，有没有一个好的理由不使用batteries included？

Answer 1

实现这一点的一个简单方法是使用^{}。 这背后的想法是利用计算discrete convolution的方式，并使用它返回滚动平均值。这可以通过卷积一个长度等于我们想要的滑动窗口长度的^{}序列来实现。

为此，我们可以定义以下函数：

def moving_average(x, w):
    return np.convolve(x, np.ones(w), 'valid') / w

这个函数将取序列x和长度为w的序列的卷积。注意，所选择的mode是valid，因此卷积积只针对序列完全重叠的点给出。

---

用例

一些例子：

x = np.array([5,3,8,10,2,1,5,1,0,2])

对于长度为2的移动平均线，我们将得到：

moving_average(x, 2)
# array([4. , 5.5, 9. , 6. , 1.5, 3. , 3. , 0.5, 1. ])

对于长度为4的窗口：

moving_average(x, 4)
# array([6.5 , 5.75, 5.25, 4.5 , 2.25, 1.75, 2.  ])

---

详细信息

让我们更深入地了解一下离散卷积的计算方式。 以下函数旨在复制np.convolve计算输出值的方式：

def mov_avg(x, w):
    for m in range(len(x)-(w-1)):
        yield sum(np.ones(w) * x[m:m+w]) / w 

对于上述相同的例子，这也会产生：

list(mov_avg(x, 2))
# [4.0, 5.5, 9.0, 6.0, 1.5, 3.0, 3.0, 0.5, 1.0]

所以每一步都要做的是，取一个数组和当前窗口之间的内积。在这种情况下，乘以np.ones(w)是多余的，因为我们直接取序列的sum。

下面是如何计算第一个输出的一个例子，这样就更清楚了。假设我们想要一个w=4的窗口：

[1,1,1,1]
[5,3,8,10,2,1,5,1,0,2]
= (1*5 + 1*3 + 1*8 + 1*10) / w = 6.5

以下输出将计算为：

  [1,1,1,1]
[5,3,8,10,2,1,5,1,0,2]
= (1*3 + 1*8 + 1*10 + 1*2) / w = 5.75

依此类推，在执行所有重叠之后返回序列的移动平均值。

Answer 2

如果您只需要一个简单的非加权移动平均值，您可以使用np.cumsum轻松实现它，它可能比基于FFT的方法快：

编辑更正了代码中Bean发现的一个off-by-one错误索引。编辑

def moving_average(a, n=3) :
    ret = np.cumsum(a, dtype=float)
    ret[n:] = ret[n:] - ret[:-n]
    return ret[n - 1:] / n

>>> a = np.arange(20)
>>> moving_average(a)
array([  1.,   2.,   3.,   4.,   5.,   6.,   7.,   8.,   9.,  10.,  11.,
        12.,  13.,  14.,  15.,  16.,  17.,  18.])
>>> moving_average(a, n=4)
array([  1.5,   2.5,   3.5,   4.5,   5.5,   6.5,   7.5,   8.5,   9.5,
        10.5,  11.5,  12.5,  13.5,  14.5,  15.5,  16.5,  17.5])

所以我想答案是：实现起来确实很容易，也许numpy已经有点臃肿了，有专门的功能。

Answer 3

NumPy缺少特定于域的函数可能是由于核心团队的纪律性和对NumPy主指令的忠实性：提供了N维数组类型，以及创建和索引这些数组的函数。像许多基本目标一样，这个目标并不小，纽比做得很出色。

更大的SciPy包含更大的领域特定库集合（SciPy devs称为子包）——例如，数值优化（优化）、信号处理（信号）和积分（积分）。

我的猜测是，您所追求的功能至少在一个SciPy子包中（SciPy.signal也许）；但是，我将首先在SciPy scikit s集合中查找，确定相关的scikit并在那里查找感兴趣的功能。

Scikits是基于NumPy/SciPy独立开发的软件包，面向特定的技术领域（例如，scikits image，scikits learn，等等），其中一些软件（特别是用于数值优化的令人敬畏的OpenOpt）受到高度重视，在选择居住在相对较新的scikits标准下之前很久就已经成熟的项目。上面喜欢的Scikits主页列出了大约30个这样的Scikits，尽管其中至少有几个已经不在积极开发之中。

遵循此建议将引导您找到scikits time series；但是，该包不再处于活动开发中；实际上，Pandas已经成为基于事实上的NumPy的时间序列库。

Pandas有几个函数可用于计算移动平均值；其中最简单的函数可能是滚动平均值，您可以这样使用：

>>> # the recommended syntax to import pandas
>>> import pandas as PD
>>> import numpy as NP

>>> # prepare some fake data:
>>> # the date-time indices:
>>> t = PD.date_range('1/1/2010', '12/31/2012', freq='D')

>>> # the data:
>>> x = NP.arange(0, t.shape[0])

>>> # combine the data & index into a Pandas 'Series' object
>>> D = PD.Series(x, t)

现在，只需调用函数rolling_mean传入Series对象和窗口大小，在下面的示例中是10天。

>>> d_mva = PD.rolling_mean(D, 10)

>>> # d_mva is the same size as the original Series
>>> d_mva.shape
    (1096,)

>>> # though obviously the first w values are NaN where w is the window size
>>> d_mva[:3]
    2010-01-01         NaN
    2010-01-02         NaN
    2010-01-03         NaN

验证它是否有效——例如，将原始序列中的值10-15与使用滚动平均值平滑的新序列进行比较

>>> D[10:15]
     2010-01-11    2.041076
     2010-01-12    2.041076
     2010-01-13    2.720585
     2010-01-14    2.720585
     2010-01-15    3.656987
     Freq: D

>>> d_mva[10:20]
      2010-01-11    3.131125
      2010-01-12    3.035232
      2010-01-13    2.923144
      2010-01-14    2.811055
      2010-01-15    2.785824
      Freq: D

函数rolling_mean，以及大约十几个其他函数在Pandas文档中非正式地分组为moving window函数；Pandas中的第二个相关函数组称为指数加权函数（例如，ewma，计算指数移动加权平均值）。第二组不包含在第一组（移动窗口函数）中，这可能是因为指数加权的变换不依赖于固定长度的窗口