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

Parece que no hay una función que simplemente calcula el promedio móvil en números y datos, lo que lleva a soluciones complicadas .

Mi pregunta es doble:

  • ¿Cuál es la forma más fácil de (correctamente) implementar un promedio móvil con números?
  • 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 solo desea un promedio móvil no ponderado directo, puede implementarlo fácilmente con np.cumsum , que puede ser más rápido que los métodos basados ​​en FFT:

EDITAR Se corrigió una indexación incorrecta off-by-one 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 muy fácil de implementar, y tal vez ya esté un poco hinchado con funcionalidad especializada.

La falta de una función específica de dominio en particular de NumPy 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 objectives fundamentales, este no es pequeño, y NumPy lo hace de manera shiny.

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

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

Los Scikits son paquetes desarrollados independientemente basados ​​en NumPy / SciPy y dirigidos a una disciplina técnica particular (por ejemplo, scikits-image , scikits-learn , etc.) Varios de estos fueron (en particular, el asombroso OpenOpt para optimización numérica) fueron altamente considerados, proyectos maduros mucho antes de elegir residir bajo la rúbrica relativamente nueva de scikits . A la página de inicio de Scikits le gustaron las listas anteriores de unos 30 de esos scikits , aunque al menos varios de ellos ya no están en desarrollo activo.

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

Las pandas tienen varias funciones que pueden usarse para calcular un promedio móvil ; El más simple de estos es probablemente rolling_mean , que usas 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 el tamaño de una ventana , que en mi ejemplo a continuación son 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 

verifique que funcionó, p. ej., valores comparados 10 – 15 en la serie original frente a 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 una docena más o menos de otra función, se agrupan de manera informal en la documentación de Pandas debajo de las funciones de la ventana de movimiento de la rúbrica; un segundo grupo relacionado de funciones en Pandas se conoce como funciones ponderadas exponencialmente (por ejemplo, ewma , que calcula el promedio ponderado exponencialmente en movimiento). El hecho de que este segundo grupo no esté incluido en el primero (funciones de la ventana móvil ) es quizás porque las transformaciones ponderadas exponencialmente no se basan en una ventana de longitud fija

Una forma sencilla de lograr esto es mediante 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 por medio de una secuencia de np.ones de una longitud igual a la longitud de la ventana deslizante que deseamos.

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 x y una secuencia de unidades de longitud w . Tenga en cuenta que el mode elegido es valid por lo que el producto de convolución solo se otorga para los puntos donde las secuencias se superponen completamente.


Caso de uso

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 2 tendríamos:

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

Y para una ventana de longitud 4 :

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

Detalles

Veamos más a fondo la forma en que se calcula la convolución discreta. La siguiente función apunta a replicar la manera en que np.convolve está calculando 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 

Que, para el mismo ejemplo anterior también daría lugar 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 la sum 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 computarí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 han realizado todas las superposiciones.

Esta respuesta utilizando Pandas se adapta desde arriba, ya que rolling_mean ya 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 que está rolling en el dataframe con un tamaño de ventana, que en mi ejemplo a continuación son 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 

En caso de que quiera cuidar las condiciones de los bordes cuidadosamente ( calcular la media solo de los elementos disponibles en los bordes ), la siguiente función hará el truco.

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

Siento que esto se puede resolver fácilmente usando un 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 una 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 que se deben tener en cuenta 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 de nan y también es muy eficiente.