Media móvil o media móvil
¿Existe una función SciPy o una función o módulo NumPy para Python que calcule la media corriente de una matriz 1D dada una ventana específica?
NOTA: Las soluciones más eficientes pueden incluir
scipy.ndimage.uniform_filter1d(consulte esta respuesta ) o el uso de bibliotecas más nuevas, incluidas las de talibtalib.MA.
Usar np.convolve:
np.convolve(x, np.ones(N)/N, mode='valid')
Explicación
La media móvil es un caso de la operación matemática de convolución . Para la media móvil, desliza una ventana a lo largo de la entrada y calcula la media del contenido de la ventana. Para señales 1D discretas, la convolución es lo mismo, excepto que en lugar de la media se calcula una combinación lineal arbitraria, es decir, se multiplica cada elemento por un coeficiente correspondiente y se suman los resultados. Esos coeficientes, uno para cada posición en la ventana, a veces se denominan núcleo de convolución . La media aritmética de N valores es (x_1 + x_2 + ... + x_N) / N, por lo que el núcleo correspondiente es (1/N, 1/N, ..., 1/N), y eso es exactamente lo que obtenemos al usar np.ones(N)/N.
Bordes
El modeargumento de np.convolveespecifica cómo manejar los bordes. Elegí valideste modo porque creo que así es como la mayoría de la gente espera que funcione el modo de carrera, pero es posible que tengas otras prioridades. Aquí hay un gráfico que ilustra la diferencia entre los modos:
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()
Solución eficiente
La convolución es mucho mejor que el enfoque sencillo, pero (supongo) utiliza FFT y, por lo tanto, es bastante lento. Sin embargo, especialmente para calcular la media en ejecución, el siguiente enfoque funciona bien
def running_mean(x, N):
cumsum = numpy.cumsum(numpy.insert(x, 0, 0))
return (cumsum[N:] - cumsum[:-N]) / float(N)
El código para comprobar
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
Tenga en cuenta que numpy.allclose(result1, result2)dos Truemétodos son equivalentes. A mayor N, mayor diferencia en el tiempo.
Advertencia: aunque el cumsum es más rápido, habrá un mayor error de punto flotante que puede provocar que los resultados no sean válidos/incorrectos/inaceptables.
Los comentarios señalaron este problema de error de punto flotante aquí, pero lo hago más obvio aquí en la respuesta. .
# 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)
- Cuantos más puntos acumule, mayor será el error de punto flotante (por lo que 1e5 puntos es notable, 1e6 puntos es más significativo, más de 1e6 y es posible que desee restablecer los acumuladores)
- puedes hacer trampa usando
np.longdouble, pero tu error de punto flotante aún será significativo para una cantidad relativamente grande de puntos (alrededor de >1e5, pero depende de tus datos) - puedes trazar el error y verlo aumentar relativamente rápido
- la solución de convolución es más lenta pero no tiene esta pérdida de precisión de punto flotante
- la solución uniform_filter1d es más rápida que esta solución cumsum Y no tiene esta pérdida de precisión de punto flotante
Puedes usar 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:
- da la salida con la misma forma numerosa (es decir, número de puntos)
- permite múltiples formas de manejar el borde donde
'reflect'está el valor predeterminado, pero en mi caso, prefería'nearest'
También es bastante rápido (casi 50 veces más rápido np.convolvey entre 2 y 5 veces más rápido que el enfoque cumsum indicado anteriormente ):
%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
Aquí hay 3 funciones que le permiten comparar el error/velocidad de diferentes implementaciones:
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)]
Actualización: el siguiente ejemplo muestra la pandas.rolling_meanfunción anterior que se ha eliminado en versiones recientes de pandas. Un equivalente moderno de esa llamada a función usaría pandas.Series.rolling :
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 es más adecuado para esto que NumPy o SciPy. Su función Rolling_mean hace el trabajo cómodamente. También devuelve una matriz NumPy cuando la entrada es una matriz.
Es difícil superar rolling_meanel rendimiento con cualquier implementación personalizada de Python puro. A continuación se muestra un ejemplo de rendimiento frente a dos de las soluciones propuestas:
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
También hay buenas opciones sobre cómo lidiar con los valores de los bordes.
