Wie berechnet man den gleitenden Durchschnitt mit NumPy?

Wenn Sie nur einen einfachen, nicht gewichteten gleitenden Durchschnitt wünschen, können Sie ihn einfach mit np.cumsum implementieren könnte sein ist schneller als FFT-basierte Methoden:

EDITEine falsche Codierung, die von Bean im Code entdeckt wurde, wurde korrigiert.EDIT

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

Ich denke, die Antwort ist: Es ist wirklich einfach zu implementieren, und vielleicht ist Numpy schon etwas aufgebläht mit spezieller Funktionalität.

Das Fehlen einer bestimmten domänenspezifischen Funktion von NumPy beruht möglicherweise auf der Disziplin des Kernteams und der Vertrauenswürdigkeit von NumPys Hauptanweisung: N-dimensionalen Array-Typ angeben sowie Funktionen zum Erstellen und Indizieren dieser Arrays. Wie viele grundlegende Ziele ist auch dieses Ziel nicht klein und NumPy macht es hervorragend.

Das (viel) größere SciPy enthält eine viel größere Sammlung domänenspezifischer Bibliotheken (genannt subpackages von SciPy devs) - zum Beispiel numerische Optimierung (Optimieren), Signalverarbeitung (Signal) und Integralrechnung (Integrieren).

Meine Vermutung ist, dass die Funktion, nach der Sie suchen, in mindestens einem der SciPy-Subpakete enthalten ist (scipy.signal vielleicht); Ich würde jedoch zuerst in der Sammlung von SciPy-Scikits suchen, die relevanten Scikits identifizieren und dort nach der Funktion suchen, die von Interesse ist.

Scikits sind unabhängig entwickelte Pakete, die auf NumPy/SciPy basieren und auf eine bestimmte technische Disziplin gerichtet sind (z. B. scikits-image, scikits-learn usw. Mehrere davon (vor allem die großartigen OpenOpt für numerische Optimierung) waren hoch angesehene, reife Projekte, lange bevor sie sich für die relativ neue Rubrik scikits entschieden haben. Auf der Homepage von Scikits wurden ungefähr 30 solcher scikits -Angaben aufgeführt, obwohl mindestens einige von ihnen nicht mehr aktiv entwickelt werden. 

Wenn Sie diesen Rat befolgen, würden Sie zu scikits-timeseries; Dieses Paket befindet sich jedoch nicht mehr in der aktiven Entwicklung. Tatsächlich wurde Pandas aus AFAIK die de factoNumPy - basierte Zeitreihenbibliothek.

Pandas hat mehrere Funktionen, mit denen ein gleitender Durchschnitt berechnet werden kann; Die einfachste davon ist wahrscheinlich rolling_mean, die Sie so verwenden:

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

Rufen Sie nun einfach die Funktion rolling_mean auf, die im Series-Objekt übergeben wird, und ein Fenstergröße, das in meinem folgenden Beispiel 10 Tage ist.

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

vergewissern Sie sich, dass es funktioniert - vergleichen Sie z. B. die Werte 10 - 15 in der Originalserie mit der neuen Serie, die mit dem rollenden Mittelwert geglättet wurde

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

Die Funktion rolling_mean sowie etwa ein Dutzend anderer Funktionen sind in der Pandas-Dokumentation informell unter der Rubrik Bewegungsfenster Funktionen zusammengefasst. Eine zweite, verwandte Gruppe von Funktionen in Pandas wird als exponentiell gewichtete Funktionen bezeichnet (z. B. ewma, die den exponentiell bewegten gewichteten Durchschnitt berechnet. Die Tatsache, dass diese zweite Gruppe nicht in der ersten (sich bewegendes Fenster Funktionen) enthalten ist, ist möglicherweise darauf zurückzuführen, dass die exponentiell gewichteten Transformationen nicht auf einem Fenster fester Länge basieren

Eine einfache Möglichkeit, dies zu erreichen, ist die Verwendung von np.convolve . Die Idee dahinter ist, die Art und Weise, wie diskrete Faltung berechnet wird, zu nutzen und a zurückzugeben rollmittel. Dies kann durch Falten mit einer Sequenz von np.ones mit einer Länge erreicht werden, die der gewünschten Länge des Schiebefensters entspricht.

Dazu können wir folgende Funktion definieren:

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

Diese Funktion nimmt die Faltung der Sequenz x und einer Sequenz von Einsen der Länge w. Beachten Sie, dass das gewählte modevalid ist, sodass das Faltungsprodukt nur für Punkte angegeben wird, bei denen sich die Sequenzen vollständig überlappen.

Anwendungsfall

Einige Beispiele:

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

Für einen gleitenden Durchschnitt mit einem Fenster der Länge 2 Hätten wir:

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

Und für ein Fenster der Länge 4:

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

Einzelheiten

Betrachten wir die Art und Weise, wie die diskrete Faltung berechnet wird. Die folgende Funktion soll die Art und Weise replizieren, wie np.convolve Die Ausgabewerte berechnet:

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

Was für dasselbe Beispiel auch ergeben würde:

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

Was also bei jedem Schritt getan wird, ist, das innere Produkt zwischen der Reihe von Einsen und dem Strom zu nehmen fenster. In diesem Fall erübrigt sich die Multiplikation mit np.ones(w), da wir direkt das sum der Sequenz nehmen.

Das Folgende ist ein Beispiel dafür, wie die ersten Ausgaben berechnet werden, damit sie ein wenig klarer werden. Nehmen wir an, wir wollen ein Fenster mit 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

Und die folgende Ausgabe würde wie folgt berechnet:

  [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

Und so weiter, indem ein gleitender Durchschnitt der Sequenz zurückgegeben wird, sobald alle Überlappungen durchgeführt wurden.

Diese Antwort mit Pandas wird von oben angepasst, da rolling_mean nicht mehr Bestandteil von Pandas ist

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

Rufen Sie nun einfach die Funktion rolling auf dem Datenrahmen mit einer Fenstergröße auf, die in meinem folgenden Beispiel 10 Tage beträgt.

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

Ich glaube, das lässt sich leicht mit Engpass lösen.

Siehe grundlegendes Beispiel unten:

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)

Dies gibt den Mittelwert der Bewegung entlang jeder Achse.

• "mm" ist der gleitende Mittelwert für "a". 

• "Fenster" ist die maximale Anzahl von Einträgen, die für den gleitenden Mittelwert berücksichtigt werden müssen. 

• "min_count" ist die minimale Anzahl von Einträgen, die beim Verschieben des Mittelwerts berücksichtigt werden müssen (z. B. für das erste Element oder wenn das Array nan-Werte hat).

Das Gute daran ist, dass Bottleneck hilft, mit Nan-Werten umzugehen, und es ist auch sehr effizient.

Wenn Sie die Randbedingungen sorgfältig berücksichtigen möchten ( Mittelwert nur aus verfügbaren Elementen an Kanten berechnen), wird die folgende Funktion den Trick ausführen. 

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

Eigentlich wollte ich ein etwas anderes Verhalten als die akzeptierte Antwort. Ich habe einen Feature-Extraktor für den gleitenden Durchschnitt für eine sklearn -Pipeline erstellt. Daher musste die Ausgabe des gleitenden Durchschnitts dieselbe Dimension wie die Eingabe haben. Ich möchte, dass der gleitende Durchschnitt annimmt, dass die Reihe konstant bleibt, dh ein gleitender Durchschnitt von [1,2,3,4,5] mit Fenster 2 würde [1.5,2.5,3.5,4.5,5.0] ergeben.

Für Spaltenvektoren (mein Anwendungsfall) erhalten wir

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

Und für Arrays

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

Natürlich muss man keine konstanten Werte für die Polsterung annehmen, aber dies sollte in den meisten Fällen ausreichend sein.