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

我的问题有两个方面:

用numpy(正确地)实现移动平均的最简单方法是什么? 既然这似乎不是小事,而且容易出错,有没有一个很好的理由不包括电池在这种情况下?


当前回答

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

其他回答

通过比较下面的解决方案与使用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

所有的答案似乎都集中在预先计算的列表的情况下。对于实际运行的用例,数字一个接一个地进来,这里有一个简单的类,它提供了对最后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()}')

实际上,我想要一个稍微不同于公认答案的行为。我正在为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

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

实现这一点的一个简单方法是使用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

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

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的元素。否则将引发以下错误

例外:实不是双的