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

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

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. 

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)

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:

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.

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)

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

Eine einfache Lösung ohne Pandas:

import numpy as np

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

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)

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 .

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.

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.