Кажется, что нет функции, которая просто вычисляет скользящее среднее на numpy/scipy, что приводит к запутанным решениям.
Мой вопрос в два раза:
Если вам просто нужна простая взвешенная скользящая средняя, вы можете легко реализовать ее с помощью 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 prime: предоставить N-мерный тип массива, а также функции для создания и индексирования этих массивов. Как и многие основополагающие цели, этот не мал, и NumPy делает это блестяще.
(намного больше) SciPy содержит гораздо больший набор библиотек, специфичных для домена (называемых подпакетами для разработчиков SciPy) - например, численная оптимизация (оптимизация), процесс обработки сигналов (сигнал) и интегральное исчисление (интеграция).
Я предполагаю, что функция, по которой вы находитесь, находится, по крайней мере, в одном из подпакетов SciPy (возможно, scipy.signal); однако я бы посмотрел сначала в коллекции SciPy scikits, идентифицировал соответствующие scikit и искал интересующую функцию.
Scikits - это независимо разработанные пакеты на основе NumPy/SciPy и направленные на определенную техническую дисциплину (например, scikits-image, scikits- узнать и т.д.). Некоторые из них (в частности, удивительный OpenOpt для численной оптимизации) были высоко оценены, зрелые проекты задолго до выбора проживать под относительно новой рубрикой scikits. На домашней странице Scikits было написано около 30 таких сценариев, хотя по крайней мере некоторые из них больше не находятся в активной разработке.
Следуя этому совету, вы попадете в scikits-timeseries; однако этот пакет больше не находится в активной разработке; Фактически, Pandas стал AFAIK, фактической библиотекой временных рядов на основе NumPy.
Pandas имеет несколько функций, которые могут быть использованы для вычисления скользящей средней; самый простой из них, вероятно, roll_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)
Теперь просто вызовите функцию roll_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
Функция roll_mean, а также около дюжины других функций неформально сгруппированы в документацию Pandas в разделе функций движущихся окон рубрики; вторая связанная группа функций в Pandas называется экспоненциально взвешенными функциями (например, ewma, которая вычисляет экспоненциально движущуюся средневзвешенную). Тот факт, что эта вторая группа не включена в первую (функции движущегося окна), возможно, потому, что экспоненциально взвешенные преобразования не полагаются на окно фиксированной длины
Простой способ добиться этого - использовать np.convolve. Идея заключается в том, чтобы использовать способ вычисления дискретной свертки и использовать его для возврата скользящего среднего. Это можно сделать, свернув с последовательностью np.ones длины, равной длине скользящего окна, которую мы хотим.
Для этого мы можем определить следующую функцию:
def moving_average(x, w):
return np.convolve(x, np.ones(w), 'valid') / w
Эта функция будет принимать свертку последовательности x и последовательности единиц длины w. Обратите внимание, что выбранный mode valid так что произведение свертки дается только для точек, где последовательности полностью перекрываются.
Случай использования
Некоторые примеры:
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.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]
Итак, что делается на каждом шаге - это поместить внутренний продукт между массивом единиц и текущим окном. В этом случае умножение на 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
И так далее, возвращая скользящее среднее последовательности после того, как все перекрытия были выполнены.
Этот ответ с использованием панд адаптирован сверху, так как 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 по информационному фрейму с размером окна, который в моем примере ниже равен 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)
Это дает среднее значение перемещения по каждой оси.
"мм" - это скользящее среднее для "а".
"Окно" - это максимальное количество записей, которые нужно учитывать для скользящего среднего.
"min_count" - это минимальное количество записей, которые нужно учитывать для скользящего среднего (например, для первого элемента или если массив имеет значения nan).
Хорошая часть заключается в том, что "Узкое место" помогает справиться со значениями наночастиц, а также очень эффективно.
Если вы хотите тщательно позаботиться о краевых условиях (вычислите среднее значение только из доступных элементов по краям), следующая функция поможет.
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])
Вот несколько способов сделать это, а также некоторые тесты. Лучшие методы - это версии, использующие оптимизированный код из других библиотек. Метод bottleneck.move_mean, вероятно, лучший во всем. Подход scipy.convolve также очень быстрый, расширяемый, синтаксически и концептуально простой, но не подходит для очень больших значений окна. Метод numpy.cumsum хорош, если вам нужен чистый подход numpy.
Примечание. Некоторые из них (например, 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(5, 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
Я действительно хотел немного отличаться от принятого ответа. Я строил средство sklearn скользящего среднего для конвейера 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
Конечно, не нужно принимать постоянные значения для заполнения, но в большинстве случаев этого достаточно.
Талиб содержит простой инструмент скользящего среднего, а также другие подобные инструменты усреднения (то есть экспоненциальное скользящее среднее). Ниже сравнивается метод с некоторыми другими решениями.
%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)
Одно предостережение состоит в том, что реальное должно иметь элементы dtype = float. В противном случае возникает следующая ошибка
Исключение: настоящее не двойное