似乎没有函数可以简单地计算numpy/scipy的移动平均值,这导致了复杂的解决方案。
我的问题有两个方面:
用numpy(正确地)实现移动平均的最简单方法是什么? 既然这似乎不是小事,而且容易出错,有没有一个很好的理由不包括电池在这种情况下?
似乎没有函数可以简单地计算numpy/scipy的移动平均值,这导致了复杂的解决方案。
我的问题有两个方面:
用numpy(正确地)实现移动平均的最简单方法是什么? 既然这似乎不是小事,而且容易出错,有没有一个很好的理由不包括电池在这种情况下?
当前回答
如果有人需要一个简单的解决方案,这里有一个
def moving_average(a,n):
N=len(a)
return np.array([np.mean(a[i:i+n]) for i in np.arange(0,N-n+1)])
你可以通过在np.arange(0,N-n+1,step)中添加step参数来改变窗口之间的重叠
其他回答
下面是一个使用numba的快速实现(注意类型)。注意它确实包含移位的nan。
import numpy as np
import numba as nb
@nb.jit(nb.float64[:](nb.float64[:],nb.int64),
fastmath=True,nopython=True)
def moving_average( array, window ):
ret = np.cumsum(array)
ret[window:] = ret[window:] - ret[:-window]
ma = ret[window - 1:] / window
n = np.empty(window-1); n.fill(np.nan)
return np.concatenate((n.ravel(), ma.ravel()))
我觉得使用瓶颈可以很容易地解决这个问题
参见下面的基本示例:
import numpy as np
import bottleneck as bn
a = np.random.randint(4, 1000, size=(5, 7))
mm = bn.move_mean(a, window=2, min_count=1)
这就给出了每个轴上的移动平均值。
“mm”是“a”的移动平均值。 “窗口”是考虑移动均值的最大条目数。 "min_count"是考虑移动平均值的最小条目数(例如,对于第一个元素或如果数组有nan值)。
好在瓶颈有助于处理nan值,而且非常高效。
实现这一点的一个简单方法是使用np.卷积。 这背后的思想是利用离散卷积的计算方式,并使用它来返回滚动平均值。这可以通过与np序列进行卷积来实现。长度等于我们想要的滑动窗口长度。
为了做到这一点,我们可以定义以下函数:
def moving_average(x, w):
return np.convolve(x, np.ones(w), 'valid') / w
该函数将对序列x和长度为w的序列进行卷积。注意,所选模式是有效的,因此卷积积只对序列完全重叠的点给出。
一些例子:
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。卷积计算输出值:
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]
所以每一步要做的就是求1数组和当前窗口之间的内积。在这种情况下,乘以np.ones(w)是多余的,因为我们直接取序列的和。
下面是一个计算第一个输出的例子,这样会更清楚一些。假设我们想要一个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
依此类推,在所有重叠完成后返回序列的移动平均值。
如果你只想要一个简单的非加权移动平均,你可以很容易地用np实现它。cumsum,可能比基于FFT的方法更快:
修正了Bean在代码中发现的偏离一的错误索引。编辑
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已经有了一些专门的功能。
这个使用Pandas的答案是从上面改编的,因为rolling_mean不再是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)
现在,只需要在窗口大小的数据框架上调用滚动函数,在下面的例子中,窗口大小是10天。
d_mva10 = D.rolling(10).mean()
# d_mva is the same size as the original Series
# though obviously the first w values are NaN where w is the window size
d_mva10[:11]
2010-01-01 NaN
2010-01-02 NaN
2010-01-03 NaN
2010-01-04 NaN
2010-01-05 NaN
2010-01-06 NaN
2010-01-07 NaN
2010-01-08 NaN
2010-01-09 NaN
2010-01-10 4.5
2010-01-11 5.5
Freq: D, dtype: float64