¿Hay alguna función de autocorrellación numpy con salida estandarizada?

Seguí el consejo de definir la función de autocorrelación en otro post:

def autocorr(x): result = np.correlate(x, x, mode = 'full') maxcorr = np.argmax(result) #print 'maximum = ', result[maxcorr] result = result / result[maxcorr] # <=== normalization return result[result.size/2:] 

sin embargo, el valor máximo no era “1.0”. por lo tanto, introduje la línea etiquetada con “<=== normalización"

Probé la función con el conjunto de datos del “Análisis de series de tiempo” (Cuadro – Jenkins), capítulo 2. Esperaba obtener un resultado como el de la fig. 2.7 en ese libro. Sin embargo conseguí lo siguiente:

introduzca la descripción de la imagen aquí

¿Alguien tiene una explicación para este extraño comportamiento inesperado de autocorrelación?

Adición (2012-09-07):

Me metí en Python – progtwigción e hice lo siguiente:

 from ClimateUtilities import * import phys # # the above imports are from RTPierrehumbert's book "principles of planetary # climate" # and the homepage of that book at "cambridge University press" ... they mostly # define the # class "Curve()" used in the below section which is not necessary in order to solve # my # numpy-problem ... :) # import numpy as np; import scipy.spatial.distance; # functions to be defined ... : # # def autocorr(x): result = np.correlate(x, x, mode = 'full') maxcorr = np.argmax(result) # print 'maximum = ', result[maxcorr] result = result / result[maxcorr] # return result[result.size/2:] ## # second try ... "Box and Jenkins" chapter 2.1 Autocorrelation Properties # of stationary models ## # from table 2.1 I get: s1 = np.array([47,64,23,71,38,64,55,41,59,48,71,35,57,40,58,44,\ 80,55,37,74,51,57,50,60,45,57,50,45,25,59,50,71,56,74,50,58,45,\ 54,36,54,48,55,45,57,50,62,44,64,43,52,38,59,\ 55,41,53,49,34,35,54,45,68,38,50,\ 60,39,59,40,57,54,23],dtype=float); # alternatively in order to test: s2 = np.array([47,64,23,71,38,64,55,41,59,48,71]) ##################################################################################3 # according to BJ, ch.2 ###################################################################################3 print '*************************************************' global s1short, meanshort, stdShort, s1dev, s1shX, s1shXk s1short = s1 #s1short = s2 # for testing take s2 meanshort = s1short.mean() stdShort = s1short.std() s1dev = s1short - meanshort #print 's1short = \n', s1short, '\nmeanshort = ', meanshort, '\ns1deviation = \n',\ # s1dev, \ # '\nstdShort = ', stdShort s1sh_len = s1short.size s1shX = np.arange(1,s1sh_len + 1) #print 'Len = ', s1sh_len, '\nx-value = ', s1shX ########################################################## # c0 to be computed ... ########################################################## sumY = 0 kk = 1 for ii in s1shX: #print 'ii-1 = ',ii-1, if ii > s1sh_len: break sumY += s1dev[ii-1]*s1dev[ii-1] #print 'sumY = ',sumY, 's1dev**2 = ', s1dev[ii-1]*s1dev[ii-1] c0 = sumY / s1sh_len print 'c0 = ', c0 ########################################################## # now compute autocorrelation ########################################################## auCorr = [] s1shXk = s1shX lenS1 = s1sh_len nn = 1 # factor by which lenS1 should be divided in order # to reduce computation length ... 1, 2, 3, 4 # should not exceed 4 #print 's1shX = ',s1shX for kk in s1shXk: sumY = 0 for ii in s1shX: #print 'ii-1 = ',ii-1, ' kk = ', kk, 'kk+ii-1 = ', kk+ii-1 if ii >= s1sh_len or ii + kk - 1>=s1sh_len/nn: break sumY += s1dev[ii-1]*s1dev[ii+kk-1] #print sumY, s1dev[ii-1], '*', s1dev[ii+kk-1] auCorrElement = sumY / s1sh_len auCorrElement = auCorrElement / c0 #print 'sum = ', sumY, ' element = ', auCorrElement auCorr.append(auCorrElement) #print '', auCorr # #manipulate s1shX # s1shX = s1shXk[:lenS1-kk] #print 's1shX = ',s1shX #print 'AutoCorr = \n', auCorr ######################################################### # # first 15 of above Values are consistent with # Box-Jenkins "Time Series Analysis", p.34 Table 2.2 # ######################################################### s1sh_sdt = s1dev.std() # Standardabweichung short #print '\ns1sh_std = ', s1sh_sdt print '#########################################' # "Curve()" is a class from RTP ClimateUtilities.py c2 = Curve() s1shXfloat = np.ndarray(shape=(1,lenS1),dtype=float) s1shXfloat = s1shXk # to make floating point from integer # might be not necessary #print 'test plotting ... ', s1shXk, s1shXfloat c2.addCurve(s1shXfloat) c2.addCurve(auCorr, '', 'Autocorr') c2.PlotTitle = 'Autokorrelation' w2 = plot(c2) ########################################################## # # now try function "autocorr(arr)" and plot it # ########################################################## auCorr = autocorr(s1short) c3 = Curve() c3.addCurve( s1shXfloat ) c3.addCurve( auCorr, '', 'Autocorr' ) c3.PlotTitle = 'Autocorr with "autocorr"' w3 = plot(c3) # # well that should it be! # 

Entonces, su problema con su bash inicial es que no restó el promedio de su señal. El siguiente código debería funcionar:

 timeseries = (your data here) mean = np.mean(timeseries) timeseries -= np.mean(timeseries) autocorr_f = np.correlate(timeseries, timeseries, mode='full') temp = autocorr_f[autocorr_f.size/2:]/autocorr_f[autocorr_f.size/2] iact.append(sum(autocorr_f[autocorr_f.size/2:]/autocorr_f[autocorr_f.size/2])) 

En mi ejemplo, temp es la variable que te interesa; Es la función de autocorrelación integrada hacia adelante. Si desea el tiempo de autocorrelación integrado, está interesado en iact .

No estoy seguro de cuál es el problema.

La autocorrelación de un vector x debe ser 1 en el retardo 0, ya que es solo la norma L2 al cuadrado dividida por sí misma, es decir, dot(x, x) / dot(x, x) == 1 .

En general, para cualquier retraso i, j in Z, where i != j la autocorrelación escalada de la unidad es dot(shift(x, i), shift(x, j)) / dot(x, x) donde shift(y, n) es una función que desplaza el vector y por n puntos temporales y Z es el conjunto de enteros ya que estamos hablando de la implementación (en teoría, los retrasos pueden estar en el conjunto de números reales).

Obtengo 1.0 como máximo con el siguiente código (comience en la línea de comandos como $ ipython --pylab ), como se esperaba:

 In[1]: n = 1000 In[2]: x = randn(n) In[3]: xc = correlate(x, x, mode='full') In[4]: xc /= xc[xc.argmax()] In[5]: xchalf = xc[xc.size / 2:] In[6]: xchalf_max = xchalf.max() In[7]: print xchalf_max Out[1]: 1.0 

El único momento en que la autocorrelación 0 no es igual a 1 es cuando x es la señal cero (todos ceros).

La respuesta a su pregunta es : no , no hay ninguna función NumPy que realice automáticamente la estandarización para usted.

Además, incluso si lo hiciera, tendría que compararlo con la salida esperada y si puede decir “Sí, esto realizó la estandarización correctamente”, entonces asumo que usted sabe cómo implementarlo usted mismo.

Voy a sugerir que podría ser el caso de que haya implementado su algoritmo incorrectamente, aunque no puedo estar seguro ya que no estoy familiarizado con él.