¿Cómo calcular la media móvil usando NumPy?

Parece que no hay una función que simplemente calcule el promedio móvil en numpy / scipy, lo que lleva a soluciones complicadas .

Mi pregunta es doble:

• ¿Cuál es la forma más fácil de implementar (correctamente) una media móvil con numpy?
• Dado que esto parece no trivial y propenso a errores, ¿hay alguna buena razón para no incluir las baterías en este caso?

Si lo que desea es un simple no ponderado promedio móvil, se puede implementar fácilmente con np.cumsum, que puede ser es métodos más rápido que la FFT en base:

EDITAR Se corrigió una indexación incorrecta de uno en uno detectada por Bean en el código. EDITAR

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

Así que supongo que la respuesta es: es realmente fácil de implementar, y tal vez numpy ya esté un poco hinchado con funcionalidad especializada.

La falta de NumPy de una función específica de dominio particular se debe quizás a la disciplina y fidelidad del Core Team a la directiva principal de NumPy: proporcionar un tipo de matriz N-dimensional , así como funciones para crear e indexar esas matrices. Como muchos objetivos fundamentales, este no es pequeño y NumPy lo hace de manera brillante.

El SciPy (mucho) más grande contiene una colección mucho mayor de bibliotecas específicas de dominio (llamadas subpaquetes por los desarrolladores de SciPy), por ejemplo, optimización numérica ( optimizar ), procesamiento de señales ( señal ) y cálculo integral ( integrar ).

Supongo que la función que busca está en al menos uno de los subpaquetes de SciPy ( quizás scipy.signal ); sin embargo, buscaría primero en la colección de scikits de ciencia ficción , identificaría los scikit relevantes y buscaría la función de interés allí.

Scikits son paquetes desarrollados de forma independiente basados ​​en NumPy / SciPy y dirigidos a una disciplina técnica en particular (por ejemplo, scikits-image , scikits-learn , etc.). Varios de ellos (en particular, el impresionante OpenOpt para optimización numérica) fueron muy apreciados, proyectos maduros mucho antes de elegir residir bajo la rúbrica de scikits relativamente nueva . La página de inicio de Scikits que le gustó arriba enumera alrededor de 30 de esos scikits , aunque al menos varios de ellos ya no están en desarrollo activo.

Seguir este consejo lo llevará a scikits-timeseries ; sin embargo, ese paquete ya no se encuentra en desarrollo activo; En efecto, Pandas se ha convertido, AFAIK, el de facto NumPy biblioteca de series de tiempo basada en .

Pandas tiene varias funciones que pueden usarse para calcular un promedio móvil ; el más simple de estos es probablemente rolling_mean , que usa así:

>>> # 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)

Ahora, simplemente llame a la función rolling_mean pasando el objeto Serie y un tamaño de ventana , que en mi ejemplo a continuación es de 10 días .

>>> 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

verificar que funcionó, por ejemplo, comparar los valores 10-15 en la serie original con la nueva Serie suavizada con media móvil

>>> 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

La función rolling_mean, junto con alrededor de una docena de otras funciones se agrupan informalmente en la documentación de Pandas bajo la rúbrica funciones de ventana móvil ; un segundo grupo de funciones relacionado en Pandas se conoce como funciones ponderadas exponencialmente (por ejemplo, ewma , que calcula el promedio ponderado móvil exponencial). El hecho de que este segundo grupo no esté incluido en el primero ( funciones de ventana móvil ) se debe quizás a que las transformaciones ponderadas exponencialmente no dependen de una ventana de longitud fija

Una forma sencilla de lograr esto es usando np.convolve. La idea detrás de esto es aprovechar la forma en que se calcula la convolución discreta y usarla para devolver una media móvil . Esto se puede hacer convolucionando con una secuencia denp.ones de una longitud igual a la longitud de la ventana deslizante que queremos.

Para ello podríamos definir la siguiente función:

def moving_average(x, w):
    return np.convolve(x, np.ones(w), 'valid') / w

Esta función tomará la convolución de la secuencia xy una secuencia de unos de longitud w. Tenga en cuenta que lo elegido modees validpara que el producto de convolución solo se proporcione para los puntos donde las secuencias se superponen por completo.

Algunos ejemplos:

x = np.array([5,3,8,10,2,1,5,1,0,2])

Para una media móvil con una ventana de longitud 2tendríamos:

moving_average(x, 2)
# array([4. , 5.5, 9. , 6. , 1.5, 3. , 3. , 0.5, 1. ])

Y para una ventana de largo 4:

moving_average(x, 4)
# array([6.5 , 5.75, 5.25, 4.5 , 2.25, 1.75, 2.  ])

¿Cómo convolvefunciona?

Echemos un vistazo más en profundidad a la forma en que se calcula la convolución discreta. La siguiente función tiene como objetivo replicar la formanp.convolve se calculan los valores de salida:

def mov_avg(x, w):
    for m in range(len(x)-(w-1)):
        yield sum(np.ones(w) * x[m:m+w]) / w 

Lo cual, para el mismo ejemplo anterior, también produciría:

list(mov_avg(x, 2))
# [4.0, 5.5, 9.0, 6.0, 1.5, 3.0, 3.0, 0.5, 1.0]

Entonces, lo que se está haciendo en cada paso es tomar el producto interno entre la matriz de unos y la ventana actual . En este caso la multiplicación por np.ones(w)es superflua dado que estamos tomando directamente lasum de la secuencia.

A continuación se muestra un ejemplo de cómo se calculan las primeras salidas para que sea un poco más claro. Supongamos que queremos una ventana de 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

Y la siguiente salida se calcularía como:

  [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

Y así sucesivamente, devolviendo un promedio móvil de la secuencia una vez que se hayan realizado todas las superposiciones.

Aquí hay una variedad de formas de hacer esto, junto con algunos puntos de referencia. Los mejores métodos son las versiones que utilizan código optimizado de otras bibliotecas. El bottleneck.move_meanmétodo probablemente sea el mejor en todos los aspectos. El scipy.convolveenfoque también es muy rápido, extensible y sintácticamente y conceptualmente simple, pero no se escala bien para valores de ventana muy grandes. El numpy.cumsummétodo es bueno si necesitas un puronumpy enfoque .

Nota: algunos de estos (p bottleneck.move_mean. Ej. ) No están centrados y cambiarán sus datos.

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 /programming/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)

Tiempo, ventana pequeña (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)

Tiempo, ventana grande (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)

Memoria, ventana pequeña (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

Memoria, ventana grande (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

Esta respuesta que usa Pandas está adaptada de arriba, ya que rolling_meanya no es parte de 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)

Ahora, simplemente llame a la función rollingen el marco de datos con un tamaño de ventana, que en mi ejemplo a continuación es de 10 días.

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

Siento que esto se puede resolver fácilmente usando el cuello de botella

Vea la muestra básica a continuación:

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)

Esto da media de movimiento a lo largo de cada eje.

• "mm" es la media móvil de "a".

• "ventana" es el número máximo de entradas a considerar para la media móvil.

• "min_count" es el número mínimo de entradas a considerar para la media móvil (por ejemplo, para el primer elemento o si la matriz tiene valores nan).

Lo bueno es que Bottleneck ayuda a lidiar con los valores nan y también es muy eficiente.

En caso de que desee cuidar las condiciones de los bordes con cuidado ( calcule la media solo a partir de los elementos disponibles en los bordes ), la siguiente función funcionará.

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

for i in range(len(Data)):
    Data[i, 1] = Data[i-lookback:i, 0].sum() / lookback

Prueba este fragmento de código. Creo que es más simple y funciona. retroceso es la ventana de la media móvil.

En el Data[i-lookback:i, 0].sum()he puesto 0para referirse a la primera columna del conjunto de datos pero puedes poner cualquier columna que quieras en caso de que tengas más de una columna.

De hecho, quería un comportamiento ligeramente diferente a la respuesta aceptada. Estaba construyendo un extractor de características de media móvil para una sklearntubería, por lo que necesitaba que la salida de la media móvil tuviera la misma dimensión que la entrada. Lo que quiero es que la media móvil asuma que la serie permanece constante, es decir, una media móvil de [1,2,3,4,5]con la ventana 2 daría[1.5,2.5,3.5,4.5,5.0] .

Para los vectores de columna (mi caso de uso) obtenemos

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

Y para matrices

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

Por supuesto, no es necesario asumir valores constantes para el relleno, pero hacerlo debería ser adecuado en la mayoría de los casos.

talib contiene una herramienta de promedio móvil simple, así como otras herramientas de promedio similares (es decir, promedio móvil exponencial). A continuación, se compara el método con algunas de las otras soluciones.

%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)

Una advertencia es que lo real debe tener elementos de dtype = float. De lo contrario, se genera el siguiente error

Excepción: real no es doble

Aquí hay una implementación rápida usando numba (tenga en cuenta los tipos). Tenga en cuenta que contiene nans donde se desplazó.

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())) 

media móvil

• invierta la matriz en i, y simplemente tome la media de ia n.

• use la comprensión de listas para generar mini matrices sobre la marcha.

x = np.random.randint(10, size=20)

def moving_average(arr, n):
    return [ (arr[:i+1][::-1][:n]).mean() for i, ele in enumerate(arr) ]
n = 5

moving_average(x, n)

Utilizo la solución de la respuesta aceptada , ligeramente modificada para tener la misma longitud para la salida que la entrada, o pandasla versión como se menciona en un comentario de otra respuesta. Resumo ambos aquí con un ejemplo reproducible para referencia futura:

import numpy as np
import pandas as pd

def moving_average(a, n):
    ret = np.cumsum(a, dtype=float)
    ret[n:] = ret[n:] - ret[:-n]
    return ret / n

def moving_average_centered(a, n):
    return pd.Series(a).rolling(window=n, center=True).mean().to_numpy()

A = [0, 0, 1, 2, 4, 5, 4]
print(moving_average(A, 3))    
# [0.         0.         0.33333333 1.         2.33333333 3.66666667 4.33333333]
print(moving_average_centered(A, 3))
# [nan        0.33333333 1.         2.33333333 3.66666667 4.33333333 nan       ]

Al comparar la solución a continuación con la que usa cumsum of numpy, esta toma casi la mitad del tiempo . Esto se debe a que no necesita pasar por toda la matriz para hacer el cumsum y luego hacer toda la resta. Además, el cumsum puede ser " peligroso " si la matriz es enorme y el número es enorme ( posible desbordamiento ). Por supuesto, también aquí existe el peligro, pero al menos se suman solo los números esenciales.

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