webentwicklung-frage-antwort-db.com.de

Wie kann ich numpy.correlate für die Autokorrelation verwenden?

Ich muss eine automatische Korrelation einer Menge von Zahlen durchführen, die, wie ich es verstehe, nur die Korrelation der Menge mit sich selbst ist. 

Ich habe es mit der Korrelationsfunktion von numpy ausprobiert, aber ich glaube nicht an das Ergebnis, da es fast immer einen Vektor gibt, bei dem die erste Zahl nicht die größte Zahl ist, wie es sein sollte.

Diese Frage besteht also aus zwei Fragen:

  1. Was genau macht numpy.correlate?
  2. Wie kann ich es (oder etwas anderes) für die automatische Korrelation verwenden?
81
Ben

Um Ihre erste Frage zu beantworten, führt numpy.correlate(a, v, mode) die Faltung von a mit der Umkehrung von v durch und gibt die Ergebnisse durch den angegebenen Modus aus. Die Definition von Faltung , C (t) = ∑ -∞ <i <∞  einichvt + i wo -∞ <t <∞, ermöglicht Ergebnisse von -∞ bis allows, aber offensichtlich können Sie kein unendlich langes Array speichern. Es muss also abgeschnitten werden, und hier kommt der Modus ins Spiel. Es gibt drei verschiedene Modi: Full, Same und Valid: 

  • "full" -Modus gibt Ergebnisse für jede t zurück, wobei sich a und v etwas überlappen. 
  • Der "gleiche" Modus gibt ein Ergebnis mit derselben Länge wie der kürzeste Vektor (a oder v) zurück. 
  • Der Modus "gültig" gibt nur dann Ergebnisse zurück, wenn sich a und v vollständig überlappen. Die Dokumentation für numpy.convolve enthält weitere Informationen zu den Modi.

Für Ihre zweite Frage denke ich, dass numpy.correlateis Ihnen die Autokorrelation gibt. Die Autokorrelation wird verwendet, um herauszufinden, wie ähnlich sich ein Signal oder eine Funktion bei einer bestimmten Zeitdifferenz ist. Bei einer Zeitdifferenz von 0 sollte die Autokorrelation am höchsten sein, da das Signal mit sich selbst identisch ist. Daher haben Sie erwartet, dass das erste Element im Ergebniskorridor der Autokorrelation das größte Element ist. Die Korrelation beginnt jedoch nicht bei einer Zeitdifferenz von 0. Sie beginnt bei einer negativen Zeitdifferenz, schließt auf 0 und wird dann positiv. Das heißt, Sie haben erwartet: 

autokorrelation (a) = ∑ -∞ <i <∞  einichvt + i wobei 0 <= t <∞ 

Aber was du hast, war: 

autokorrelation (a) = ∑ -∞ <i <∞  einichvt + i wo -∞ <t <∞

Was Sie tun müssen, ist die letzte Hälfte Ihres Korrelationsergebnisses zu nehmen, und dies sollte die Autokorrelation sein, nach der Sie suchen. Eine einfache Python-Funktion dazu wäre:

def autocorr(x):
    result = numpy.correlate(x, x, mode='full')
    return result[result.size/2:]

Sie müssen natürlich eine Fehlerprüfung durchführen, um sicherzustellen, dass x tatsächlich ein 1-d-Array ist. Auch diese Erklärung ist wahrscheinlich nicht die mathematisch strengste. Ich habe Unendlichkeiten herumgeworfen, weil die Definition der Faltung sie verwendet, aber das gilt nicht unbedingt für die Autokorrelation. Der theoretische Teil dieser Erklärung mag etwas ungewöhnlich sein, aber hoffentlich sind die praktischen Ergebnisse hilfreich. DieseSeiten auf Autokorrelation sind sehr hilfreich und können Ihnen einen viel besseren theoretischen Hintergrund bieten, wenn Sie nichts dagegen haben, durch die Notation und die schweren Konzepte zu blättern.

93
A. Levy

Auto-Korrelation gibt es in zwei Versionen: statistisch und Faltung. Sie tun beide dasselbe, bis auf ein kleines Detail: Die statistische Version ist normalisiert und liegt im Intervall [-1,1]. Hier ist ein Beispiel, wie Sie das statistische Verfahren durchführen:

def acf(x, length=20):
    return numpy.array([1]+[numpy.corrcoef(x[:-i], x[i:])[0,1]  \
        for i in range(1, length)])
15
jonathf

Verwenden Sie die Funktion numpy.corrcoef anstelle von numpy.correlate , um die statistische Korrelation für eine Verzögerung von t zu berechnen:

def autocorr(x, t=1):
    return numpy.corrcoef(numpy.array([x[:-t], x[t:]]))

Da ich gerade auf dasselbe Problem gestoßen bin, möchte ich gerne ein paar Zeilen Code mit Ihnen teilen. In der Tat gibt es mehrere ähnliche Beiträge zur Autokorrelation in stackoverflow. Wenn Sie die Autokorrelation als a(x, L) = sum(k=0,N-L-1)((xk-xbar)*(x(k+L)-xbar))/sum(k=0,N-1)((xk-xbar)**2) definieren [dies ist die Definition in der a_correlate-Funktion von IDL und stimmt mit dem überein, was ich in Antwort 2 von question # 12269834 ] sehe, dann scheint das Folgende die richtigen Ergebnisse zu geben:

import numpy as np
import matplotlib.pyplot as plt

# generate some data
x = np.arange(0.,6.12,0.01)
y = np.sin(x)
# y = np.random.uniform(size=300)
yunbiased = y-np.mean(y)
ynorm = np.sum(yunbiased**2)
acor = np.correlate(yunbiased, yunbiased, "same")/ynorm
# use only second half
acor = acor[len(acor)/2:]

plt.plot(acor)
plt.show()

Wie Sie sehen, habe ich dies mit einer Sinuskurve und einer gleichmäßigen Zufallsverteilung getestet, und beide Ergebnisse sehen so aus, als würde ich sie erwarten. Beachten Sie, dass ich mode="same" anstelle von mode="full" wie die anderen verwendet habe. 

10
maschu

Ihre Frage 1 wurde bereits in mehreren hervorragenden Antworten ausführlich diskutiert.

Ich dachte, ich möchte Ihnen einige Codezeilen mitteilen, mit denen Sie die Autokorrelation eines Signals nur auf der Grundlage der mathematischen Eigenschaften der Autokorrelation berechnen können. Das heißt, die Autokorrelation kann auf folgende Weise berechnet werden:

  1. subtrahieren Sie den Mittelwert vom Signal und erhalten Sie ein unverzerrtes Signal

  2. berechnen Sie die Fourier-Transformation des unverzerrten Signals

  3. berechnen Sie die spektrale Leistungsdichte des Signals, indem Sie die quadratische Norm jedes Werts der Fouriertransformation des unverzerrten Signals verwenden

  4. berechnen Sie die inverse Fourier-Transformation der spektralen Leistungsdichte

  5. normalisieren Sie die inverse Fourier-Transformation der spektralen Leistungsdichte durch die Summe der Quadrate des unverzerrten Signals und nehmen Sie nur die Hälfte des resultierenden Vektors

Der Code dafür ist folgender:

def autocorrelation (x) :
    """
    Compute the autocorrelation of the signal, based on the properties of the
    power spectral density of the signal.
    """
    xp = x-np.mean(x)
    f = np.fft.fft(xp)
    p = np.array([np.real(v)**2+np.imag(v)**2 for v in f])
    pi = np.fft.ifft(p)
    return np.real(pi)[:x.size/2]/np.sum(xp**2)
8
Ruggero

Ich denke, es gibt zwei Dinge, die dieses Thema verwirren:

  1. statistische vs. Signalverarbeitungsdefinition: Wie andere darauf hingewiesen haben, normalisieren wir in der Statistik die Autokorrelation in [-1,1]. 
  2. teilweise v.s. Nicht-partieller Mittelwert/Abweichung: Verschiebt sich die Zeitreihe um eine Verzögerung> 0, wird ihre Überlappungsgröße immer <ursprüngliche Länge sein. Verwenden wir den Mittelwert und den Standardwert des Originals (nicht partiell), oder berechnen Sie immer einen neuen Mittelwert und Standardwert, indem Sie die sich ständig ändernde Überlappung (teilweise) verwenden. (Es gibt wahrscheinlich eine formale Bezeichnung dafür, aber ich werde jetzt "partiell" verwenden).

Ich habe 5 Funktionen erstellt, die die automatische Korrelation eines 1d-Arrays mit partiellen Vs berechnen. nicht partielle Unterscheidungen. Einige verwenden Formeln aus Statistiken, andere verwenden Korrelate im Sinne der Signalverarbeitung, was auch über FFT erfolgen kann. Bei allen Ergebnissen handelt es sich jedoch um automatische Korrelationen in der Definition von statistics, sodass sie veranschaulichen, wie sie miteinander verknüpft sind. Code unten:

import numpy
import matplotlib.pyplot as plt

def autocorr1(x,lags):
    '''numpy.corrcoef, partial'''

    corr=[1. if l==0 else numpy.corrcoef(x[l:],x[:-l])[0][1] for l in lags]
    return numpy.array(corr)

def autocorr2(x,lags):
    '''manualy compute, non partial'''

    mean=numpy.mean(x)
    var=numpy.var(x)
    xp=x-mean
    corr=[1. if l==0 else numpy.sum(xp[l:]*xp[:-l])/len(x)/var for l in lags]

    return numpy.array(corr)

def autocorr3(x,lags):
    '''fft, pad 0s, non partial'''

    n=len(x)
    # pad 0s to 2n-1
    ext_size=2*n-1
    # nearest power of 2
    fsize=2**numpy.ceil(numpy.log2(ext_size)).astype('int')

    xp=x-numpy.mean(x)
    var=numpy.var(x)

    # do fft and ifft
    cf=numpy.fft.fft(xp,fsize)
    sf=cf.conjugate()*cf
    corr=numpy.fft.ifft(sf).real
    corr=corr/var/n

    return corr[:len(lags)]

def autocorr4(x,lags):
    '''fft, don't pad 0s, non partial'''
    mean=x.mean()
    var=numpy.var(x)
    xp=x-mean

    cf=numpy.fft.fft(xp)
    sf=cf.conjugate()*cf
    corr=numpy.fft.ifft(sf).real/var/len(x)

    return corr[:len(lags)]

def autocorr5(x,lags):
    '''numpy.correlate, non partial'''
    mean=x.mean()
    var=numpy.var(x)
    xp=x-mean
    corr=numpy.correlate(xp,xp,'full')[len(x)-1:]/var/len(x)

    return corr[:len(lags)]


if __name__=='__main__':

    y=[28,28,26,19,16,24,26,24,24,29,29,27,31,26,38,23,13,14,28,19,19,\
            17,22,2,4,5,7,8,14,14,23]
    y=numpy.array(y).astype('float')

    lags=range(15)
    fig,ax=plt.subplots()

    for funcii, labelii in Zip([autocorr1, autocorr2, autocorr3, autocorr4,
        autocorr5], ['np.corrcoef, partial', 'manual, non-partial',
            'fft, pad 0s, non-partial', 'fft, no padding, non-partial',
            'np.correlate, non-partial']):

        cii=funcii(y,lags)
        print(labelii)
        print(cii)
        ax.plot(lags,cii,label=labelii)

    ax.set_xlabel('lag')
    ax.set_ylabel('correlation coefficient')
    ax.legend()
    plt.show()

Hier ist die Ausgabefigur:

 enter image description here

Wir sehen nicht alle 5 Zeilen, da sich 3 von ihnen überlappen (am Purpur). Die Überlappungen sind alle nicht partiellen Autokorrelationen. Dies liegt daran, dass Berechnungen aus den Signalverarbeitungsmethoden (np.correlate, FFT) nicht für jede Überlappung einen anderen Mittelwert/Standardwert berechnen.

Beachten Sie auch, dass das Ergebnis fft, no padding, non-partial (rote Linie) anders ist, da die Zeitreihe vor der FFT nicht mit 0s aufgefüllt wurde. Es handelt sich also um eine kreisförmige FFT. Ich kann nicht im Detail erklären warum, das habe ich von anderswo gelernt.

7
Jason

Ich verwende talib.CORREL für eine solche Autokorrelation. Ich vermute, Sie könnten dasselbe mit anderen Paketen machen:

def autocorrelate(x, period):

    # x is a deep indicator array 
    # period of sample and slices of comparison

    # oldest data (period of input array) may be nan; remove it
    x = x[-np.count_nonzero(~np.isnan(x)):]
    # subtract mean to normalize indicator
    x -= np.mean(x)
    # isolate the recent sample to be autocorrelated
    sample = x[-period:]
    # create slices of indicator data
    correls = []
    for n in range((len(x)-1), period, -1):
        alpha = period + n
        slices = (x[-alpha:])[:period]
        # compare each slice to the recent sample
        correls.append(ta.CORREL(slices, sample, period)[-1])
    # fill in zeros for sample overlap period of recent correlations    
    for n in range(period,0,-1):
        correls.append(0)
    # oldest data (autocorrelation period) will be nan; remove it
    correls = np.array(correls[-np.count_nonzero(~np.isnan(correls)):])      

    return correls

# CORRELATION OF BEST FIT
# the highest value correlation    
max_value = np.max(correls)
# index of the best correlation
max_index = np.argmax(correls)
1
litepresence

Mit der Fouriertransformation und dem Faltungssatz 

Die zeitliche Komplexität ist N * log (N)  

def autocorr1(x):
    r2=np.fft.ifft(np.abs(np.fft.fft(x))**2).real
    return r2[:len(x)//2]

Hier ist eine normalisierte und unvoreingenommene Version, es ist auch N * log (N)  

def autocorr2(x):
    r2=np.fft.ifft(np.abs(np.fft.fft(x))**2).real
    c=(r2/x.shape-np.mean(x)**2)/np.std(x)**2
    return c[:len(x)//2]

Die Methode, die von A. Levy bereitgestellt wird, funktioniert, aber ich habe sie in meinem PC getestet. Die Komplexität der Zeit scheint N * N zu sein.

def autocorr(x):
    result = numpy.correlate(x, x, mode='full')
    return result[result.size/2:]
0
wwwjjj

Eine einfache Lösung ohne Pandas:

import numpy as np

def auto_corrcoef(x):
   return np.corrcoef(x[1:-1], x[2:])[0,1]
0
dignitas

Zeichnen Sie die statistische Korrelation für eine Pandas-Datums-Zeitreihe auf:

import matplotlib.pyplot as plt

def plot_autocorr(returns, lags):
    autocorrelation = []
    for lag in range(lags+1):
        corr_lag = returns.corr(returns.shift(-lag)) 
        autocorrelation.append(corr_lag)
    plt.plot(range(lags+1), autocorrelation, '--o')
    plt.xticks(range(lags+1))
    return np.array(autocorrelation)
0

Eine Alternative zu numpy.correlate ist in statsmodels.tsa.stattools.acf () verfügbar. Dies ergibt eine kontinuierlich abnehmende Autokorrelationsfunktion wie die von OP beschriebene. Die Implementierung ist ziemlich einfach:

from statsmodels.tsa import stattools
# x = 1-D array
# Yield normalized autocorrelation function of number lags
autocorr = stattools.acf( x )

# Get autocorrelation coefficient at lag = 1
autocorr_coeff = autocorr[1]

Das Standardverhalten ist, bei 40 nlags anzuhalten, dies kann jedoch mit der Option nlag= für Ihre spezifische Anwendung angepasst werden. Am Ende der Seite befindet sich ein Zitat für die Statistiken hinter der Funktion .

0
fisherp

Ich bin ein Computerbiologe und als ich die Auto/Kreuz-Korrelationen zwischen Paaren von Zeitreihen stochastischer Prozesse berechnen musste, wurde mir klar, dass np.correlate nicht die Arbeit erledigte, die ich brauchte.

Tatsächlich scheint in np.correlate dieMittelung über alle möglichen Paare von Zeitpunktenim Abstand ???? zu fehlen.

So habe ich eine Funktion definiert, die das tut, was ich brauchte:

def autocross(x, y):
    c = np.correlate(x, y, "same")
    v = [c[i]/( len(x)-abs( i - (len(x)/2)  ) ) for i in range(len(c))]
    return v

Es scheint mir, dass keine der vorherigen Antworten diesen Fall von Auto/Kreuzkorrelation abdeckt: Ich hoffe, diese Antwort kann für jemanden nützlich sein, der an stochastischen Prozessen wie mir arbeitet.

0
Orso

Ich denke, die wirkliche Antwort auf die Frage des OP ist in diesem Auszug aus der Numpy.correlate-Dokumentation kurz zusammengefasst:

mode : {'valid', 'same', 'full'}, optional
    Refer to the `convolve` docstring.  Note that the default
    is `valid`, unlike `convolve`, which uses `full`.

Dies impliziert, dass die Numpy.correlate-Funktion einen Skalar zurückgibt, wenn derselbe Vektor für die beiden Eingabeargumente angegeben wird (d. H.

0
dbanas