Moyenne mobile ou moyenne glissante

Question

Existe-t-il une fonction SciPy, une fonction NumPy ou un module Python permettant de calculer la moyenne mobile d'un tableau 1D en fonction d'une fenêtre spécifique ?

Answer 1

UPDATE : des solutions plus efficaces ont été proposées, uniform_filter1d de scipy est probablement la meilleure parmi les bibliothèques tierces "standard", et certaines bibliothèques plus récentes ou spécialisées sont également disponibles.

---

Vous pouvez utiliser np.convolve pour ça :

np.convolve(x, np.ones(N)/N, mode='valid')

Explication

La moyenne mobile est un cas de l'opération mathématique de convolution . Pour la moyenne mobile, vous faites glisser une fenêtre le long de l'entrée et vous calculez la moyenne du contenu de la fenêtre. Pour les signaux 1D discrets, la convolution est la même chose, sauf qu'au lieu de la moyenne, vous calculez une combinaison linéaire arbitraire, c'est-à-dire que vous multipliez chaque élément par le coefficient correspondant et additionnez les résultats. Ces coefficients, un pour chaque position dans la fenêtre, sont parfois appelés la convolution. Noyau . La moyenne arithmétique de N valeurs est de (x_1 + x_2 + ... + x_N) / N donc le noyau correspondant est (1/N, 1/N, ..., 1/N) et c'est exactement ce que nous obtenons en utilisant np.ones(N)/N .

Bords

El mode argument de np.convolve spécifie comment traiter les bords. J'ai choisi le valid parce que je pense que c'est ainsi que la plupart des gens s'attendent à ce que le moyen de transport fonctionne, mais vous avez peut-être d'autres priorités. Voici un graphique qui illustre la différence entre les modes :

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

Answer 2

Une solution efficace

La convolution est bien meilleure que l'approche directe, mais (je suppose) elle utilise la FFT et est donc assez lente. Cependant, l'approche suivante fonctionne parfaitement pour le calcul de la moyenne mobile

def running_mean(x, N):
    cumsum = numpy.cumsum(numpy.insert(x, 0, 0)) 
    return (cumsum[N:] - cumsum[:-N]) / float(N)

Le code à vérifier

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

Notez que numpy.allclose(result1, result2) es True les deux méthodes sont équivalentes. Plus N est grand, plus la différence de temps est importante.

avertissement : bien que le cumsum soit plus rapide, il y aura une augmentation des erreurs de virgule flottante qui peuvent rendre vos résultats invalides/incorrects/inacceptables.

Les commentaires ont mis en évidence ce problème d'erreur de virgule flottante ici, mais je le rends plus évident ici dans la réponse. .

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

• plus vous accumulez de points, plus l'erreur en virgule flottante est importante (ainsi 1e5 points est notable, 1e6 points est plus significatif, plus de 1e6 et vous voudrez peut-être remettre les accumulateurs à zéro)
• vous pouvez tricher en utilisant np.longdouble mais votre erreur en virgule flottante sera toujours significative pour un nombre relativement important de points (environ >1e5 mais cela dépend de vos données).
• vous pouvez tracer l'erreur et voir qu'elle augmente relativement vite
• la solution convolve est plus lent mais n'a pas cette perte de précision en virgule flottante
• la solution uniform_filter1d est plus rapide que cette solution cumsum ET n'a pas cette perte de précision en virgule flottante

Answer 3

Mise à jour : L'exemple ci-dessous montre l'ancien pandas.rolling_mean qui a été supprimée dans les versions récentes de pandas. Un équivalent moderne de cet appel de fonction utiliserait 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 est plus adapté à cette tâche que NumPy ou SciPy. Sa fonction moyenne mobile fait le travail de manière pratique. Il retourne également un tableau NumPy lorsque l'entrée est un tableau.

Il est difficile de battre rolling_mean en termes de performances par rapport à toute implémentation personnalisée en Python pur. Voici un exemple de performance par rapport à deux des solutions proposées :

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

Il existe également des options intéressantes quant à la manière de traiter les valeurs de bord.

Answer 4

Vous pouvez utiliser scipy.ndimage.filters.uniform_filter1d :

import numpy as np
from scipy.ndimage.filters import uniform_filter1d
N = 1000
x = np.random.random(100000)
y = uniform_filter1d(x, size=N)

uniform_filter1d :

• donne la sortie avec la même forme numpy (i.e. le nombre de points)
• permet de gérer de multiples façons la frontière où 'reflect' est la valeur par défaut, mais dans mon cas, j'ai plutôt voulu que 'nearest'

Il est également assez rapide (près de 50 fois plus rapide que les np.convolve et de 2 à 5 fois plus rapide que l'approche cumsum donnée ci-dessus ) :

%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

Voici 3 fonctions qui vous permettent de comparer l'erreur/la vitesse de différentes implémentations :

from __future__ import division
import numpy as np
import scipy.ndimage.filters as ndif
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 ndif.uniform_filter1d(x, N, mode='constant', origin=-(N//2))[:-(N-1)]

Answer 5

Vous pouvez calculer une moyenne mobile avec :

import numpy as np

def runningMean(x, N):
    y = np.zeros((len(x),))
    for ctr in range(len(x)):
         y[ctr] = np.sum(x[ctr:(ctr+N)])
    return y/N

Mais c'est lent.

Heureusement, numpy inclut un convolve que nous pouvons utiliser pour accélérer les choses. La moyenne mobile est équivalente à la convolution x avec un vecteur qui est N de long, dont tous les membres sont égaux à 1/N . L'implémentation numpy de convolve inclut le transitoire de départ, vous devez donc supprimer les N-1 premiers points :

def runningMeanFast(x, N):
    return np.convolve(x, np.ones((N,))/N)[(N-1):]

Sur ma machine, la version rapide est 20 à 30 fois plus rapide, en fonction de la longueur du vecteur d'entrée et de la taille de la fenêtre de calcul de la moyenne.

Notez que convolve comprend un 'same' qui semble devoir résoudre le problème du transitoire de départ, mais qui le divise entre le début et la fin.