Como posso usar numpy.correlate fazer autocorrelação?

https://stackoverflow.com/questions/643699

22-07-2019
|

Pergunta

Eu preciso fazer auto-correlação de um conjunto de números, que como eu entendo que é apenas a correlação do conjunto com a própria.

Já fiz usando a função de correlação de numpy, mas eu não acredito que o resultado, já que quase sempre dá um vector em que o primeiro número é não o maior, como deveria ser .

Assim, esta questão é realmente duas perguntas:

O que é exatamente numpy.correlate fazendo?
Como posso usá-lo (ou outra coisa) para fazer auto-correlação?

Solução

Para responder à sua primeira pergunta, numpy.correlate(a, v, mode) está realizando a convolução de a com o reverso de v e dando os resultados cortadas pelo modo especificado. A definição de convolução , C (t) = S _{-8 a _i v _{t + i}, onde -8}

"cheio" de modo retorna resultados para cada t onde ambos a e v têm alguma sobreposição.

"mesmo" modo retorna um resultado com o mesmo comprimento que o vector mais curto (a ou v).

"válido" modo retorna resultados apenas quando a e v completamente sobrepõem uns aos outros. A documentação para numpy.convolve dá mais detalhes sobre os modos.

Para sua segunda pergunta, eu acho numpy.correlate é dando-lhe a autocorrelação, está apenas dando-lhe um pouco mais também. A autocorrelação é usado para encontrar sinal de uma forma semelhante, ou função, é em si a uma certa diferença de tempo. Com uma diferença de tempo de 0, a autocorrelação deve ser o mais alto porque o sinal é idêntico a si mesmo, para que você esperar que o primeiro elemento na matriz resultado autocorrelação seria o maior. No entanto, a correlação não é a partir de uma diferença de tempo de 0. Ela começa com uma diferença de tempo negativo, fecha-se a 0, e, em seguida, vai positivo. Ou seja, você estava esperando:

autocorrelação (a) = S _{-8 _i v _{t + i} em que 0 <= t <8}

Mas o que você conseguiu foi:

autocorrelação (a) = S _{-8 _i v _{t + i} onde -8}
O que você precisa fazer é tomar a última metade do seu resultado de correlação, e que deve ser a autocorrelação que você está procurando. A função de python simples de fazer isso seria:

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

Você vai, naturalmente, erro necessidade de verificar para se certificar de que x é realmente uma matriz 1-d. Além disso, esta explicação provavelmente não é o mais matematicamente rigorosa. Eu estive jogando em torno de infinitos porque a definição de convolução usa-los, mas isso não se aplica necessariamente para autocorrelação. Assim, a parte teórica desta explicação pode ser um pouco instável, mas espero que os resultados práticos são úteis. Estes páginas em autocorrelação são muito úteis, e pode dar-lhe um fundo muito melhor teórica, se você não se importa vasculhar a notação e os conceitos pesados.

Outras dicas

Auto-correlação vem em duas versões: estatística e convolução. Ambos fazem a mesma, exceto por um pequeno detalhe: A versão estatística é normalizado para estar no intervalo [-1,1]. Aqui está um exemplo de como você faz a estatística:

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

Use a numpy.corrcoef funcionar em vez de < a href = "https://docs.scipy.org/doc/numpy/reference/generated/numpy.correlate.html" rel = "nofollow noreferrer"> numpy.correlate para calcular a correlação estatística para um atraso de t :

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

Como eu apenas corri para o mesmo problema, eu gostaria de compartilhar algumas linhas de código com você. Na verdade, existem vários posts em vez semelhantes sobre autocorrelação em stackoverflow por agora. Se você definir a autocorrelação como a(x, L) = sum(k=0,N-L-1)((xk-xbar)*(x(k+L)-xbar))/sum(k=0,N-1)((xk-xbar)**2) [esta é a definição dada em função a_correlate do IDL e concorda com o que vejo em resposta 2 de pergunta # 12269834 ], então o seguinte parece dar os resultados corretos:

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

Como você vê eu testei isso com uma curva de pecado e uma distribuição aleatória uniforme, e ambos os resultados olhar como eu esperaria. Note que eu usei mode="same" vez de mode="full" como os outros fizeram.

Eu acho que existem 2 coisas que acrescentam confusão a este tópico:

V.S. estatística definição de processamento de sinal: como outros já apontaram, nas estatísticas normalizamos auto-correlação em [-1,1].

V.S. parciais não parcial média / variância: quando os TimeSeries turnos em um lag> 0, seu tamanho sobreposição será sempre

Eu criei 5 funções que auto-correlação de computação de uma matriz 1d, com V.S. parcial distinções não parciais. Alguns fórmula uso de estatísticas, algum correlato uso no sentido de processamento de sinal, que também pode ser feito através de FFT. Mas todos os resultados são auto-correlações nos estatísticas definição, para que eles mostram como eles estão ligados uns aos outros. Código abaixo:

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

Aqui está a figura de saída:

Nós não vemos todas as 5 linhas por 3 deles se sobrepõem (no roxo). As sobreposições são todas as auto-correlações não parciais. Isto é porque os cálculos dos métodos de processamento de sinal (np.correlate, FFT) não calcular um / std significativo diferente para cada sobreposição.

Além disso, note que o fft, no padding, non-partial (linha vermelha) resultado é diferente, porque ele não fez pad os timeseries com 0s antes de fazer FFT, por isso é FFT circular. Eu não posso explicar em detalhes o porquê, isso é o que eu aprendi de outros lugares.

A sua questão 1 já foi amplamente discutido em várias respostas excelentes aqui.

Eu pensei para compartilhar com vocês algumas linhas de código que permitem calcular a autocorrelação de um sinal apenas com base nas propriedades matemáticas de autocorrelação. Ou seja, a autocorrelação pode ser calculado da seguinte forma:

subtrair a média do sinal e obter um sinal imparcial

calcular a transformada de Fourier do sinal imparcial

calcular a densidade espectral de energia do sinal, levando a norma quadrado de cada valor da transformada de Fourier do sinal imparcial

calcular a transformada de Fourier inversa da densidade espectral de potência

normalizar a transformada de Fourier inversa da densidade espectral de potência pela soma dos quadrados do sinal imparcial, e ter apenas metade do vector que resulta

O código para fazer isso é o seguinte:

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)

Eu uso talib.CORREL para autocorrelação assim, eu suspeito que você poderia fazer o mesmo com outros pacotes:

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)

Usando transformação de Fourier e a convolução teorema

O complexicity tempo é N * log (N)

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

Aqui está uma versão normalizada e imparcial, é também 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]

O método fornecido por A. Levy funciona, mas eu testei no meu PC, seu complexicity tempo parece ser N * N

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

Eu acho que a verdadeira resposta à pergunta do OP é sucintamente contida neste trecho da documentação Numpy.correlate:

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

Isto implica que, quando utilizado com nenhuma definição 'mode', a função Numpy.correlate irá retornar um escalar, quando dado o mesmo vector de seus dois argumentos de entrada (ou seja - quando usado para executar autocorrelação).

Uma solução simples, sem pandas:

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

Traçar a autocorrelação estatística dada uma pandas datetime série de retornos:

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)

Eu sou um biólogo computacional, e quando eu tinha para calcular as auto / correlações cruzadas entre casais de séries temporais de processos estocásticos percebi que np.correlate não estava fazendo o trabalho que eu precisava.

Na verdade, o que parece estar faltando np.correlate é o média sobre todos os possíveis pares de pontos de tempo em t distância.

Aqui está como eu definida uma função de fazer o que eu precisava:

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

Parece-me nenhuma das respostas anteriores cobrir esta instância de auto / correlação cruzada:. Esperamos que esta resposta pode ser útil a alguém trabalhando em processos estocásticos como eu

Uma alternativa para numpy.correlate está disponível em statsmodels. tsa.stattools.acf () . Isso resulta uma função de autocorrelação continuamente decrescente como o descrito por OP. Implementá-lo é bastante simples:

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]

O comportamento padrão é parar em 40 nlags, mas isso pode ser ajustado com a opção nlag= para sua aplicação específica. Há uma citação na parte inferior da página para os estatísticas por trás da função .

Licenciado em: CC-BY-SA com atribuição

Não afiliado a StackOverflow