似乎没有函数可以简单地计算numpy/scipy的移动平均值,这导致了复杂的解决方案。
我的问题有两个方面:
用numpy(正确地)实现移动平均的最简单方法是什么? 既然这似乎不是小事,而且容易出错,有没有一个很好的理由不包括电池在这种情况下?
似乎没有函数可以简单地计算numpy/scipy的移动平均值,这导致了复杂的解决方案。
我的问题有两个方面:
用numpy(正确地)实现移动平均的最简单方法是什么? 既然这似乎不是小事,而且容易出错,有没有一个很好的理由不包括电池在这种情况下?
当前回答
实现这一点的一个简单方法是使用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.卷积。 这背后的思想是利用离散卷积的计算方式,并使用它来返回滚动平均值。这可以通过与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
依此类推,在所有重叠完成后返回序列的移动平均值。
实际上,我想要一个稍微不同于公认答案的行为。我正在为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
当然,不必假设填充值为常数,但在大多数情况下这样做应该足够了。
从Numpy 1.20开始,sliding_window_view提供了一种在元素窗口中滑动/滚动的方法。然后你可以分别取平均值。
例如,对于一个4元素的窗口:
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]])
下面是一个使用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()))
Talib包含一个简单的移动平均工具,以及其他类似的平均工具(即指数移动平均)。下面将该方法与其他一些解决方案进行比较。
%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的元素。否则将引发以下错误
例外:实不是双的