Comment calculer la moyenne mobile en utilisant NumPy?

Il ne semble y avoir aucune fonction qui calcule simplement la moyenne mobile sur numpy / scipy, conduisant à des solutions alambiquées .

Ma question est double:

• Quel est le moyen le plus simple d'implémenter (correctement) une moyenne mobile avec numpy?
• Étant donné que cela semble non trivial et sujet aux erreurs, y a-t-il une bonne raison de ne pas inclure les piles dans ce cas?

Si vous voulez juste une moyenne, vous pouvez facilement non pondérée simple mouvement de mettre en œuvre avec np.cumsum, ce qui peut être est des méthodes basées plus rapide que FFT:

EDIT Correction d'une mauvaise indexation ponctuelle repérée par Bean dans le code. ÉDITER

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

Donc, je suppose que la réponse est: c'est vraiment facile à mettre en œuvre, et peut-être que numpy est déjà un peu gonflé de fonctionnalités spécialisées.

L'absence de fonction spécifique à un domaine par NumPy est peut-être due à la discipline de l'équipe principale et à la fidélité à la directive principale de NumPy: fournir un type de tableau à N dimensions , ainsi que des fonctions pour créer et indexer ces tableaux. Comme de nombreux objectifs fondamentaux, celui-ci n'est pas petit et NumPy le fait avec brio.

Le SciPy (beaucoup) plus grand contient une collection beaucoup plus grande de bibliothèques spécifiques à un domaine (appelées sous- packages par les développeurs de SciPy) - par exemple, l'optimisation numérique ( optimiser ), le traitement du signal ( signal ) et le calcul intégral ( intégrer ).

Je suppose que la fonction que vous recherchez se trouve dans au moins l'un des sous-packages SciPy ( scipy.signal peut-être); Cependant, je chercherais d'abord dans la collection de scikits SciPy , identifierais le (s) scikit (s) pertinent (s) et chercherais la fonction qui l'intéresse.

Les scikits sont des packages développés indépendamment basés sur NumPy / SciPy et dirigés vers une discipline technique particulière (par exemple, scikits-image , scikits-learn , etc.) Plusieurs d'entre eux étaient (en particulier, l'impressionnant OpenOpt pour l'optimisation numérique) étaient très appréciés, projets mûrs bien avant de choisir de résider sous la rubrique relativement nouvelle des scikits . La page d'accueil de Scikits comme ci-dessus répertorie environ 30 de ces scikits , bien qu'au moins plusieurs d'entre eux ne soient plus en développement actif.

Suivre ce conseil vous conduirait à des scikits-timeseries ; cependant, ce paquet n'est plus en cours de développement actif; En effet, Pandas est devenu, AFAIK, la bibliothèque de séries chronologiques de facto basée sur NumPy .

Pandas a plusieurs fonctions qui peuvent être utilisées pour calculer une moyenne mobile ; le plus simple d'entre eux est probablement rolling_mean , que vous utilisez comme ceci:

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

Maintenant, appelez simplement la fonction rolling_mean en passant l'objet Series et une taille de fenêtre , qui dans mon exemple ci-dessous est de 10 jours .

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

vérifier que cela a fonctionné - par exemple, comparez les valeurs de 10 à 15 dans la série d'origine par rapport à la nouvelle série lissée avec une moyenne de roulement

>>> 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 fonction rolling_mean, ainsi qu'une douzaine d'autres fonctions sont regroupées de manière informelle dans la documentation Pandas sous la rubrique Fonctions de la fenêtre mobile ; un deuxième groupe de fonctions apparentées dans Pandas est appelé fonctions à pondération exponentielle (par exemple, ewma , qui calcule une moyenne pondérée en mouvement exponentiel). Le fait que ce deuxième groupe ne soit pas inclus dans le premier ( fonctions de fenêtre mobile ) est peut-être dû au fait que les transformations à pondération exponentielle ne reposent pas sur une fenêtre de longueur fixe

Un moyen simple d'y parvenir est d'utiliser np.convolve. L'idée derrière cela est d'exploiter la façon dont la convolution discrète est calculée et de l'utiliser pour renvoyer une moyenne mobile . Cela peut être fait en convoluant avec une séquence denp.ones une longueur égale à la longueur de la fenêtre glissante que nous voulons.

Pour ce faire, nous pourrions définir la fonction suivante:

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

Cette fonction prendra la convolution de la séquence xet une séquence de uns de longueur w. Notez que le choix modeest validtel que le produit de convolution n'est donné que pour les points où les séquences se chevauchent complètement.

Quelques exemples:

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

Pour une moyenne mobile avec une fenêtre de longueur, 2nous aurions:

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

Et pour une fenêtre de longueur 4:

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

Comment ça convolvemarche?

Regardons plus en profondeur la façon dont la convolution discrète est calculée. La fonction suivante vise à reproduire la méthode de np.convolvecalcul des valeurs de sortie:

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

Ce qui, pour le même exemple ci-dessus, donnerait également:

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

Donc ce qui est fait à chaque étape est de prendre le produit interne entre le tableau de uns et la fenêtre courante . Dans ce cas, la multiplication par np.ones(w)est superflue étant donné que nous prenons directement lesum de la séquence.

Ci-dessous est un exemple de la façon dont les premières sorties sont calculées pour que ce soit un peu plus clair. Supposons que nous voulons une fenêtre 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

Et la sortie suivante serait calculée comme:

  [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

Et ainsi de suite, renvoyer une moyenne mobile de la séquence une fois que tous les chevauchements ont été effectués.

Voici une variété de façons de le faire, ainsi que quelques repères. Les meilleures méthodes sont les versions utilisant du code optimisé d'autres bibliothèques. La bottleneck.move_meanméthode est probablement la meilleure partout. L' scipy.convolveapproche est également très rapide, extensible et syntaxiquement et conceptuellement simple, mais ne s'adapte pas bien aux très grandes valeurs de fenêtre. La numpy.cumsumméthode est bonne si vous avez besoin d'une numpyapproche pure .

Remarque: certains d'entre eux (par exemple bottleneck.move_mean) ne sont pas centrés et déplaceront vos données.

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)

Calendrier, petite fenêtre (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)

Calendrier, grande fenêtre (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)

Mémoire, petite fenêtre (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

Mémoire, grande fenêtre (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

Cette réponse utilisant Pandas est adaptée d'en haut, car elle rolling_meanne fait plus partie 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)

Maintenant, appelez simplement la fonction rollingsur le dataframe avec une taille de fenêtre, qui dans mon exemple ci-dessous est de 10 jours.

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

Je pense que cela peut être facilement résolu en utilisant un goulot d'étranglement

Voir l'exemple de base ci-dessous:

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)

Cela donne une moyenne de déplacement le long de chaque axe.

• "mm" est la moyenne mobile de "a".

• "window" est le nombre maximum d'entrées à considérer pour la moyenne mobile.

• "min_count" est le nombre minimal d'entrées à considérer pour la moyenne mobile (par exemple pour le premier élément ou si le tableau a des valeurs nan).

La bonne partie est que Bottleneck aide à gérer les valeurs nanométriques et il est également très efficace.

Dans le cas où vous voulez prendre soin des conditions de bord avec soin ( calculer la moyenne uniquement à partir des éléments disponibles sur les bords ), la fonction suivante fera l'affaire.

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

Essayez ce morceau de code. Je pense que c'est plus simple et fait le travail. lookback est la fenêtre de la moyenne mobile.

Dans le Data[i-lookback:i, 0].sum()J'ai mis 0pour faire référence à la première colonne de l'ensemble de données, mais vous pouvez mettre n'importe quelle colonne que vous aimez au cas où vous auriez plus d'une colonne.

Je voulais en fait un comportement légèrement différent de la réponse acceptée. Je construisais un extracteur de caractéristiques de moyenne mobile pour un sklearnpipeline, j'ai donc exigé que la sortie de la moyenne mobile ait la même dimension que l'entrée. Ce que je veux, c'est que la moyenne mobile suppose que la série reste constante, c'est-à-dire qu'une moyenne mobile de [1,2,3,4,5]avec la fenêtre 2 donnerait[1.5,2.5,3.5,4.5,5.0] .

Pour les vecteurs de colonne (mon cas d'utilisation), nous obtenons

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

Et pour les tableaux

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

Bien sûr, il n'est pas nécessaire de supposer des valeurs constantes pour le remplissage, mais cela devrait être adéquat dans la plupart des cas.

talib contient un outil de moyenne mobile simple, ainsi que d'autres outils de calcul de moyenne similaires (c.-à-d. moyenne mobile exponentielle). Ci-dessous compare la méthode à certaines des autres solutions.

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

Une mise en garde est que le réel doit avoir des éléments de dtype = float. Sinon, l'erreur suivante est générée

Exception: le réel n'est pas double

Voici une implémentation rapide utilisant numba (attention aux types). Notez qu'il contient des nans lorsqu'ils sont décalés.

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

moyenne mobile

• inverser le tableau en i, et prendre simplement la moyenne de i à n.

• Utilisez la compréhension de liste pour générer des mini-tableaux à la volée.

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)

J'utilise soit la solution de la réponse acceptée , légèrement modifiée pour avoir la même longueur pour la sortie que l'entrée, soit pandasla version comme mentionné dans un commentaire d'une autre réponse. Je résume les deux ici avec un exemple reproductible pour référence future:

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       ]

En comparant la solution ci-dessous avec celle qui utilise du sperme de numpy, celle-ci prend presque la moitié du temps . En effet, il n'est pas nécessaire de parcourir tout le tableau pour faire le sperme, puis de faire toute la soustraction. De plus, le sperme peut être " dangereux " si le tableau est énorme et le nombre est énorme ( débordement possible ). Bien sûr, ici aussi le danger existe, mais au moins ne sont additionnés que les nombres essentiels.

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