似乎没有函数可以简单地计算numpy/scipy的移动平均值,这导致了复杂的解决方案。
我的问题有两个方面:
用numpy(正确地)实现移动平均的最简单方法是什么? 既然这似乎不是小事,而且容易出错,有没有一个很好的理由不包括电池在这种情况下?
似乎没有函数可以简单地计算numpy/scipy的移动平均值,这导致了复杂的解决方案。
我的问题有两个方面:
用numpy(正确地)实现移动平均的最简单方法是什么? 既然这似乎不是小事,而且容易出错,有没有一个很好的理由不包括电池在这种情况下?
当前回答
我觉得使用瓶颈可以很容易地解决这个问题
参见下面的基本示例:
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值,而且非常高效。
其他回答
所有的答案似乎都集中在预先计算的列表的情况下。对于实际运行的用例,数字一个接一个地进来,这里有一个简单的类,它提供了对最后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()}')
如果你想仔细考虑边缘条件(只从边缘的可用元素计算平均值),下面的函数可以解决这个问题。
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])
NumPy缺乏特定领域的函数可能是由于核心团队的纪律和对NumPy主要指令的忠实:提供n维数组类型,以及用于创建和索引这些数组的函数。像许多基本目标一样,这个目标并不小,NumPy出色地完成了它。
更大的SciPy包含更大的特定于领域的库集合(被SciPy开发人员称为子包)——例如,数值优化(optimize)、信号处理(signal)和积分(integrate)。
我的猜测是,您要查找的函数至少在SciPy子包中的一个(SciPy。也许信号);然而,我将首先在SciPy scikit集合中查找,确定相关的scikit并在其中寻找感兴趣的函数。
Scikits是基于NumPy/SciPy独立开发的包,并针对特定的技术规程(例如,Scikits -image, Scikits -learn等),其中几个(特别是用于数值优化的令人钦佩的OpenOpt)在选择位于相对较新的Scikits主题之下很久以前就得到了高度重视,成熟的项目。Scikits主页上列出了大约30个这样的Scikits,尽管其中至少有几个已经不再处于积极的开发中。
按照这个建议,你会发现scikits-timeseries;但是,该软件包已不再处于积极开发阶段;实际上,Pandas已经成为AFAIK,事实上的基于numpy的时间序列库。
Pandas有几个函数可以用来计算移动平均线;其中最简单的可能是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
The function rolling_mean, along with about a dozen or so other function are informally grouped in the Pandas documentation under the rubric moving window functions; a second, related group of functions in Pandas is referred to as exponentially-weighted functions (e.g., ewma, which calculates exponentially moving weighted average). The fact that this second group is not included in the first (moving window functions) is perhaps because the exponentially-weighted transforms don't rely on a fixed-length window
这里有许多实现这一点的方法,以及一些基准测试。最好的方法是使用来自其他库的优化代码。瓶颈。Move_mean方法可能是最好的方法。scipy。卷积方法也非常快,可扩展,并且语法和概念简单,但是对于非常大的窗口值不能很好地扩展。numpy。如果你需要一个纯numpy方法,Cumsum方法是很好的。
注意:其中一些(例如:瓶颈。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
这个使用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