使用 Python 估计自相关

Posted 2023-02-14

【中文标题】使用 Python 估计自相关

【英文标题】：Estimate Autocorrelation using Python

【发布时间】：2012-12-27 03:32:49

【问题描述】：

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

这是我想使用的估计自相关的方程（取自http://en.wikipedia.org/wiki/Autocorrelation，估计部分）：

在 python 中查找我的数据的估计自相关的最简单方法是什么？我可以使用类似于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()

【问题讨论】：

***.com/questions/4503325/…

我想更具体地谈谈估计的自相关方程。

另见：***.com/questions/643699/… 和 ***.com/questions/12269834/…

【参考方案1】：

我发现这只是稍微改变了预期的结果：

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 的自相关结果进行测试。

【参考方案2】：

我在最近一次编辑时编写的方法现在甚至比 scipy.statstools.acf 和 fft=True 还要快，直到样本量变得非常大。

错误分析如果您想调整偏差并获得高度准确的错误估计：查看我的代码 here，它实现了 Ulli Wolff 的 this paper(or original by UW in Matlab)

功能测试

a = correlatedData(n=10000) 来自here 的例程 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：我改变了我的例程并重新测试了 n=10000 和 n=20000 的 FFT

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

健全性检查

示例错误分析

这有点超出范围，但如果没有集成自相关时间或集成窗口计算，我不会费心重做这个数字。与错误的自相关在底部图中很清楚

【讨论】：

您的数据样本太小了，这怎么能比 fft 方法更快呢？您在这里将 n^2 与 nlogn 进行比较。

随意重复 - 代码都在上面。实际上，我确实声明 FFT 例程在非常大的 n 上成为更快的方法，所以我怀疑你只是扫描了第一句话。有可能是FFT方法每次调用都有大量的python开销和错误检查

抱歉，这样说只是误导，您只计算 40 个不同时移的相关性。通常，如果您有一个 200000 点的数据集，那么在性能方面有趣的部分是查看整个时间的相关函数。在这种情况下，您要处理 200000^2 数量级的操作，而 FFT 方法将处理 5*200000 左右。顺便说一句，我的幼稚实现与您的实现基本相同，在相同的输入上每个循环需要 18.1 毫秒...

首先，感谢您抽出宝贵时间查看此代码 - 我经常使用此代码，因此非常欢迎任何改进！随着相关时间的增加，样本量会减少，因此测量所有时间并没有那么有用，因为在较长时间内缺乏样本。您能否澄清您正在计时的方法并用另一种方法对其进行基准测试？ （不同的机器等）为了澄清，我同意 FFT 在规模上更快。然而，对于很多日常案例，样本n<10000 和我们查看一个小的自相关时间（典型的财务数据），这种方法似乎相当不错？

让我们continue this discussion in chat.

【参考方案3】：

我从 pandas 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

【参考方案4】：

statsmodels 包添加了一个内部使用np.correlate 的自相关函数（根据statsmodels 文档）。

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

【参考方案5】：

我认为这个特定的计算没有 NumPy 函数。以下是我的写法：

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

断言语句用于检查计算并记录其意图。

当您确信此函数按预期运行时，您可以注释掉 assert 语句，或使用 python -O 运行您的脚本。 （-O 标志告诉 Python 忽略断言语句。）

【讨论】：

谢谢。我认为这是计算自相关估计值的唯一方法。这可以很容易地用于查找使用x = loadtxt('fn.txt',unpack=True,usecols=[0]) 加载并绘制pylab.plot(autoCorr, t) 的数据的自相关吗？

是的，应该可以。也许试试pylab.plot(x, estimate_autocorrelation(x))...

如果数据是复值，则应在断言assert np.allclose(r, np.array([(x.conj()[:n-k]*x[-(n-k):]).sum() for k in range(n)]))中添加共轭