Moving average or running mean

Существует ли функция SciPy или функция NumPy или модуль для Python, который вычисляет текущее среднее значение одномерного массива для определенного окна?

^{ПРИМЕЧАНИЕ: Более эффективные решения могут включать scipy.ndimage.uniform_filter1d (см. Этот ответ) или использовать более новые библиотеки, включая talib talib.MA.}

Использование np.convolve:

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

Объяснение

Текущее среднее - это случай математической операции свертки. Для текущего среднего вы перемещаете окно вдоль входных данных и вычисляете среднее значение содержимого окна. Для дискретных одномерных сигналов свертка - это то же самое, за исключением того, что вместо среднего вы вычисляете произвольную линейную комбинацию, т. Е. Умножаете каждый элемент на соответствующий коэффициент и суммируете результаты. Эти коэффициенты, по одному для каждой позиции в окне, иногда называют ядром свертки. Среднее арифметическое из N значений равно (x_1 + x_2 + ... + x_N) / N, поэтому соответствующее ядро равно (1/N, 1/N, ..., 1/N), и это именно то, что мы получаем, используя np.ones(N)/N.

Ребра

В mode аргументе np.convolve указывается, как обрабатывать ребра. Я выбрал здесь режим valid, потому что думаю, что большинство людей ожидают, что именно так будет работать running mean, но у вас могут быть другие приоритеты. Вот график, иллюстрирующий разницу между режимами:

import numpy as np
import matplotlib.pyplot as plt
modes = ['full', 'same', 'valid']
for m in modes:
    plt.plot(np.convolve(np.ones(200), np.ones(50)/50, mode=m));
plt.axis([-10, 251, -.1, 1.1]);
plt.legend(modes, loc='lower center');
plt.show()

Эффективное решение

Свертка намного лучше, чем простой подход, но (я думаю) она использует БПФ и, следовательно, довольно медленная. Однако специально для вычисления текущего среднего следующий подход работает нормально

def running_mean(x, N):
    cumsum = numpy.cumsum(numpy.insert(x, 0, 0)) 
    return (cumsum[N:] - cumsum[:-N]) / float(N)

Код для проверки

In[3]: x = numpy.random.random(100000)
In[4]: N = 1000
In[5]: %timeit result1 = numpy.convolve(x, numpy.ones((N,))/N, mode='valid')
10 loops, best of 3: 41.4 ms per loop
In[6]: %timeit result2 = running_mean(x, N)
1000 loops, best of 3: 1.04 ms per loop

Обратите внимание, что numpy.allclose(result1, result2) есть True, два метода эквивалентны.
Чем больше N, тем больше разница во времени.

предупреждение: хотя cumsum работает быстрее, будет увеличена ошибка с плавающей запятой, которая может привести к тому, что ваши результаты будут недействительными / некорректными / неприемлемыми

в комментариях здесь указывалась на эту проблему с ошибкой с плавающей запятой, но я делаю это более очевидным здесь, в ответе..

# demonstrate loss of precision with only 100,000 points
np.random.seed(42)
x = np.random.randn(100000)+1e6
y1 = running_mean_convolve(x, 10)
y2 = running_mean_cumsum(x, 10)
assert np.allclose(y1, y2, rtol=1e-12, atol=0)

• чем больше очков вы накапливаете, тем больше ошибка с плавающей запятой (так что 1e5 баллов заметны, 1e6 баллов более значимы, больше, чем 1e6, и вы можете захотеть сбросить накопители)


• вы можете схитрить, используя np.longdouble но ваша ошибка с плавающей запятой по-прежнему будет значимой для относительно большого количества точек (около > 1e5, но зависит от ваших данных)


• вы можете отобразить ошибку на графике и увидеть, как она относительно быстро увеличивается


• решение convolve работает медленнее, но не приводит к потере точности с плавающей запятой


• решение uniform_filter1d быстрее, чем это решение cumsum, И не имеет такой потери точности с плавающей запятой

Вы можете использовать scipy.ndimage.uniform_filter1d:

import numpy as np
from scipy.ndimage import uniform_filter1d
N = 1000
x = np.random.random(100000)
y = uniform_filter1d(x, size=N)

uniform_filter1d:

• выдает результат с той же числовой формой (т. Е. Количеством точек)


• позволяет несколькими способами обрабатывать границу, где по умолчанию используется 'reflect', но в моем случае я скорее хотел 'nearest'

Это также довольно быстро (почти в 50 раз быстрее, чем np.convolve и в 2-5 раз быстрее, чем подход cumsum, приведенный выше):

%timeit y1 = np.convolve(x, np.ones((N,))/N, mode='same')
100 loops, best of 3: 9.28 ms per loop

%timeit y2 = uniform_filter1d(x, size=N)
10000 loops, best of 3: 191 µs per loop

вот 3 функции, которые позволяют сравнивать ошибки / скорость различных реализаций:

from __future__ import division
import numpy as np
import scipy.ndimage as ndi
def running_mean_convolve(x, N):
    return np.convolve(x, np.ones(N) / float(N), 'valid')
def running_mean_cumsum(x, N):
    cumsum = np.cumsum(np.insert(x, 0, 0))
    return (cumsum[N:] - cumsum[:-N]) / float(N)
def running_mean_uniform_filter1d(x, N):
    return ndi.uniform_filter1d(x, N, mode='constant', origin=-(N//2))[:-(N-1)]

Обновление: Приведенный ниже пример показывает старую pandas.rolling_mean функцию, которая была удалена в последних версиях pandas. Современный эквивалент вызова этой функции использовал бы pandas.Серия.переходящая:

In [8]: pd.Series(x).rolling(window=N).mean().iloc[N-1:].values
Out[8]: 
array([ 0.49815397,  0.49844183,  0.49840518, ...,  0.49488191,
        0.49456679,  0.49427121])

Для этого больше подходитpandas, чем NumPy или SciPy. Его функция rolling_mean выполняет эту работу удобно. Он также возвращает массив NumPy, когда входными данными является массив.

Трудно превзойти rolling_mean по производительности любую пользовательскую реализацию на чистом Python. Вот пример производительности по сравнению с двумя предлагаемыми решениями:

In [1]: import numpy as np

In [2]: import pandas as pd

In [3]: def running_mean(x, N):
   ...:     cumsum = np.cumsum(np.insert(x, 0, 0)) 
   ...:     return (cumsum[N:] - cumsum[:-N]) / N
   ...:

In [4]: x = np.random.random(100000)

In [5]: N = 1000

In [6]: %timeit np.convolve(x, np.ones((N,))/N, mode='valid')
10 loops, best of 3: 172 ms per loop

In [7]: %timeit running_mean(x, N)
100 loops, best of 3: 6.72 ms per loop

In [8]: %timeit pd.rolling_mean(x, N)[N-1:]
100 loops, best of 3: 4.74 ms per loop

In [9]: np.allclose(pd.rolling_mean(x, N)[N-1:], running_mean(x, N))
Out[9]: True

Есть также хорошие варианты того, как обращаться со значениями границ.