2009-03-13 40 views
69

Necesito hacer la autocorrelación de un conjunto de números, que según entiendo, es solo la correlación del conjunto consigo mismo.¿Cómo puedo usar numpy.correlate para hacer la autocorrelación?

Lo he intentado usando la función de correlación de numpy, pero no creo el resultado, ya que casi siempre da un vector donde el primer número es no el más grande, como debería ser.

Por lo tanto, esta cuestión es en realidad dos preguntas:

  1. ¿Qué es exactamente está haciendo numpy.correlate?
  2. ¿Cómo puedo usarlo (u otra cosa) para hacer una autocorrelación?
+0

Consulte también: http://stackoverflow.com/questions/12269834/is-there-any-numpy-autocorrellation-function-with-standardized-output para obtener información acerca de autocorrelación normalizada. – amcnabb

Respuesta

79

para responder a su primera pregunta, numpy.correlate(a, v, mode) está realizando la convolución de a con el reverso de v y dando los resultados recortadas por el modo especificado. El definition of convolution, C (t) = Σ -∞ ∞ < i < un i v t + i donde -∞ < t < ∞, permite obtener resultados de -∞ a ∞, pero es obvio que no puede almacenar una matriz infinitamente larga. . Así que tiene que ser cortada, y que es donde entra en juego el modo en el Hay 3 modos diferentes: por completo, mismo, & válida:

  • "Full" devuelve resultados para cada t donde tanto a y v tienen algunos se superponen
  • "modo mismo" devuelve un resultado con la misma longitud que el vector más corto (a o v).
  • modo "válido" devuelve resultados solo cuando a y v se superponen por completo. El documentation para numpy.convolve brinda más detalles sobre los modos.

a su segunda pregunta, creo que es numpy.correlate que le da la autocorrelación, sólo se le está dando un poco más también. La autocorrelación se usa para encontrar cuán similar es una señal, o función, a sí misma en una determinada diferencia de tiempo. En una diferencia de tiempo de 0, la autocorrelación debería ser la más alta porque la señal es idéntica a sí misma, por lo que esperaba que el primer elemento en la matriz de resultados de autocorrelación fuera el más grande. Sin embargo, la correlación no comienza en una diferencia de tiempo de 0. Comienza en una diferencia de tiempo negativa, se cierra a 0 y luego se vuelve positiva.Es decir, que se esperaba:

autocorrelación (a) = Σ -∞ ∞ < i < un i v t + i donde 0 = t < < ∞

Pero lo que tengo era:

autocorrelación (a) = Σ -∞ < i < ∞ un i v t + i donde -∞ ∞ < t <

Lo que hay que hacer es tomar la última mitad de su resultado de correlación, y que debe ser la autocorrelación que busca. Una función de Python sencilla de hacerlo sería:

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

Usted, por supuesto, la necesidad de comprobación de errores para asegurarse de que x es en realidad una matriz 1-d. Además, esta explicación probablemente no sea la más matemáticamente rigurosa. He estado lanzando infinitos porque la definición de convolución los usa, pero eso no se aplica necesariamente a la autocorrelación. Por lo tanto, la parte teórica de esta explicación puede ser un poco inestable, pero es de esperar que los resultados prácticos sean útiles. Thesepages en la autocorrelación son muy útiles, y pueden proporcionarle una base teórica mucho mejor si no le importa vadear la notación y los conceptos pesados.

+4

En compilaciones actuales de numpy, el modo 'mismo' se puede especificar para lograr exactamente lo que propuso A. Levy. El cuerpo de la función podría leer 'return numpy.correlate (x, x, mode = 'same')' –

+10

@DavidZwicker pero los resultados son diferentes. 'np.correlate (x, x, mode = 'lleno') [len (x) // 2:]! = np.correlate (x, x, mode = 'igual')'. Por ejemplo, 'x = [1,2,3,1,2]; np.correlate (x, x, mode = 'full'); '{' >>> array ([2, 5, 11, 13, 19, 13, 11, 5, 2]) '}' np.correlate (x, x, mode = 'mismo'); '{' >>> matriz ([11, 13, 19, 13, 11]) '}. El correcto es: 'np.correlate (x, x, mode = 'full') [len (x) -1:];' {'>>> array ([19, 13, 11, 5, 2]) '} ver ** el primer elemento ** es ** el más grande **. – Developer

+7

Tenga en cuenta que esta respuesta proporciona la autocorrelación no normalizada. – amcnabb

12

La autocorrelación viene en dos versiones: estadística y convolución. Ambos hacen lo mismo, excepto por un pequeño detalle: el primero está normalizado para estar en el intervalo [-1,1]. Aquí está un ejemplo de cómo se hace la estadística:

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

Desea '' numpy.corrcoef [x: -i], x [i:]) [0,1] '' en la segunda línea, ya que el valor de retorno de '' corrcoef'' es una matriz de 2x2 – luispedro

+0

¿Cuál es la diferencia? entre las autocorrelaciones estadísticas y las convolucionales? –

9

Como lo acabo de encontrar con el mismo problema, me gustaría compartir unas pocas líneas de código con usted. De hecho, hay varias publicaciones bastante similares sobre autocorrelación en stackoverflow por el momento. Si se define la autocorrelación 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 es la definición dada en función a_correlate de IDL y está de acuerdo con lo que veo en respuesta 2 de pregunta #12269834], a continuación, la siguiente parece dar los resultados correctos:

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 verá que he probado esto con una curva de sincrón y una distribución aleatoria uniforme, y ambos resultados parecen como los esperaría. Tenga en cuenta que utilicé mode="same" en lugar de mode="full" como lo hicieron los demás.

13

Uso de la función en lugar de numpy.correlatenumpy.corrcoef para calcular la correlación estadística para un retraso de t:

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

¿Los "coeficientes de correlación" no se refieren a la autocorrelación utilizada en el procesamiento de la señal y no a la autocorrelación utilizada en las estadísticas? https://en.wikipedia.org/wiki/Autocorrelation#Signal_processing –

+0

@DanielPendergast No estoy tan familiarizado con el procesamiento de la señal. De los documentos numpy: "Devolver los coeficientes de correlación producto-momento de Pearson". ¿Es esa la versión de procesamiento de señal? –

0

creo que la verdadera respuesta a la pregunta del OP está contenida sucintamente en este extracto de la documentación Numpy.correlate :

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

Esto implica que, cuando se utiliza con ninguna definición 'modo', la función Numpy.correlate devolverá un escalar, cuando se les da el mismo vector para sus dos argumentos de entrada (es decir, - cuando utilizado para realizar la autocorrelación).

1

Yo uso talib.COEF.DE.CORREL la autocorrelación como esto, sospecho que podría hacer lo mismo con otros paquetes:

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

Su pregunta 1 ya ha sido ampliamente discutido en varias excelentes respuestas aquí.

Pensé compartir con ustedes algunas líneas de código que le permiten calcular la autocorrelación de una señal basada únicamente en las propiedades matemáticas de la autocorrelación. Es decir, la autocorrelación puede ser calculada de la siguiente manera:

  1. restar la media de la señal y obtener una señal imparcial

  2. calcular la transformada de Fourier de la señal imparcial

  3. de cómputo la densidad espectral de potencia de la señal, al tomar la norma cuadrada de cada valor de la transformada de Fourier de la señal no sesgada

  4. calcular la transformada de Fourier inversa de la densidad

  5. espectral de potencia normalizar la transformada de Fourier inversa de la densidad espectral de potencia por la suma de los cuadrados de la señal no sesgada, y tomar sólo la mitad del vector resultante

El código para hacer esto es la siguiente:

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) 
+0

¿es posible que haya algún problema con esto? No puedo hacer coincidir sus resultados con otras funciones de correlación automática. La función se ve similar pero parece algo aplastada. – pindakaas

+0

@pindakaas ¿podría ser más específico? por favor brinde información sobre cuáles son las discrepancias que encuentra, con qué funciones. – Ruggero

+0

¿Por qué no usar 'p = np.abs (f)'? – dylnan

Cuestiones relacionadas