如何使用numpy.correlate进行自相关？

Question

如何使用numpy.correlate进行自相关？

pythonmathnumpynumerical-methods

133

我需要对一组数字进行自相关，据我所知这只是该集合与其自身的相关性。

我已经尝试使用NumPy的correlate函数，但我不相信结果，因为它几乎总是给出一个向量，其中第一个数字不是最大值，而它应该是。

因此，这个问题实际上有两个问题：

numpy.correlate 到底在做什么？
我如何使用它（或其他东西）来进行自相关？

- Ben

请参阅以下链接以获取有关标准化自相关的NumPy自相关函数的信息：https://dev59.com/0WjWa4cB1Zd3GeqPmw_0。 - amcnabb

14个回答

37

自相关有两种版本：统计和卷积。它们都做相同的事情，除了一个小细节：统计版本被标准化到区间[-1,1]。以下是如何执行统计版本的示例:

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

- jonathf

9

第二行中你希望返回numpy.corrcoef[x:-i], x[i:])[0,1]，因为corrcoef的返回值是一个2x2的矩阵。 - luispedro

统计自相关和卷积自相关有什么区别？ - Daniel Node.js

2

@DanielPendergast：第二句话回答了这个问题：_它们都是一样的，除了一个小细节：前者 [统计学] 被规范化到区间[-1，1]上_。 - n1k31t4

1

统计自相关意味着您可以建立时间序列与先前协变量（例如参数回归）在某个时间点的统计关系：$E(y_t|y_{t-k}) = fn(y_{t-k}, \beta)$ 和 $Var(y_t|y_{t-k}) = f(...) $。由于期望和方差本质上是统计量，因此得名。卷积自相关通常用于信号处理，例如通过使用卷积（滑动）窗口（如正弦波）在信号上进行平滑/过滤，并在每个点上逐步相乘和求和。 - Sid

28

我认为这个主题会有两个混淆的地方：

统计学与信号处理定义：正如其他人所指出的，在统计学中，我们将自相关归一化为[-1,1]。
部分和非部分均值/方差：当时间序列在滞后>0时发生变化时，它们的重叠大小总是<原始长度。我们使用原始（非部分）的平均值和标准差，还是始终使用不断变化的重叠（部分）计算新的平均值和标准差会有所不同。（可能有一个正式术语，但我现在要使用“部分”）。

我创建了5个函数来计算1d数组的自相关，包括部分和非部分的区别。有些使用统计学公式，有些使用信号处理中的相关性，也可以通过FFT完成。但所有结果都是自相关的统计学定义，因此它们说明了它们之间的联系。以下是代码：

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

这是输出的图像：

我们看不到所有5条线，因为它们中的3条在紫色处重叠。这些重叠都是非部分自相关。这是因为信号处理方法（np.correlate、FFT）的计算没有为每个重叠计算不同的均值/标准差。

此外，请注意fft、无填充、非部分（红线）的结果是不同的，因为它在进行FFT之前没有用0对时序数据进行填充，所以它是循环FFT。我不能详细解释为什么，这是我从其他地方学到的。

- Jason

有没有一种使用基于FFT的方法执行“部分”自相关的方式？ - Daniel Whettam

@DanielWhettam 我不这么认为。 - Jason

在autocorr2中，为什么不用len(x) - l来除，而是用len(x)？现在这样会低估相关性，对吧？好吧，你可能会说只要l << len(x)就可以了，这是一个鲁棒估计所必需的条件。 - Mahé

@Mahé len(x)是“非部分”的规范，因此使用完整长度。 - Jason

26

使用numpy.corrcoef函数代替numpy.correlate函数来计算延迟t的统计相关性：

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

- Ramón J Romero y Vigil

“相关系数”不是指统计学中使用的自相关，而是信号处理中使用的自相关吗？https://en.wikipedia.org/wiki/Autocorrelation#Signal_processing - Daniel Node.js

1

@DanielPendergast 我对信号处理不太熟悉。从NumPy文档中可以看到：“返回Pearson积矩相关系数”。这是否是信号处理版本？ - Ramón J Romero y Vigil

12

你的问题1已经在这里的一些优秀回答中得到了广泛讨论。

我想与您分享几行代码，可以基于自相关的数学性质计算信号的自相关。也就是说，可以按以下方式计算自相关：

从信号中减去均值并获得无偏信号
计算无偏信号的傅里叶变换
通过对无偏信号的傅里叶变换的每个值求平方范数来计算信号的功率谱密度
计算功率谱密度的反傅里叶变换
通过无偏信号的平方和来规范化功率谱密度的反傅里叶变换，并仅取结果向量的一半

实现此操作的代码如下：

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)

- Ruggero

@pindakaas，你能具体说明一下吗？请提供有关哪些函数存在差异的信息。 - Ruggero

@dylnan 是的，它本来可以工作，但是它会不必要地计算许多平方根，然后再返回到平方模块。 - Ruggero

1

是的，但你有没有意识到使用列表推导式可能会更慢。 - Jason

你的结果（大致上）等同于 numpy.correlate(x, x, mode='same')，但是完整函数 numpy.correlate(x, x, mode='full') 呢？我们能否也使用FFT来实现它呢？ - Feodoran

顺便提一下：np.abs(f)**2 比列表推导式要快得多。原则上，通过使用 f.real**2 + f.imag**2 避免不必要的平方甚至更快，但是在这里效果微不足道，因为FFT当然比一些基本的乘法/加法慢。 - Feodoran

显示剩余4条评论

12

由于我刚刚遇到同样的问题，我想与你分享几行代码。实际上，stackoverflow 上已经有几篇关于自相关性的相似帖子。如果您将自相关定义为 a(x, L) = sum(k=0,N-L-1)((xk-xbar)*(x(k+L)-xbar))/sum(k=0,N-1)((xk-xbar)**2) [这是IDL的a_correlate函数中给出的定义，并且它与问题#12269834的答案2一致]，那么以下似乎会给出正确的结果：

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

你可以看到，我已经用正弦曲线和均匀随机分布进行了测试，两个结果都符合我的预期。请注意，我使用了mode="same"而不是其他人使用的mode="full"。

- maschu

10

在statsmodels.tsa.stattools.acf()中提供了numpy.correlate的替代方法。它产生一个连续递减的自相关函数，就像OP所描述的那样。实现它相当简单：

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]

默认行为是在40个滞后值处停止，但可以使用nlag=选项来调整您特定应用程序的滞后值。该函数背后的统计数据有一个引用，请参考页面底部的文献。

- fisherp

4

我是一名计算生物学家，在计算随机过程时间序列之间的自相关/互相关时，我发现np.correlate不能满足我的需求。

事实上，np.correlate似乎缺少距离为d的所有可能时间点对之间的平均值。

以下是我定义实现所需功能的函数的方法：

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

似乎之前的回答都没有涵盖到自相关/互相关的这种情况：希望这个回答对像我一样从事随机过程研究的人有用。

- Orso

对于那些想要使用这个答案的人，最后一行可以通过NumPy加速：N = len(x); v = c/(N - np.abs(np.arange(N) - N/2))。 - jared

2

使用傅里叶变换和卷积定理

时间复杂度为N*log(N)

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

这是一个标准化且无偏的版本，它也是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]

A. Levy提供的方法可行，但我在我的电脑上测试过，它的时间复杂度似乎是N*N

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

- Jian

2

没有使用pandas的简单解决方案:

import numpy as np

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

- dignitas

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- A. Levy · Accepted Answer

回答你的第一个问题，numpy.correlate(a, v, mode)执行的是将a与v的反向卷积，并将结果剪裁为指定模式。卷积的定义C(t)=∑_{-∞ < i < ∞}a_iv_t+i，其中-∞ < t < ∞，允许结果从-∞到∞，但显然无法存储无限长的数组。因此必须剪裁，这就是模式发挥作用的地方。有3种不同的模式：full、same和valid：

"full"模式返回每个t都有一些重叠的a和v的结果。
"same"模式返回与最短向量（a或v）长度相同的结果。
"valid"模式仅在a和v完全重叠时返回结果。关于模式的更多细节，请参见numpy.convolve的文档。

对于您的第二个问题，我认为numpy.correlate实际上正在给您提供自相关性，它只是给您更多。自相关性用于查找信号或函数在某个时间差异处与其本身的相似程度。在时间差异为0时，自相关性应该最高，因为信号与其本身完全相同，因此您预期自相关结果数组中的第一个元素将是最大的。但是，相关性并没有从时间差为0开始。它从负时间差开始，接近0，然后变为正数。也就是说，您希望得到的是：

自相关(a) = ∑ _{-∞ < i < ∞} a_iv_t+i 其中 0 <= t < ∞

但您得到的是：

自相关(a) = ∑ _{-∞ < i < ∞} a_iv_t+i 其中 -∞ < t < ∞

您需要做的是取相关结果的后一半，那就是您要找的自相关性。一个简单的Python函数可以这样做：

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

当然，你需要进行错误检查以确保x实际上是一个一维数组。此外，这个解释可能不是最数学严谨的。我一直在使用无限大，因为卷积的定义使用它们，但这并不一定适用于自相关。因此，这个解释的理论部分可能有些奇怪，但希望实际结果有所帮助。这些页面关于自相关非常有帮助，如果你不介意穿过符号和重概念，可以给你更好的理论背景。