python - iir - scipy help

Estimar la autocorrelación usando Python

Question

Me gustaría realizar la autocorrelación en la señal que se muestra a continuación. El tiempo entre dos puntos consecutivos es 2.5ms (o una tasa de repetición de 400Hz).

Esta es la ecuación para estimar la autoacrrelación que me gustaría usar (tomada de http://en.wikipedia.org/wiki/Autocorrelation , sección Estimación)

¿Cuál es el método más simple para encontrar la autocorrelación estimada de mis datos en python? ¿Hay algo similar a numpy.correlate que pueda usar?

¿O debería simplemente calcular la media y la varianza?

Editar:

Con la ayuda de unutbu , he escrito:

from numpy import *
import numpy as N
import pylab as P

fn = ''data.txt''
x = loadtxt(fn,unpack=True,usecols=[1])
time = loadtxt(fn,unpack=True,usecols=[0]) 

def estimated_autocorrelation(x):
    n = len(x)
    variance = x.var()
    x = x-x.mean()
    r = N.correlate(x, x, mode = ''full'')[-n:]
    #assert N.allclose(r, N.array([(x[:n-k]*x[-(n-k):]).sum() for k in range(n)]))
    result = r/(variance*(N.arange(n, 0, -1)))
    return result

P.plot(time,estimated_autocorrelation(x))
P.xlabel(''time (s)'')
P.ylabel(''autocorrelation'')
P.show()

Answer 1

El método que escribí a partir de mi última edición ahora es más rápido que incluso scipy.statstools.acf con fft=True hasta que el tamaño de la muestra es muy grande.

Análisis de errores Si desea ajustar los sesgos y obtener estimaciones de error altamente precisas: vea mi código here que implementa este documento de Ulli Wolff ( u original de UW en Matlab )

Funciones probadas

• a = correlatedData(n=10000) es de una rutina que se encuentra here
• gamma() es del mismo lugar que correlated_data()
• acorr() es mi función a continuación
• estimated_autocorrelation se encuentra en otra respuesta
• acf() es from statsmodels.tsa.stattools import acf

Tiempos

%timeit a0, junk, junk = gamma(a, f=0)                            # puwr.py
%timeit a1 = [acorr(a, m, i) for i in range(l)]                   # my own
%timeit a2 = acf(a)                                               # statstools
%timeit a3 = estimated_autocorrelation(a)                         # numpy
%timeit a4 = acf(a, fft=True)                                     # stats FFT

## -- End pasted text --
100 loops, best of 3: 7.18 ms per loop
100 loops, best of 3: 2.15 ms per loop
10 loops, best of 3: 88.3 ms per loop
10 loops, best of 3: 87.6 ms per loop
100 loops, best of 3: 3.33 ms per loop

Editar ... Verifiqué nuevamente manteniendo l=40 y cambiando n=10000 a n=200000 muestras. Los métodos de FFT comienzan a obtener un poco de tracción y las implementaciones de statsmodels simplemente lo limitan ... (el orden es el mismo)

## -- End pasted text --
10 loops, best of 3: 86.2 ms per loop
10 loops, best of 3: 69.5 ms per loop
1 loops, best of 3: 16.2 s per loop
1 loops, best of 3: 16.3 s per loop
10 loops, best of 3: 52.3 ms per loop

Edición 2: cambié mi rutina y volví a probar frente a la FFT para n=10000 n=20000

a = correlatedData(n=200000); b=correlatedData(n=10000)
m = a.mean(); rng = np.arange(40); mb = b.mean()
%timeit a1 = map(lambda t:acorr(a, m, t), rng)
%timeit a1 = map(lambda t:acorr.acorr(b, mb, t), rng)
%timeit a4 = acf(a, fft=True)
%timeit a4 = acf(b, fft=True)

10 loops, best of 3: 73.3 ms per loop   # acorr below
100 loops, best of 3: 2.37 ms per loop  # acorr below
10 loops, best of 3: 79.2 ms per loop   # statstools with FFT
100 loops, best of 3: 2.69 ms per loop # statstools with FFT

Implementación

def acorr(op_samples, mean, separation, norm = 1):
    """autocorrelation of a measured operator with optional normalisation
    the autocorrelation is measured over the 0th axis

    Required Inputs
        op_samples  :: np.ndarray :: the operator samples
        mean        :: float :: the mean of the operator
        separation  :: int :: the separation between HMC steps
        norm        :: float :: the autocorrelation with separation=0
    """
    return ((op_samples[:op_samples.size-separation] - mean)*(op_samples[separation:]- mean)).ravel().mean() / norm

4x aceleración se puede lograr a continuación. Debe tener cuidado de solo pasar op_samples=a.copy() ya que modificará la matriz a por a-=mean contrario:

op_samples -= mean
return (op_samples[:op_samples.size-separation]*op_samples[separation:]).ravel().mean() / norm

Prueba de cordura

Ejemplo de análisis de error

Esto está un poco fuera de alcance, pero no puedo molestarme en rehacer la figura sin el cálculo del tiempo de autocorrelación integrado o de la ventana de integración. Las autocorrelaciones con errores son claras en la gráfica inferior.

Answer 2

El paquete statsmodels agrega una función de autocorrelación que utiliza internamente np.correlate (de acuerdo con la documentación de statsmodels ).

Consulte: http://statsmodels.sourceforge.net/stable/generated/statsmodels.tsa.stattools.acf.html#statsmodels.tsa.stattools.acf

Answer 3

Encontré que esto obtuvo los resultados esperados con solo un ligero cambio:

def estimated_autocorrelation(x):
    n = len(x)
    variance = x.var()
    x = x-x.mean()
    r = N.correlate(x, x, mode = ''full'')
    result = r/(variance*n)
    return result

Pruebas contra los resultados de autocorrelación de Excel.

Answer 4

No creo que haya una función NumPy para este cálculo en particular. Así es como lo escribiría:

def estimated_autocorrelation(x):
    """
    http://.com/q/14297012/190597
    http://en.wikipedia.org/wiki/Autocorrelation#Estimation
    """
    n = len(x)
    variance = x.var()
    x = x-x.mean()
    r = np.correlate(x, x, mode = ''full'')[-n:]
    assert np.allclose(r, np.array([(x[:n-k]*x[-(n-k):]).sum() for k in range(n)]))
    result = r/(variance*(np.arange(n, 0, -1)))
    return result

La declaración de afirmación está ahí para verificar el cálculo y documentar su intención.

Cuando esté seguro de que esta función se comporta como se espera, puede comentar la declaración de afirmación o ejecutar su script con python -O . (La -O le dice a Python que ignore las aseveraciones).

Answer 5

Tomé parte del código de la función pandas autocorrelation_plot (). Comprobé las respuestas con R y los valores coinciden exactamente.

import numpy
def acf(series):
    n = len(series)
    data = numpy.asarray(series)
    mean = numpy.mean(data)
    c0 = numpy.sum((data - mean) ** 2) / float(n)

    def r(h):
        acf_lag = ((data[:n - h] - mean) * (data[h:] - mean)).sum() / float(n) / c0
        return round(acf_lag, 3)
    x = numpy.arange(n) # Avoiding lag 0 calculation
    acf_coeffs = map(r, x)
    return acf_coeffs