¿Cómo calcular el promedio móvil/móvil usando python + NumPy/SciPy?
Parece que no hay ninguna función que simplemente calcule la media móvil en numpy/scipy, lo que lleva a soluciones complicadas .
Mi pregunta es doble:
- ¿Cuál es la forma más sencilla de implementar (correctamente) una media móvil con numpy?
- Dado que esto no parece trivial y propenso a errores, ¿existe una buena razón para no incluir las baterías en este caso?
Si solo desea un promedio móvil no ponderado sencillo, 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 uno por 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])
Entonces supongo que la respuesta es: es realmente fácil de implementar, y tal vez numpy ya esté un poco repleto de funcionalidades especializadas.
Una forma sencilla de lograrlo es mediante el uso de np.convolve. La idea detrás de esto es aprovechar la forma en que se calcula la convolución discreta y utilizarla para devolver una media móvil . Esto se puede hacer convolucionando con una secuencia de np.onesuna 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 la opción elegida modees validque 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 longitud 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 profundo a la forma en que se calcula la convolución discreta. La siguiente función tiene como objetivo replicar la forma en que np.convolvese 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 hace 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 ya que estamos tomando directamente el sumde la secuencia.
A continuación se muestra un ejemplo de cómo se calculan las primeras salidas para que quede 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 el siguiente resultado 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 una media móvil de la secuencia una vez que se hayan realizado todas las superposiciones.
La falta de NumPy de una función específica de dominio particular se debe quizás a la disciplina y fidelidad del equipo central 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 (mucho) más grande SciPy 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 ( tal vez scipy.signal ); 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 de forma independiente basados en NumPy/SciPy y dirigidos a una disciplina técnica 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 relativamente nueva rúbrica scikits . La página de inicio de Scikits que aparece arriba enumera alrededor de 30 scikits de este tipo , aunque al menos varios de ellos ya no están en desarrollo activo.
Seguir este consejo le llevará a scikits-timeseries ; sin embargo, ese paquete ya no se encuentra en desarrollo activo; De hecho, Pandas se ha convertido, AFAIK, en la biblioteca de series temporales de facto basada en NumPy .
Pandas tiene varias funciones que se pueden utilizar para calcular una media móvil ; El más simple de ellos es probablemente Rolling_mean , que se 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 versus 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 aproximadamente 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 relacionado de funciones 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.
