如何计算滚动/移动平均使用python + NumPy / SciPy?

似乎没有函数可以简单地计算numpy/scipy的移动平均值,从而导致复杂的解决方案

我的问题有两个方面:

  • 用numpy(正确地)实现移动平均的最简单方法是什么?
  • 既然这看起来很重要且容易出错,那么在这种情况下是否有一个好的理由不使用电池包括呢?
348775 次浏览

如果你只想要一个简单的非加权移动平均,你可以很容易地用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已经有了一些专门的功能。

NumPy缺乏特定领域的函数可能是由于核心团队的纪律和对NumPy的主要指令:提供n维数组类型的忠实,以及用于创建和索引这些数组的函数。像许多基本目标一样,这个目标并不小,NumPy出色地完成了它。

更大的< em > SciPy < / em >包含更大的领域特定库集合(被SciPy devs称为< em >子包< / em >)——例如,数值优化(优化),信号处理(信号)和积分(集成)。

我猜你要找的函数至少在SciPy子包中的一个(scipy.signal可能);然而,我将首先在< em > SciPy scikits < / em >的集合中查找,确定相关的scikit(s)并在那里寻找感兴趣的函数。

Scikits是基于NumPy/SciPy独立开发的包,并针对特定的技术规程(例如,scikits-imagescikits-learn等)。其中有几个(特别是用于数值优化的令人敬畏的OpenOpt)早在选择位于相对较新的scikits标题之下之前就被高度重视,成熟的项目。Scikits主页像上面一样列出了大约30个这样的scikits,尽管至少有几个已经不再处于积极的开发中。

按照这个建议,你会得到scikits-timeseries;但是,该软件包已不再处于积极开发阶段;实际上,熊猫<强> < em > < / em > < / >强已经成为,AFAIK,基于事实上的的时间序列库。

熊猫有几个函数可以用来计算移动平均线;其中最简单的可能是rolling_mean,你可以这样使用它:

>>> # 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文档中被非正式地分组在移动窗口 functions的标题下;Pandas中第二组相关的函数被称为指数加权函数(例如,ewma,它计算指数移动加权平均)。第二组不包含在第一个(移动窗口函数)中,这可能是因为指数加权变换不依赖于固定长度的窗口

如果你想仔细处理边缘条件(仅从边缘的可用元素计算平均值),下面的函数就可以了。

import numpy as np


def running_mean(x, N):
out = np.zeros_like(x, dtype=np.float64)
dim_len = x.shape[0]
for i in range(dim_len):
if N%2 == 0:
a, b = i - (N-1)//2, i + (N-1)//2 + 2
else:
a, b = i - (N-1)//2, i + (N-1)//2 + 1


#cap indices to min and max indices
a = max(0, a)
b = min(dim_len, b)
out[i] = np.mean(x[a:b])
return out


>>> running_mean(np.array([1,2,3,4]), 2)
array([1.5, 2.5, 3.5, 4. ])


>>> running_mean(np.array([1,2,3,4]), 3)
array([1.5, 2. , 3. , 3.5])

这个使用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)

现在,只需在dataframe上调用函数rolling,窗口大小,在下面的例子中是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

我觉得这可以很容易地解决使用瓶颈

参见下面的基本示例:

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.convolve。 这背后的思想是利用计算离散卷积的方式,并使用它返回滚的意思。这可以通过与长度等于我们想要的滑动窗口长度的np.ones序列进行卷积来实现

为了做到这一点,我们可以定义以下函数:

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

此函数将对序列x和长度为w的序列进行卷积。注意,所选的modevalid,因此卷积积只对序列完全重叠的点给出。


一些例子:

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.  ])

convolve是如何工作的?

让我们更深入地看看离散卷积是如何计算的。 下面的函数旨在复制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]

所以每一步要做的是取1数组和当前窗口之间的内积。在这种情况下,乘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

依此类推,在所有重叠完成后返回序列的移动平均值。

实际上,我想要一个稍微不同于公认答案的行为。我正在为sklearn管道构建一个移动平均特征提取器,因此我要求移动平均的输出与输入具有相同的维度。我想要的是让移动平均假设级数保持不变,即[1,2,3,4,5]与窗口2的移动平均将得到[1.5,2.5,3.5,4.5,5.0]

对于列向量(我的用例)我们得到

def moving_average_col(X, n):
z2 = np.cumsum(np.pad(X, ((n,0),(0,0)), 'constant', constant_values=0), axis=0)
z1 = np.cumsum(np.pad(X, ((0,n),(0,0)), 'constant', constant_values=X[-1]), axis=0)
return (z1-z2)[(n-1):-1]/n

对于数组

def moving_average_array(X, n):
z2 = np.cumsum(np.pad(X, (n,0), 'constant', constant_values=0))
z1 = np.cumsum(np.pad(X, (0,n), 'constant', constant_values=X[-1]))
return (z1-z2)[(n-1):-1]/n

当然,不必假设填充值为常数,但在大多数情况下这样做应该足够了。

这里有许多实现这一点的方法,以及一些基准测试。最好的方法是使用来自其他库的优化代码。bottleneck.move_mean方法可能是最好的方法。scipy.convolve方法也非常快,可扩展,并且在语法和概念上简单,但是对于非常大的窗口值不能很好地扩展。如果你需要一个纯粹的numpy方法,numpy.cumsum方法是很好的。

注意:其中一些(例如bottleneck.move_mean)不是居中的,并且会转移你的数据。

import numpy as np
import scipy as sci
import scipy.signal as sig
import pandas as pd
import bottleneck as bn
import time as time


def rollavg_direct(a,n):
'Direct "for" loop'
assert n%2==1
b = a*0.0
for i in range(len(a)) :
b[i]=a[max(i-n//2,0):min(i+n//2+1,len(a))].mean()
return b


def rollavg_comprehension(a,n):
'List comprehension'
assert n%2==1
r,N = int(n/2),len(a)
return np.array([a[max(i-r,0):min(i+r+1,N)].mean() for i in range(N)])


def rollavg_convolve(a,n):
'scipy.convolve'
assert n%2==1
return sci.convolve(a,np.ones(n,dtype='float')/n, 'same')[n//2:-n//2+1]


def rollavg_convolve_edges(a,n):
'scipy.convolve, edge handling'
assert n%2==1
return sci.convolve(a,np.ones(n,dtype='float'), 'same')/sci.convolve(np.ones(len(a)),np.ones(n), 'same')


def rollavg_cumsum(a,n):
'numpy.cumsum'
assert n%2==1
cumsum_vec = np.cumsum(np.insert(a, 0, 0))
return (cumsum_vec[n:] - cumsum_vec[:-n]) / n


def rollavg_cumsum_edges(a,n):
'numpy.cumsum, edge handling'
assert n%2==1
N = len(a)
cumsum_vec = np.cumsum(np.insert(np.pad(a,(n-1,n-1),'constant'), 0, 0))
d = np.hstack((np.arange(n//2+1,n),np.ones(N-n)*n,np.arange(n,n//2,-1)))
return (cumsum_vec[n+n//2:-n//2+1] - cumsum_vec[n//2:-n-n//2]) / d


def rollavg_roll(a,n):
'Numpy array rolling'
assert n%2==1
N = len(a)
rolling_idx = np.mod((N-1)*np.arange(n)[:,None] + np.arange(N), N)
return a[rolling_idx].mean(axis=0)[n-1:]


def rollavg_roll_edges(a,n):
# see https://stackoverflow.com/questions/42101082/fast-numpy-roll
'Numpy array rolling, edge handling'
assert n%2==1
a = np.pad(a,(0,n-1-n//2), 'constant')*np.ones(n)[:,None]
m = a.shape[1]
idx = np.mod((m-1)*np.arange(n)[:,None] + np.arange(m), m) # Rolling index
out = a[np.arange(-n//2,n//2)[:,None], idx]
d = np.hstack((np.arange(1,n),np.ones(m-2*n+1+n//2)*n,np.arange(n,n//2,-1)))
return (out.sum(axis=0)/d)[n//2:]


def rollavg_pandas(a,n):
'Pandas rolling average'
return pd.DataFrame(a).rolling(n, center=True, min_periods=1).mean().to_numpy()


def rollavg_bottlneck(a,n):
'bottleneck.move_mean'
return bn.move_mean(a, window=n, min_count=1)


N = 10**6
a = np.random.rand(N)
functions = [rollavg_direct, rollavg_comprehension, rollavg_convolve,
rollavg_convolve_edges, rollavg_cumsum, rollavg_cumsum_edges,
rollavg_pandas, rollavg_bottlneck, rollavg_roll, rollavg_roll_edges]


print('Small window (n=3)')
%load_ext memory_profiler
for f in functions :
print('\n'+f.__doc__+ ' : ')
%timeit b=f(a,3)


print('\nLarge window (n=1001)')
for f in functions[0:-2] :
print('\n'+f.__doc__+ ' : ')
%timeit b=f(a,1001)


print('\nMemory\n')
print('Small window (n=3)')
N = 10**7
a = np.random.rand(N)
%load_ext memory_profiler
for f in functions[2:] :
print('\n'+f.__doc__+ ' : ')
%memit b=f(a,3)


print('\nLarge window (n=1001)')
for f in functions[2:-2] :
print('\n'+f.__doc__+ ' : ')
%memit b=f(a,1001)

定时,小窗口(n=3)

Direct "for" loop :


4.14 s ± 23.7 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)


List comprehension :
3.96 s ± 27.9 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)


scipy.convolve :
1.07 ms ± 26.7 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)


scipy.convolve, edge handling :
4.68 ms ± 9.69 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)


numpy.cumsum :
5.31 ms ± 5.11 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)


numpy.cumsum, edge handling :
8.52 ms ± 11.1 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)


Pandas rolling average :
9.85 ms ± 9.63 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)


bottleneck.move_mean :
1.3 ms ± 12.2 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)


Numpy array rolling :
31.3 ms ± 91.9 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)


Numpy array rolling, edge handling :
61.1 ms ± 55.9 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

定时,大窗口(n=1001)

Direct "for" loop :
4.67 s ± 34 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)


List comprehension :
4.46 s ± 14.6 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)


scipy.convolve :
103 ms ± 165 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)


scipy.convolve, edge handling :
272 ms ± 1.23 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)


numpy.cumsum :
5.19 ms ± 12.4 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)


numpy.cumsum, edge handling :
8.7 ms ± 11.5 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)


Pandas rolling average :
9.67 ms ± 199 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)


bottleneck.move_mean :
1.31 ms ± 15.7 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

内存,小窗口(n=3)

The memory_profiler extension is already loaded. To reload it, use:
%reload_ext memory_profiler


scipy.convolve :
peak memory: 362.66 MiB, increment: 73.61 MiB


scipy.convolve, edge handling :
peak memory: 510.24 MiB, increment: 221.19 MiB


numpy.cumsum :
peak memory: 441.81 MiB, increment: 152.76 MiB


numpy.cumsum, edge handling :
peak memory: 518.14 MiB, increment: 228.84 MiB


Pandas rolling average :
peak memory: 449.34 MiB, increment: 160.02 MiB


bottleneck.move_mean :
peak memory: 374.17 MiB, increment: 75.54 MiB


Numpy array rolling :
peak memory: 661.29 MiB, increment: 362.65 MiB


Numpy array rolling, edge handling :
peak memory: 1111.25 MiB, increment: 812.61 MiB

内存,大窗口(n=1001)

scipy.convolve :
peak memory: 370.62 MiB, increment: 71.83 MiB


scipy.convolve, edge handling :
peak memory: 521.98 MiB, increment: 223.18 MiB


numpy.cumsum :
peak memory: 451.32 MiB, increment: 152.52 MiB


numpy.cumsum, edge handling :
peak memory: 527.51 MiB, increment: 228.71 MiB


Pandas rolling average :
peak memory: 451.25 MiB, increment: 152.50 MiB


bottleneck.move_mean :
peak memory: 374.64 MiB, increment: 75.85 MiB

塔利班成员包含一个简单的移动平均工具,以及其他类似的平均工具(即指数移动平均)。下面将该方法与其他一些解决方案进行比较。


%timeit pd.Series(np.arange(100000)).rolling(3).mean()
2.53 ms ± 40.5 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)


%timeit talib.SMA(real = np.arange(100000.), timeperiod = 3)
348 µs ± 3.5 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)


%timeit moving_average(np.arange(100000))
638 µs ± 45.1 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

需要注意的是,real必须有dtype = float的元素。否则将引发以下错误

例外:实不是双的

下面是一个使用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()))
for i in range(len(Data)):
Data[i, 1] = Data[i-lookback:i, 0].sum() / lookback
试一下这段代码。我认为这样更简单,也能达到目的。 回望是移动平均线的窗口

Data[i-lookback:i, 0].sum()中,我放了0来引用数据集的第一列,但如果你有多个列,你可以放任何你喜欢的列。

移动平均线

迭代器方法

  • 在i处反转数组,简单地取i到n的平均值。

  • 使用列表理解生成动态的迷你数组。

x = np.random.randint(10, size=20)


def moving_average(arr, n):
return [ (arr[:i+1][::-1][:n]).mean() for i, ele in enumerate(arr) ]
d = 5


moving_average(x, d)

张量卷积

moving_average = np.convolve(x, np.ones(d)/d, mode='valid')

我要么使用接受的答案的解决方案,稍微修改为输出与输入的长度相同,要么使用另一个答案的注释中提到的pandas'版本。我在这里用一个可重复的例子来总结两者,以供将来参考:

import numpy as np
import pandas as pd


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


def moving_average_centered(a, n):
return pd.Series(a).rolling(window=n, center=True).mean().to_numpy()


A = [0, 0, 1, 2, 4, 5, 4]
print(moving_average(A, 3))
# [0.         0.         0.33333333 1.         2.33333333 3.66666667 4.33333333]
print(moving_average_centered(A, 3))
# [nan        0.33333333 1.         2.33333333 3.66666667 4.33333333 nan       ]

通过将下面的解决方案与使用cumsum of numpy的解决方案进行比较,这个解决方案几乎需要一半的时间。这是因为它不需要遍历整个数组来做cumsum,然后做所有的减法。此外,如果数组很大且数量很大,则cumsum可以为"危险的" (可能的溢出)。当然,这里也存在危险,但至少我们只把重要的数字加在一起。

def moving_average(array_numbers, n):
if n > len(array_numbers):
return []
temp_sum = sum(array_numbers[:n])
averages = [temp_sum / float(n)]
for first_index, item in enumerate(array_numbers[n:]):
temp_sum += item - array_numbers[first_index]
averages.append(temp_sum / float(n))
return averages

Numpy 1.20开始,sliding_window_view提供了一种在元素窗口中滑动/滚动的方法。然后你可以分别取平均值。

例如,对于4-element窗口:

from numpy.lib.stride_tricks import sliding_window_view


# values = np.array([5, 3, 8, 10, 2, 1, 5, 1, 0, 2])
np.average(sliding_window_view(values, window_shape = 4), axis=1)
# array([6.5, 5.75, 5.25, 4.5, 2.25, 1.75, 2])

注意sliding_window_view的中间结果:

# values = np.array([5, 3, 8, 10, 2, 1, 5, 1, 0, 2])
sliding_window_view(values, window_shape = 4)
# array([[ 5,  3,  8, 10],
#        [ 3,  8, 10,  2],
#        [ 8, 10,  2,  1],
#        [10,  2,  1,  5],
#        [ 2,  1,  5,  1],
#        [ 1,  5,  1,  0],
#        [ 5,  1,  0,  2]])

所有的答案似乎都集中在预先计算的列表的情况下。对于实际运行的用例,数字一个接一个地进来,这里有一个简单的类,它提供了对最后N个值求平均的服务:

import numpy as np
class RunningAverage():
def __init__(self, stack_size):
self.stack = [0 for _ in range(stack_size)]
self.ptr = 0
self.full_cycle = False
def add(self,value):
self.stack[self.ptr] = value
self.ptr += 1
if self.ptr == len(self.stack):
self.full_cycle = True
self.ptr = 0
def get_avg(self):
if self.full_cycle:
return np.mean(self.stack)
else:
return np.mean(self.stack[:self.ptr])

用法:

N = 50  # size of the averaging window
run_avg = RunningAverage(N)
for i in range(1000):
value = <my computation>
run_avg.add(value)
if i % 20 ==0: # print once in 20 iters:
print(f'the average value is {run_avg.get_avg()}')

如果有人需要一个简单的解决方案,这里有一个

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参数来改变窗口之间的重叠

你也可以编写自己的Python C扩展

这当然不是最简单的方法,但是会让你比使用np.cumsum作为构建块更快,更有效地使用内存。

// moving_average.c
#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
#include <Python.h>
#include <numpy/arrayobject.h>


static PyObject *moving_average(PyObject *self, PyObject *args) {
PyObject *input;
int64_t window_size;
PyArg_ParseTuple(args, "Ol", &input, &window_size);
if (PyErr_Occurred()) return NULL;
if (!PyArray_Check(input) || !PyArray_ISNUMBER((PyArrayObject *)input)) {
PyErr_SetString(PyExc_TypeError, "First argument must be a numpy array with numeric dtype");
return NULL;
}
    

int64_t input_size = PyObject_Size(input);
double *input_data;
if (PyArray_AsCArray(&input, &input_data, (npy_intp[]){ [0] = input_size }, 1, PyArray_DescrFromType(NPY_DOUBLE)) != 0) {
PyErr_SetString(PyExc_TypeError, "Failed to simulate C array of type double");
return NULL;
}
    

int64_t output_size = input_size - window_size + 1;
PyObject *output = PyArray_SimpleNew(1, (npy_intp[]){ [0] = output_size }, NPY_DOUBLE);
double *output_data = PyArray_DATA((PyArrayObject *)output);
    

double cumsum_before = 0;
double cumsum_after = 0;
for (int i = 0; i < window_size; ++i) {
cumsum_after += input_data[i];
}
for (int i = 0; i < output_size - 1; ++i) {
output_data[i] = (cumsum_after - cumsum_before) / window_size;
cumsum_after += input_data[i + window_size];
cumsum_before += input_data[i];
}
output_data[output_size - 1] = (cumsum_after - cumsum_before) / window_size;


return output;
}


static PyMethodDef methods[] = {
{
"moving_average",
moving_average,
METH_VARARGS,
"Rolling mean of numpy array with specified window size"
},
{NULL, NULL, 0, NULL}
};


static struct PyModuleDef moduledef = {
PyModuleDef_HEAD_INIT,
"moving_average",
"C extension for finding the rolling mean of a numpy array",
-1,
methods
};


PyMODINIT_FUNC PyInit_moving_average(void) {
PyObject *module = PyModule_Create(&moduledef);
import_array();
return module;
}

构建C扩展

# setup.py
from setuptools import setup, Extension
import numpy


setup(
ext_modules=[
Extension(
'moving_average',
['moving_average.c'],
include_dirs=[numpy.get_include()]
)
]
)


# python setup.py build_ext --build-lib=.

基准

import numpy as np


# Our compiled C extension:
from moving_average import moving_average as moving_average_c


# Answer by Jaime using npcumsum
def moving_average_cumsum(a, n) :
ret = np.cumsum(a, dtype=float)
ret[n:] = ret[n:] - ret[:-n]
return ret[n - 1:] / n


# Answer by yatu using np.convolve
def moving_average_convolve(a, n):
return np.convolve(a, np.ones(n), 'valid') / n


a = np.random.rand(1_000_000)
print('window_size = 3')
%timeit moving_average_c(a, 3)
%timeit moving_average_cumsum(a, 3)
%timeit moving_average_convolve(a, 3)


print('\nwindow_size = 100')
%timeit moving_average_c(a, 100)
%timeit moving_average_cumsum(a, 100)
%timeit moving_average_convolve(a, 100)
window_size = 3
958 µs ± 4.68 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
4.52 ms ± 15.4 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
809 µs ± 463 ns per loop (mean ± std. dev. of 7 runs, 1,000 loops each)


window_size = 100
977 µs ± 937 ns per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
6.16 ms ± 19.1 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
14.2 ms ± 12.4 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

如果你已经有一个已知大小的数组

import numpy as np
M=np.arange(12)
                                                               

avg=[]
i=0
while i<len(M)-2: #for n point average len(M) - (n-1)
avg.append((M[i]+M[i+1]+M[i+2])/3) #n is denominator
i+=1
                                                                                                    

print(avg)