¿Cómo puedo utilizar numpy.se correlacionan para hacer de autocorrelación?

https://stackoverflow.com/questions/643699

22-07-2019
|

Pregunta

Necesito hacer auto-correlación de un conjunto de números, que según entiendo es sólo la correlación de la serie con el mismo.

Yo lo he probado usando numpy la función de correlacionar, pero no creo que el resultado, como casi siempre da un vector donde el primer número es no el más grande, como debe ser.

Así, esta pregunta es, en realidad, dos preguntas:

¿Qué es exactamente numpy.correlate haciendo?
Cómo puedo usar (o algo más) para hacer auto-correlación?

Solución

Para responder a su primera pregunta, numpy.correlate(a, v, mode) se realiza la convolución de a con el reverso de v y dando los resultados recortadas por el especificado modo.El definición de la convolución, C(t)=∑ _{-∞ < yo < ∞} un_yov_t+i donde -∞ < t < ∞, permite resultados desde -∞ a ∞, pero obviamente no puede almacenar un infinitamente larga de la matriz.Así que tiene que ser recortado, y que es donde entra en juego el modo.Hay 3 modos diferentes:completo, mismo, y válido:

"full" en el modo de devolución de los resultados para cada t donde tanto a y v tiene algunas coincidencias.
"del mismo modo" devuelve un resultado con la misma longitud que el vector más corto (a o v).
"válida" a modo de devolución de los resultados sólo cuando a y v completamente se superponen unos a otros.El documentación para numpy.convolve da más detalles sobre los modos.

Para tu segunda pregunta, creo que numpy.correlate es dándole la autocorrelación, es sólo darle un poco más así.La autocorrelación se utiliza para encontrar la forma similar a una señal, o una función, es en sí mismo una cierta diferencia de tiempo.En una diferencia de tiempo de 0, la función de auto-correlación debe ser mayor debido a que la señal es idéntica a sí misma, por lo que se espera que el primer elemento de la autocorrelación de la matriz de resultado sería el mayor.Sin embargo, la correlación no es partir de una diferencia de tiempo de 0.Se inicia en un negativo de la diferencia de tiempo, cierra a 0, y luego se va positivos.Es decir, que estaban esperando:

de autocorrelación(a) = ∑ _{-∞ < yo < ∞} un_yov_t+i donde 0 <= t < ∞

Pero lo que recibí fue:

de autocorrelación(a) = ∑ _{-∞ < yo < ∞} un_yov_t+i donde -∞ < t < ∞

Lo que usted necesita hacer es tomar la última mitad de su correlación resultado, y que debe ser la autocorrelación que usted está buscando.Una simple función de python para hacer, que serían:

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

Usted, por supuesto, la necesidad de la comprobación de errores para asegurarse de que x en realidad es un 1-d de la matriz.Además, esta explicación probablemente no es el más matemáticamente rigurosa.He estado lanzando alrededor de los infinitos debido a que la definición de la convolución los usa, pero que no se aplican necesariamente a la autocorrelación.Así, el teórico parte de la explicación puede ser un poco flojo, pero esperemos que los resultados prácticos son útiles. Estos páginas en autocorrelación son bastante útiles, y pueden darle una mejor base teórica si no te importa que vadear a través de la notación y pesado conceptos.

Otros consejos

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

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

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

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

Como me encontré con el mismo problema, me gustaría compartir algunas líneas de código con usted. De hecho, hay varias publicaciones bastante similares sobre autocorrelación en stackoverflow por ahora. Si 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 la función a_correlate de IDL y está de acuerdo con lo que veo en la respuesta 2 de la pregunta # 12269834 ], entonces lo 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 puede ver, he probado esto con una curva sin y una distribución aleatoria uniforme, y parece que ambos resultados los esperaría. Tenga en cuenta que usé mode="same" en lugar de mode="full" como lo hicieron los demás.

Creo que hay 2 cosas que agregan confusión a este tema:

estadística v.s. definición del procesamiento de la señal: como otros han señalado, en estadística normalizamos la autocorrelación en [-1,1].
parcial v.s. media / varianza no parcial: cuando la serie de tiempo cambia a un retraso > 0, su tamaño de superposición siempre será < longitud original ¿Usamos la media y estándar del original (no parcial), o siempre calculamos una nueva media y estándar usando la superposición siempre cambiante (parcial) hace la diferencia. (Probablemente haya un término formal para esto, pero voy a usar & Quot; parcial & Quot; por ahora).

He creado 5 funciones que calculan la autocorrelación de una matriz 1d, con v.s. parciales. distinciones no parciales Algunos usan fórmulas de estadísticas, otros usan correlación en el sentido del procesamiento de la señal, que también se puede hacer a través de FFT. Pero todos los resultados son correlaciones automáticas en la definición de estadísticas , por lo que ilustran cómo están vinculados entre sí. Código a continuación:

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

Aquí está la cifra de salida:

No vemos las 5 líneas porque 3 de ellas se superponen (en el morado). Las superposiciones son todas correlaciones automáticas no parciales. Esto se debe a que los cálculos de los métodos de procesamiento de señales (np.correlate, FFT) no calculan una media / estándar diferente para cada superposición.

También tenga en cuenta que el resultado fft, no padding, non-partial (línea roja) es diferente, porque no rellena la serie de tiempo con 0s antes de hacer FFT, por lo que es FFT circular. No puedo explicar en detalle por qué, eso es lo que aprendí de otra parte.

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

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

resta la media de la señal y obtiene una señal imparcial
calcula la transformada de Fourier de la señal imparcial
calcula la densidad espectral de potencia de la señal, tomando la norma cuadrada de cada valor de la transformada de Fourier de la señal imparcial
calcula la transformada inversa de Fourier de la densidad espectral de potencia
normaliza la transformada inversa de Fourier de la densidad espectral de potencia por la suma de los cuadrados de la señal imparcial, y toma solo la mitad del vector resultante

El código para hacer esto es el 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)

Uso talib.CORREL para autocorrelación como esta, 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)

Usando la transformación de Fourier y el teorema de convolución

La complejidad del tiempo es N * log (N)

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

Aquí hay una versión normalizada e imparcial, también es 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]

El método proporcionado por A. Levy funciona, pero lo probé en mi PC, su complejidad temporal parece ser N*N

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

Creo que la respuesta real a la pregunta del OP está contenida sucintamente en este extracto de la documentación de 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 usa sin definición de 'modo', la función Numpy.correlate devolverá un escalar, cuando se le dé el mismo vector para sus dos argumentos de entrada (es decir, cuando se usa para realizar la autocorrelación).

Una solución simple sin pandas:

import numpy as np

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

Trace la autocorrelación estadística dada una serie de retornos de fecha y hora pandas:

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)

Soy biólogo computacional, y cuando tuve que calcular las correlaciones automáticas / cruzadas entre parejas de series temporales de procesos estocásticos, me di cuenta de que np.correlate no estaba haciendo el trabajo que necesitaba.

De hecho, lo que parece faltar en <=> es el promedio de todos los posibles pares de puntos de tiempo a distancia & # 120591;.

Así es como definí una función haciendo lo que necesitaba:

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

Me parece que ninguna de las respuestas anteriores cubre esta instancia de correlación automática / cruzada: espero que esta respuesta pueda ser útil para alguien que trabaja en procesos estocásticos como yo.

Una alternativa a numpy.correlate está disponible en statsmodels. tsa.stattools.acf () . Esto produce una función de autocorrelación continuamente decreciente como la descrita por OP. Implementarlo es bastante 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]

El comportamiento predeterminado es detenerse en 40 nlags, pero esto se puede ajustar con la opción nlag= para su aplicación específica. Hay una cita en la parte inferior de la página para las estadísticas detrás de la función .

Licenciado bajo: CC-BY-SA con atribución

No afiliado a StackOverflow