Jak mogę użyć numpy.correlate do wykonania autokorelacji?

Muszę zrobić autokorelację zbioru liczb, co, jak rozumiem, jest po prostu korelacją zbioru ze sobą.

Wypróbowałem to za pomocą funkcji korelacji numpy, ale nie wierzę w wynik, ponieważ prawie zawsze daje wektor, w którym pierwsza liczba nie jest największa, jak powinna.

Tak więc to pytanie to tak naprawdę dwa pytania:

1. Co dokładnie numpy.correlaterobi?
2. Jak mogę tego użyć (lub czegoś innego) do wykonania autokorelacji?

Aby odpowiedzieć na pierwsze pytanie, numpy.correlate(a, v, mode)należy wykonać splot az odwrotnością vi podać wyniki obcięte przez określony tryb. Definicji splotu c (t) = Ď _{-∞ <i <∞ i} v _{t + i,} gdzie -∞ <t <∞ pozwala na wyniki z -∞ do ∞, ale oczywiście nie można przechowywać nieskończenie długa szyk. Musi więc zostać przycięty i właśnie tam pojawia się tryb. Istnieją 3 różne tryby: pełny, taki sam i ważny:

• „full” zwraca wyniki dla każdego trybu t, w którym zarówno ai vpewne nakładanie.
• Tryb „ten sam” zwraca wynik o tej samej długości, co najkrótszy wektor ( alub v).
• „ważne” zwraca wyniki tylko wtedy, gdy tryb ai vcałkowicie zachodzą na siebie. Dokumentacja dla numpy.convolvedaje więcej szczegółów na temat trybów.

Na drugie pytanie, myślę, że numpy.correlate to daje autokorelacji, to jest po prostu daje trochę więcej, jak również. Autokorelacja służy do ustalenia, jak podobny jest do siebie sygnał lub funkcja przy określonej różnicy czasu. Przy różnicy czasu równej 0 autokorelacja powinna być najwyższa, ponieważ sygnał jest identyczny ze sobą, więc spodziewano się, że pierwszy element w tablicy wyników autokorelacji będzie największy. Jednak korelacja nie zaczyna się przy różnicy czasu równej 0. Zaczyna się od ujemnej różnicy czasu, zamyka się do 0, a następnie staje się dodatnia. Oznacza to, że spodziewałeś się:

autokorelacja (a) = ∑ _{-∞ <i <∞} a _i v _{t + i} gdzie 0 <= t <∞

Ale otrzymałeś:

autokorelacja (a) = ∑ _{-∞ <i <∞} a _i v _{t + i} gdzie -∞ <t <∞

To, co musisz zrobić, to wziąć ostatnią połowę wyniku korelacji i to powinna być autokorelacja, której szukasz. Prosta funkcja w Pythonie do zrobienia to:

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

Będziesz oczywiście potrzebował sprawdzenia błędów, aby upewnić się, że xfaktycznie jest to tablica 1-wymiarowa. Poza tym to wyjaśnienie prawdopodobnie nie jest najbardziej rygorystyczne matematycznie. Rzucałem się wokół nieskończoności, ponieważ definicja splotu używa ich, ale niekoniecznie dotyczy to autokorelacji. Zatem teoretyczna część tego wyjaśnienia może być nieco niepewna, ale miejmy nadzieję, że praktyczne wyniki będą pomocne. Te strony dotyczące autokorelacji są bardzo pomocne i mogą dać ci znacznie lepsze podstawy teoretyczne, jeśli nie masz nic przeciwko przebrnięciu przez notację i ciężkie koncepcje.

Autokorelacja występuje w dwóch wersjach: statystycznej i splotowej. Obaj robią to samo, z wyjątkiem małego szczegółu: Wersja statystyczna jest znormalizowana do przedziału [-1,1]. Oto przykład, jak wykonujesz statystyczny:

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

Użyj numpy.corrcoeffunkcji zamiast numpy.correlateobliczyć korelację statystyczną dla opóźnienia t:

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

Myślę, że są dwie rzeczy, które wprowadzają zamieszanie w tym temacie:

1. definicja statystyczna a przetwarzanie sygnału: jak zauważyli inni, w statystykach normalizujemy autokorelację do [-1,1].
2. częściowa vs nieczęściowa średnia / wariancja: gdy szeregi czasowe przesuwają się z opóźnieniem> 0, ich rozmiar nakładania się zawsze będzie <pierwotnej długości. Czy używamy średniej i standardowej wartości oryginalnej (nie częściowej), czy też zawsze obliczamy nową średnią i standardową przy użyciu ciągle zmieniającego się nakładania się (częściowego). (Prawdopodobnie jest na to formalny termin, ale na razie użyję „częściowego”).

Stworzyłem 5 funkcji, które obliczają autokorelację tablicy 1d, z rozróżnieniami częściowymi i nie częściowymi. Niektórzy używają wzoru ze statystyk, inni używają korelacji w sensie przetwarzania sygnału, co można również zrobić za pomocą FFT. Ale wszystkie wyniki są autokorelacjami w definicji statystyki , więc ilustrują, w jaki sposób są ze sobą powiązane. Kod poniżej:

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

Oto liczba wyjściowa:

Nie widzimy wszystkich 5 linii, ponieważ 3 z nich nakładają się (na fioletowo). Wszystkie nakładające się elementy są nie częściowymi autokorelacjami. Dzieje się tak, ponieważ obliczenia z metod przetwarzania sygnału ( np.correlateFFT) nie obliczają innej średniej / standardowej dla każdego nakładania się.

Zauważ również, że wynik fft, no padding, non-partial(czerwona linia) jest inny, ponieważ nie dopełnił on szeregu czasowego zerami przed wykonaniem FFT, więc jest to cykliczna FFT. Nie potrafię szczegółowo wyjaśnić, dlaczego, tego nauczyłem się gdzie indziej.

Ponieważ właśnie napotkałem ten sam problem, chciałbym podzielić się z wami kilkoma wierszami kodu. W rzeczywistości istnieje już kilka dość podobnych postów dotyczących autokorelacji w przepływie stosu. Jeśli zdefiniujesz autokorelację jako a(x, L) = sum(k=0,N-L-1)((xk-xbar)*(x(k+L)-xbar))/sum(k=0,N-1)((xk-xbar)**2)[jest to definicja podana w funkcji a_correlate IDL i zgadza się z tym, co widzę w odpowiedzi 2 na pytanie nr 12269834 ], to poniższe wydaje się dawać prawidłowe wyniki:

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

Jak widać, przetestowałem to z krzywą grzechu i jednolitym losowym rozkładem i oba wyniki wyglądają tak, jakbym się ich spodziewał. Zauważ, że użyłem mode="same"zamiast tego, mode="full"co zrobili inni.

Twoje pytanie 1 zostało już obszernie omówione w kilku doskonałych odpowiedziach tutaj.

Pomyślałem, że podzielę się z wami kilkoma wierszami kodu, które pozwolą obliczyć autokorelację sygnału na podstawie wyłącznie matematycznych właściwości autokorelacji. Oznacza to, że autokorelację można obliczyć w następujący sposób:

1. odejmij średnią od sygnału i uzyskaj nieobciążony sygnał

2. obliczyć transformatę Fouriera nieobciążonego sygnału

3. obliczyć gęstość widmową mocy sygnału, biorąc normę kwadratową każdej wartości transformaty Fouriera nieobciążonego sygnału

4. obliczyć odwrotną transformatę Fouriera gęstości widmowej mocy

5. znormalizuj odwrotną transformatę Fouriera gęstości widmowej mocy przez sumę kwadratów nieobciążonego sygnału i weź tylko połowę wynikowego wektora

Kod służący do tego jest następujący:

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)

Jestem biologiem obliczeniowym i kiedy musiałem obliczyć korelacje auto / krzyżowe między parami szeregów czasowych procesów stochastycznych, zdałem sobie sprawę, że np.correlatenie spełniało to moich oczekiwań.

Rzeczywiście, wydaje się, że brakuje np.correlatew nim uśredniania wszystkich możliwych par punktów czasowych na odległość 𝜏.

Oto jak zdefiniowałem funkcję wykonującą to, czego potrzebowałem:

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

Wydaje mi się, że żadna z poprzednich odpowiedzi nie obejmuje tego przypadku korelacji automatycznej / krzyżowej: mam nadzieję, że ta odpowiedź może być przydatna dla kogoś, kto pracuje nad procesami stochastycznymi, tak jak ja.

Używam talib.CORREL do autokorelacji w ten sposób, podejrzewam, że możesz zrobić to samo z innymi pakietami:

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)

Zastosowanie transformacji Fouriera i twierdzenia o splotach

Złożoność czasowa to N * log (N)

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

Oto znormalizowana i nieobciążona wersja, jest to również 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]

Metoda dostarczona przez A. Levy'ego działa, ale przetestowałem ją na swoim komputerze, jej złożoność czasowa wydaje się być N * N

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

Alternatywa dla numpy.correlate jest dostępna w statsmodels.tsa.stattools.acf () . Daje to stale malejącą funkcję autokorelacji, taką jak ta opisana przez OP. Wdrożenie jest dość proste:

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]

Domyślnym zachowaniem jest zatrzymanie przy 40 nlagach, ale można to dostosować za pomocą nlag=opcji dla konkretnej aplikacji. Na dole strony znajduje się cytat ze statystykami funkcji .

Myślę, że prawdziwa odpowiedź na pytanie PO jest zwięźle zawarta w tym fragmencie dokumentacji Numpy.correlate:

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

Oznacza to, że gdy jest używana bez definicji „trybu”, funkcja Numpy.correlate zwróci wartość skalarną, gdy otrzyma ten sam wektor dla swoich dwóch argumentów wejściowych (tj. - gdy zostanie użyta do wykonania autokorelacji).

Proste rozwiązanie bez pand:

import numpy as np

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

Wykreśl autokorelację statystyczną, biorąc pod uwagę dane pandy Seria zwrotów:

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)