Comment utiliser numpy.correlate pour effectuer une autocorrélation?

https://stackoverflow.com/questions/643699

22-07-2019
|

Question

Je dois effectuer une auto-corrélation d'un ensemble de nombres, qui, si j'ai bien compris, n'est que la corrélation de l'ensemble avec lui-même.

Je l'ai essayé en utilisant la fonction de corrélation de numpy, mais je ne crois pas au résultat, car il donne presque toujours un vecteur dont le premier nombre n'est pas le plus grand, comme il se doit .

Donc, cette question est vraiment deux questions:

Que fait exactement numpy.correlate ?
Comment puis-je l'utiliser (ou quelque chose d'autre) pour effectuer une corrélation automatique?

La solution

Pour répondre à votre première question, numpy.correlate(a, v, mode) effectue la convolution de a avec l'inverse de v et donne les résultats tronqués selon le mode spécifié. La définition de la convolution , C (t) = & # 8721; _{- & # 8734; < je < & # 8734;} a _i v _{t + i} où - & # 8734; < t < & # 8734 ;, permet d'obtenir les résultats de - & # 8734; à & # 8734 ;, mais vous ne pouvez évidemment pas stocker un tableau infiniment long. Il doit donc être coupé, et c’est là que le mode entre en jeu. Il existe 3 modes différents: complet, identique, & Amp; valide:

" complet " Le mode renvoie des résultats pour chaque t où numpy.convolve et numpy.correlate présentent des chevauchements.
" même " Le mode renvoie un résultat de la même longueur que le vecteur le plus court (x ou <=>).
" valide " Le mode renvoie les résultats uniquement lorsque <=> et <=> se chevauchent complètement. La documentation de <=> donne plus de détails sur la modes.

Pour votre deuxième question, je pense que <=> est vous donne l'autocorrélation, elle vous en donne juste un peu plus. L'autocorrélation est utilisée pour rechercher dans quelle mesure un signal, ou une fonction, est similaire à lui-même à une certaine différence de temps. À une différence de temps de 0, l'auto-corrélation doit être la plus grande, car le signal est identique à lui-même. Vous vous attendiez donc à ce que le premier élément du tableau de résultats d'autocorrélation soit le plus grand. Toutefois, la corrélation ne commence pas à une différence de temps de 0. Elle commence à une différence de temps négative, se ferme à 0, puis devient positive. C’est-à-dire que vous attendiez:

autocorrélation (a) = & # 8721; _{- & # 8734; < je < & # 8734;} a _i v _{t + i} où 0 < = t < & # 8734;

Mais ce que vous avez eu était:

autocorrélation (a) = & # 8721; _{- & # 8734; < je < & # 8734;} a _i v _{t + i} où - & # 8734; < t < & # 8734;

Ce que vous devez faire, c'est prendre la dernière moitié de votre résultat de corrélation, et cela devrait être l'autocorrélation que vous recherchez. Une simple fonction python à faire serait:

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

Vous aurez bien sûr besoin d’une vérification d’erreur pour vous assurer que <=> est en réalité un tableau à une dimension. En outre, cette explication n’est probablement pas la plus rigoureuse mathématiquement. J'ai jeté autour d'infinis parce que la définition de convolution les utilise, mais cela ne s'applique pas nécessairement à l'autocorrélation. Donc, la partie théorique de cette explication peut être un peu confuse, mais j'espère que les résultats pratiques sont utiles. Ces les pages sur l’autocorrélation sont très utiles et peuvent vous fournir une base théorique bien meilleure si cela ne vous dérange pas de parcourir la notation et les concepts lourds.

Autres conseils

L’auto-corrélation existe en deux versions: statistique et convolution. Ils font tous les deux la même chose, à l'exception d'un petit détail: la version statistique est normalisée sur l'intervalle [-1,1]. Voici un exemple de la méthode statistique utilisée:

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

Utilisez plutôt la fonction numpy.corrcoef de numpy.correlate pour calculer la corrélation statistique de un décalage de t:

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

Comme je viens de rencontrer le même problème, je voudrais partager quelques lignes de code avec vous. En fait, il existe plusieurs publications assez similaires sur l'autocorrélation dans stackoverflow. Si vous définissez l'autocorrélation comme a(x, L) = sum(k=0,N-L-1)((xk-xbar)*(x(k+L)-xbar))/sum(k=0,N-1)((xk-xbar)**2) [c'est la définition donnée dans la fonction a_correlate d'IDL et elle est conforme à ce que je vois dans la réponse 2 de la question # 12269834 ], alors ce qui suit semble donner les résultats corrects:

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

Comme vous le voyez, j’ai testé cela avec une courbe de péché et une distribution aléatoire uniforme, et les deux résultats ressemblent à ceux que j’attendais. Notez que j'ai utilisé mode="same" au lieu de mode="full" comme les autres.

Je pense qu'il y a 2 choses qui ajoutent de la confusion à ce sujet:

statistiques v.s. définition du traitement du signal: comme d'autres l'ont souligné, dans les statistiques, nous normalisons l'auto-corrélation dans [-1,1].
v. partiel moyenne / variance non partielle: lorsque la série temporelle se décale avec un décalage > 0, leur taille de chevauchement sera toujours < longueur d'origine. Utilisons-nous les valeurs moyenne et standard de l'original (non-partiel) ou calculons-nous toujours une nouvelle moyenne et standard en utilisant le chevauchement en constante évolution (partiel) qui fait la différence? (Il y a probablement un terme officiel pour cela, mais je vais utiliser & "Partiel &"; Pour l'instant).

J'ai créé 5 fonctions qui calculent l'auto-corrélation d'un tableau 1d, avec v.s. partiel. distinctions non partielles. Certains utilisent des formules issues de statistiques, d’autres utilisent corréler dans le sens du traitement du signal, ce qui peut également se faire via FFT. Cependant, tous les résultats sont des auto-corrélations dans la définition des statistiques . Ils illustrent donc la manière dont ils sont liés. Code ci-dessous:

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

Voici le chiffre de sortie:

Nous ne voyons pas les 5 lignes car 3 d'entre elles se chevauchent (en violet). Les chevauchements sont tous des corrélations automatiques non partielles. En effet, les calculs issus des méthodes de traitement du signal (np.correlate, FFT) ne calculent pas une moyenne / std différente pour chaque chevauchement.

Notez également que le résultat fft, no padding, non-partial (ligne rouge) est différent, car il n'a pas rempli les séries de temps avec des 0 avant de réaliser la FFT, il est donc circulaire. Je ne peux pas expliquer en détail pourquoi, c'est ce que j'ai appris ailleurs.

Votre question 1 a déjà été longuement discutée dans plusieurs excellentes réponses ici.

J'ai pensé partager avec vous quelques lignes de code vous permettant de calculer l'autocorrélation d'un signal en vous basant uniquement sur les propriétés mathématiques de l'autocorrélation. Autrement dit, l’autocorrélation peut être calculée de la manière suivante:

soustrayez la moyenne du signal et obtenez un signal non biaisé
calculez la transformée de Fourier du signal non biaisé
calcule la densité spectrale de puissance du signal en prenant la norme carrée de chaque valeur de la transformée de Fourier du signal non polarisé
calcule la transformée de Fourier inverse de la densité spectrale de puissance
normalisez la transformée de Fourier inverse de la densité spectrale de puissance par la somme des carrés du signal non polarisé et ne prenez que la moitié du vecteur résultant

Le code à suivre est le suivant:

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)

J'utilise talib.CORREL pour l'autocorrélation de cette manière, je suppose que vous pourriez faire la même chose avec d'autres packages:

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)

Utilisation de la transformation de Fourier et du théorème de convolution

La complexité temporelle est N * log (N)

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

Voici une version normalisée et non biaisée, il s'agit également de 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]

La méthode fournie par A. Levy fonctionne, mais je l’ai testée sur mon PC. Sa complexité temporelle semble être N * N

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

Je pense que la véritable réponse à la question du PO figure succinctement dans cet extrait de la documentation de Numpy.correlate:

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

Cela implique que, lorsqu'elle n'est utilisée avec aucune définition de "mode", la fonction Numpy.correlate renverra un scalaire si le même vecteur lui est attribué pour ses deux arguments d'entrée (c'est-à-dire - lorsqu'elle est utilisée pour effectuer une autocorrélation).

Une solution simple sans pandas:

import numpy as np

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

Tracez l'autocorrélation statistique en fonction d'un datas de pandas. Série de retours:

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)

Je suis un biologiste informaticien et lorsque je devais calculer les corrélations auto / croisées entre des couples de séries chronologiques de processus stochastiques, je compris que np.correlate ne faisait pas le travail dont j'avais besoin.

En effet, il semble manquer à <=> la moyenne sur tous les couples possibles de points de temps à distance & # 120591;.

Voici comment j'ai défini une fonction faisant ce dont j'avais besoin:

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

Il me semble qu'aucune des réponses précédentes ne couvre cet exemple de corrélation automatique / croisée: espérons que cette réponse sera utile à quelqu'un qui travaille sur des processus stochastiques comme moi.

Une alternative à numpy.correlate est disponible dans statsmodels. tsa.stattools.acf () . Ceci donne une fonction d'autocorrélation décroissante de façon constante comme celle décrite par OP. La mise en œuvre est assez simple:

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]

Le comportement par défaut est d’arrêter à 40 nlags, mais cela peut être ajusté avec l’option nlag= pour votre application spécifique. Il y a une citation au bas de la page pour la statistiques derrière la fonction .

Licencié sous: CC-BY-SA avec attribution

Non affilié à StackOverflow