估计自相关使用Python(Estimate Autocorrelation using Pytho

我想在下面示出的信号执行自相关。 两个连续点之间的时间是2.5ms的（或400Hz的重复率）。

这是用于估计autoacrrelation等式，我想用（来自http://en.wikipedia.org/wiki/Autocorrelation ，部分估计）：

什么是找到我的数据的估计自相关蟒蛇的最简单的方法？ 是否有类似的东西numpy.correlate ，我可以使用？

或者我应该只是计算的均值和方差？

---

编辑：

从帮助unutbu ，我已经写了：

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:

我不认为有这种特别的计算一个NumPy的功能。 这是我会怎么写：

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

断言语句是有这两个检查计算并记录它的意图。

当你相信这个功能是预期行为，你可以注释掉的assert语句，或运行脚本python -O 。 （该-O标志告诉Python来忽略断言语句。）

Answer 2:

我把从大熊猫autocorrelation_plot（）函数的代码的一部分。 我检查的答案与R和值精确匹配。

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

Answer 3:

所述statsmodels包补充说，在内部使用一个自相关函数np.correlate （根据statsmodels文档）。

请参阅： http://statsmodels.sourceforge.net/stable/generated/statsmodels.tsa.stattools.acf.html#statsmodels.tsa.stattools.acf

Answer 4:

我写我的最新编辑的方法现在快甚至比scipy.statstools.acf与fft=True ，直到试样尺寸变得非常大。

错误分析如果要调整的偏见和得到高度精确的错误估计：看看我的代码在这里 ，它实现这一纸 （由乌利·沃尔夫在由UW或原Matlab ）

功能测试

• a = correlatedData(n=10000)是从找到的例程这里
• gamma()是从同一个地方correlated_data()
• acorr()是我下面的函数
• estimated_autocorrelation在另一个答案发现
• acf()是from statsmodels.tsa.stattools import acf

计时

%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

编辑...我检查再次保持l=40 ，改变n=10000至n=200000样本的FFT方法开始变得有点牵引和statsmodels FFT实现只是边缘吧...（顺序是一样的）

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

编辑2：我改变了我的常规和重新测试主场迎战FFT为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

履行

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的加速可以在下面实现的。 你必须小心，只通过op_samples=a.copy()因为它会修改数组a由a-=mean否则：

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

完整性检查

例如误差分析

这是一个有点出入的范围，但我不能被人打扰重做没有集成的自相关时间或积分窗口计算的身影。 有错误的自相关是在底部情节清晰

Answer 5:

我发现这有只略有变化的预期结果：

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

测试针对Excel的自相关结果。

文章来源: Estimate Autocorrelation using Python