Come posso usare numpy.correlate per eseguire l'autocorrelazione?

https://stackoverflow.com/questions/643699

22-07-2019
|

Domanda

Devo fare l'auto-correlazione di un insieme di numeri, che a quanto ho capito è solo la correlazione dell'insieme con se stesso.

L'ho provato usando la funzione correl di numpy, ma non credo al risultato, dato che dà quasi sempre un vettore in cui il primo numero è non il più grande, come dovrebbe essere .

Quindi, questa domanda è in realtà due domande:

Che cosa sta esattamente facendo numpy.correlate ?
Come posso usarlo (o qualcos'altro) per fare l'auto-correlazione?

Soluzione

Per rispondere alla tua prima domanda, numpy.correlate(a, v, mode) sta eseguendo la convoluzione di a con il contrario di v e dando i risultati ritagliati dalla modalità specificata. La definizione di convoluzione , C (t) = & # 8721; & - _{# 8734; lt &; i < & # 8734;} a _i v _{t + i} dove - & # 8734; lt &; t < & # 8734 ;, consente di ottenere risultati da - & # 8734; su & # 8734 ;, ma ovviamente non è possibile memorizzare un array infinitamente lungo. Quindi deve essere tagliato, ed è qui che entra in gioco la modalità. Ci sono 3 diverse modalità: full, same, & Amp; valido:

quot; & Pieno quot; mode restituisce risultati per ogni t in cui sia numpy.convolve che numpy.correlate presentano alcune sovrapposizioni.
quot; & Stessa quot; mode restituisce un risultato con la stessa lunghezza del vettore più corto (x o <=>).
quot; quot valido &; mode restituisce risultati solo quando <=> e <=> si sovrappongono completamente. La documentazione per <=> fornisce ulteriori dettagli sulla le modalità.

Per la tua seconda domanda, penso che <=> ti stia dandoti l'autocorrelazione, ma ti sta dando anche un po 'di più. L'autocorrelazione viene utilizzata per scoprire quanto un segnale, o una funzione, sia simile a se stesso in una certa differenza di tempo. Con una differenza temporale di 0, l'auto-correlazione dovrebbe essere la più alta perché il segnale è identico a se stesso, quindi ci si aspettava che il primo elemento nella matrice dei risultati di autocorrelazione fosse il massimo. Tuttavia, la correlazione non inizia con una differenza oraria di 0. Inizia con una differenza oraria negativa, si chiude su 0 e quindi diventa positiva. Cioè, ti aspettavi:

autocorrelazione (a) = & # 8721; & - _{# 8734; lt &; i < & # 8734;} a _i v _{t + i} dove 0 < = t < & # 8734;

Ma quello che hai ottenuto è stato:

autocorrelazione (a) = & # 8721; & - _{# 8734; lt &; i < & # 8734;} a _i v _{t + i} dove - & # 8734; lt &; t < & # 8734;

Quello che devi fare è prendere l'ultima metà del risultato di correlazione, e quella dovrebbe essere l'autocorrelazione che stai cercando. Una semplice funzione Python per farlo sarebbe:

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

Naturalmente avrai bisogno di un controllo degli errori per assicurarti che <=> sia effettivamente un array 1-d. Inoltre, questa spiegazione probabilmente non è la più matematicamente rigorosa. Ho gettato infinito perché la definizione di convoluzione li usa, ma ciò non si applica necessariamente all'autocorrelazione. Quindi, la parte teorica di questa spiegazione può essere leggermente traballante, ma speriamo che i risultati pratici siano utili. Questi pagine sull'autocorrelazione sono piuttosto utili e possono darti un background teorico molto migliore se non ti dispiace andare in giro con la notazione e concetti pesanti.

Altri suggerimenti

L'auto-correlazione è disponibile in due versioni: statistica e convoluzione. Entrambi fanno lo stesso, tranne per un piccolo dettaglio: la versione statistica è normalizzata per essere sull'intervallo [-1,1]. Ecco un esempio di come si fa quello statistico:

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

Utilizza invece la funzione numpy.corrcoef di numpy.correlate per calcolare la correlazione statistica per un ritardo di t:

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

Dato che ho riscontrato lo stesso problema, vorrei condividere alcune righe di codice con te. In effetti, ci sono molti post piuttosto simili sull'autocorrelazione nello stackoverflow ormai. Se definisci l'autocorrelazione come a(x, L) = sum(k=0,N-L-1)((xk-xbar)*(x(k+L)-xbar))/sum(k=0,N-1)((xk-xbar)**2) [questa è la definizione fornita nella funzione a_correlate di IDL e concorda con ciò che vedo nella risposta 2 della domanda # 12269834 ], quindi quanto segue sembra dare i risultati corretti:

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

Come vedi, l'ho testato con una curva sin e una distribuzione casuale uniforme, ed entrambi i risultati sembrano come me li aspettassi. Nota che ho usato mode="same" invece di mode="full" come hanno fatto gli altri.

Penso che ci siano 2 cose che aggiungono confusione a questo argomento:

statistiche v.s. definizione di elaborazione del segnale: come altri hanno sottolineato, nelle statistiche normalizziamo l'auto-correlazione in [-1,1].
v.s. parziali media / varianza non parziale: quando la serie temporale si sposta con un ritardo > 0, la loro dimensione di sovrapposizione sarà sempre < lunghezza originale. Usiamo la media e lo std dell'originale (non parziale) o calcoliamo sempre una nuova media e lo std usando la sovrapposizione sempre variabile (parziale) fa la differenza. (Probabilmente esiste un termine formale per questo, ma userò & Quot; parziale & Quot; per ora).

Ho creato 5 funzioni che calcolano l'auto-correlazione di un array 1d, con v.s. parziali distinzioni non parziali. Alcuni usano la formula delle statistiche, altri usano la correlazione nel senso dell'elaborazione del segnale, che può anche essere fatto tramite FFT. Ma tutti i risultati sono auto-correlazioni nella definizione statistiche , quindi illustrano come sono collegati tra loro. Codice seguente:

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

Ecco la figura di output:

Non vediamo tutte e 5 le righe perché 3 di esse si sovrappongono (in viola). Le sovrapposizioni sono tutte autocorrelazioni non parziali. Questo perché i calcoli dei metodi di elaborazione del segnale (np.correlate, FFT) non calcolano una media / std diversa per ogni sovrapposizione.

Nota anche che il risultato fft, no padding, non-partial (linea rossa) è diverso, perché non ha riempito la serie con 0s prima di fare FFT, quindi è FFT circolare. Non posso spiegare in dettaglio perché, è quello che ho imparato altrove.

La tua domanda 1 è già stata ampiamente discussa in diverse eccellenti risposte qui.

Ho pensato di condividere con te alcune righe di codice che ti consentono di calcolare l'autocorrelazione di un segnale basato solo sulle proprietà matematiche dell'autocorrelazione. Cioè, l'autocorrelazione può essere calcolata nel modo seguente:

sottrae la media dal segnale e ottiene un segnale imparziale
calcola la trasformata di Fourier del segnale imparziale
calcola la densità spettrale di potenza del segnale, prendendo la norma quadrata di ciascun valore della trasformata di Fourier del segnale imparziale
calcola la trasformata inversa di Fourier della densità spettrale di potenza
normalizza la trasformata inversa di Fourier della densità spettrale di potenza per la somma dei quadrati del segnale imparziale e prende solo la metà del vettore risultante

Il codice per farlo è il seguente:

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)

Uso talib.CORREL per autocorrelazione come questa, sospetto che potresti fare lo stesso con altri pacchetti:

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)

Uso della trasformazione di Fourier e del teorema di convoluzione

La complessità temporale è N * log (N)

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

Ecco una versione normalizzata e imparziale, è anche 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]

Il metodo fornito da A. Levy funziona, ma l'ho provato sul mio PC, la sua complessità temporale sembra essere N*N

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

Penso che la vera risposta alla domanda del PO sia sinteticamente contenuta in questo estratto della documentazione di Numpy.correlate:

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

Ciò implica che, se usata senza una definizione 'mode', la funzione Numpy.correlate restituirà uno scalare, quando viene dato lo stesso vettore per i suoi due argomenti di input (cioè - quando usato per eseguire l'auto-correlazione).

Una soluzione semplice senza panda:

import numpy as np

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

Tracciare l'autocorrelazione statistica data una serie di resi datetime panda:

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)

Sono un biologo computazionale e quando ho dovuto calcolare le correlazioni automatiche / incrociate tra coppie di serie temporali di processi stocastici mi sono reso conto che np.correlate non stava facendo il lavoro di cui avevo bisogno.

In effetti, ciò che sembra mancare da <=> è la media su tutte le possibili coppie di punti temporali a distanza & # 120591;.

Ecco come ho definito una funzione che fa ciò di cui avevo bisogno:

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

Mi sembra che nessuna delle risposte precedenti copra questa istanza di auto / correlazione incrociata: spero che questa risposta possa essere utile a qualcuno che lavora su processi stocastici come me.

Un'alternativa a numpy.correlate è disponibile in statsmodels. tsa.stattools.acf () . Ciò produce una funzione di autocorrelazione in continuo calo come quella descritta da OP. L'implementazione è abbastanza semplice:

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]

Il comportamento predefinito è fermarsi a 40 nlags, ma questo può essere regolato con l'opzione nlag= per l'applicazione specifica. C'è una citazione in fondo alla pagina per le dietro la funzione .

Autorizzato sotto: CC-BY-SA insieme a attribuzione

Non affiliato a StackOverflow