webentwicklung-frage-antwort-db.com.de

Schätzen Sie die Autokorrelation mit Python

Ich möchte eine Autokorrelation mit dem unten gezeigten Signal durchführen. Die Zeit zwischen zwei aufeinander folgenden Punkten beträgt 2,5 ms (oder eine Wiederholrate von 400 Hz).

enter image description here

Dies ist die Gleichung zum Schätzen der Autoacrrelation, die ich verwenden möchte (entnommen aus http://en.wikipedia.org/wiki/Autocorrelation , Abschnitt Estimation):

enter image description here

Was ist die einfachste Methode, um die geschätzte Autokorrelation meiner Daten in Python zu ermitteln? Gibt es etwas Ähnliches zu numpy.correlate, das ich verwenden kann? 

Oder soll ich nur den Mittelwert und die Varianz berechnen?


Bearbeiten:

Mit der Hilfe von unutbu habe ich geschrieben:

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

Ich glaube nicht, dass es für diese Berechnung eine NumPy-Funktion gibt. So würde ich es schreiben:

def estimated_autocorrelation(x):
    """
    http://stackoverflow.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

Die assert-Anweisung dient dazu, die Berechnung zu überprüfen und ihre Absicht zu dokumentieren. 

Wenn Sie sicher sind, dass sich diese Funktion wie erwartet verhält, können Sie die assert-Anweisung auskommentieren oder Ihr Skript mit python -O ausführen. (Das Flag -O weist Python an, Assert-Anweisungen zu ignorieren.)

29
unutbu

Ich habe einen Teil des Codes aus der Pandas autocorrelation_plot () - Funktion entnommen. Ich habe die Antworten mit R überprüft und die Werte stimmen genau überein.

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
16

Das Paket statsmodels fügt eine Autokorrelationsfunktion hinzu, die intern np.correlate verwendet (gemäß der statsmodels-Dokumentation).

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

11
hamogu

Die Methode, die ich in meiner letzten Bearbeitung geschrieben habe, ist jetzt schneller als scipy.statstools.acf mit fft=True, bis die Stichprobengröße sehr groß wird.

Fehleranalyse Wenn Sie Abweichungen anpassen und hochpräzise Fehlerschätzungen erhalten möchten: Sehen Sie sich meinen Code here an, der dieses Dokument von Ulli Wolff ( oder Original von UW) implementiert in Matlab )

Getestete Funktionen

  • a = correlatedData(n=10000) stammt aus einer Routine hier
  • gamma() ist vom selben Ort wie correlated_data() 
  • acorr() ist meine Funktion unten
  • estimated_autocorrelation wird in einer anderen Antwort gefunden
  • acf() ist von from statsmodels.tsa.stattools import acf

Timings

%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

Bearbeiten ... Ich überprüfte erneut, l=40 und n=10000 in n=200000-Samples ändern. Die FFT-Methoden beginnen, ein wenig Traktion zu erhalten, und statsmodels fft-Implementierung fasst sie einfach zusammen ... (Reihenfolge ist gleich)

## -- 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

Edit 2: Ich habe meine Routine geändert und die FFT für n=10000 und n=20000 erneut getestet.

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

Implementierung

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 speedup kann unten erreicht werden. Sie müssen nur op_samples=a.copy() übergeben, da das Array a um a-=mean geändert wird. Andernfalls:

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

Gesundheitsüberprüfung

enter image description here

Beispiel Fehleranalyse

Dies ist ein wenig außerhalb des Bereichs, aber ich kann mich nicht darum kümmern, die Figur ohne die integrierte Autokorrelationszeit oder die Berechnung des Integrationsfensters erneut auszuführen. Die fehlerbehafteten Autokorrelationen sind in der unteren Grafik sichtbar enter image description here

7

Ich fand, dass dies mit einer geringfügigen Änderung die erwarteten Ergebnisse erzielte:

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

Testen gegen die Autokorrelationsergebnisse von Excel.

1
Corin Dawson