我如何使用 numpy.correlate 来做自相关？

为了回答你的第一个问题，numpy.correlate(a, v, mode)正在执行 a与 v相反的卷积，并给出由指定模式剪切的结果。其中-∞ < t < ∞的 卷积的定义卷积的定义，C (t) = ∑ _{- ∞ < i < ∞} a_我v_{T + i}允许从-∞到∞的结果，但显然不能存储无限长的数组。所以它必须被剪裁，这就是模式出现的原因。有三种不同的模式: 完整、相同和有效:

• “ full”模式返回每个 t的结果，其中 a和 v都有一些重叠。
• “ same”模式返回与最短向量(a或 v)具有相同长度的结果。
• “有效”模式只有在 a和 v完全重叠时才返回结果。文件 for numpy.convolve提供了关于模式的更多细节。

对于你的第二个问题，我认为 numpy.correlate 是给你的自相关，它只是给你一点点以及。自相关是用来找出一个信号，或函数，是多么相似的本身在一定的时间差。在时差为0时，自相关应该是最高的，因为信号与它本身是相同的，所以您期望自相关结果阵列中的第一个元素将是最大的。但是，相关性不是从0的时差开始的。它从负的时差开始，接近0，然后变为正的。也就是说，你期望:

自相关(a) = ∑ _{- ∞ < i < ∞} a_我v_{T + i}，其中0 < = t < ∞

但你得到的是:

自相关(a) = ∑ _{- ∞ < i < ∞} a_我v_{T + i}，其中-∞ < t < ∞

您需要做的是取出相关性结果的后半部分，这应该就是您正在寻找的自相关性。一个简单的 python 函数可以这样做:

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

当然，您需要进行错误检查，以确保 x实际上是一个1-d 数组。此外，这种解释在数学上可能不是最严谨的。我一直在讨论无穷大，因为卷积的定义使用了它们，但这并不一定适用于自相关。因此，这个解释的理论部分可能有点站不住脚，但希望实际结果是有帮助的。关于自相关的 这些 页数是非常有帮助的，如果你不介意涉水通过符号和繁重的概念，可以给你一个更好的理论背景。

自相关有两个版本: 统计和卷积。它们都执行相同的操作，除了一些细节: 统计版本被标准化为在区间[-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)])

由于我刚刚遇到了同样的问题，我想与您分享几行代码。事实上，现在已经有几篇相当类似的关于堆栈溢出中自相关性的文章。如果你将自相关定义为 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 _ correlevant 函数中给出的定义，它与我在问题 # 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"。

使用 numpy.corrcoef函数而不是 numpy.correlate来计算 t 滞后的统计相关性:

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

我认为 OP 问题的真正答案简洁地包含在 Numpy 相关文档的这段摘录中:

mode : {'valid', 'same', 'full'}, optional
Refer to the `convolve` docstring.  Note that the default
is `valid`, unlike `convolve`, which uses `full`.

这意味着，在没有“ mode”定义的情况下使用时，当为其两个输入参数(即-用于执行自相关性时)提供相同的向量时，Numpy.correlate 函数将返回一个标量。

我使用 talib.CORREL 来实现这样的自相关，我怀疑你也可以对其他包做同样的事情:

def autocorrelate(x, period):


# x is a deep indicator array
# period of sample and slices of comparison


# oldest data (period of input array) may be nan; remove it
x = x[-np.count_nonzero(~np.isnan(x)):]
# subtract mean to normalize indicator
x -= np.mean(x)
# isolate the recent sample to be autocorrelated
sample = x[-period:]
# create slices of indicator data
correls = []
for n in range((len(x)-1), period, -1):
alpha = period + n
slices = (x[-alpha:])[:period]
# compare each slice to the recent sample
correls.append(ta.CORREL(slices, sample, period)[-1])
# fill in zeros for sample overlap period of recent correlations
for n in range(period,0,-1):
correls.append(0)
# oldest data (autocorrelation period) will be nan; remove it
correls = np.array(correls[-np.count_nonzero(~np.isnan(correls)):])


return correls


# CORRELATION OF BEST FIT
# the highest value correlation
max_value = np.max(correls)
# index of the best correlation
max_index = np.argmax(correls)

你的问题1已经在这里的几个出色的回答中得到了广泛的讨论。

我想和你们分享几行代码，这些代码允许你们仅仅基于自相关的数学性质来计算信号的自相关性。也就是说，自相关可以通过以下方式计算:

1. 从信号中减去平均值，得到一个无偏的信号

2. 计算无偏信号的傅里叶变换

3. 计算信号的功率谱密度，方法是取无偏信号傅里叶变换的每个值的平方范数

4. 计算功率傅里叶变换的反谱密度

5. 用无偏信号的平方和来标准化功率傅里叶变换的反向谱密度，然后只取得所得矢量的一半

执行此操作的代码如下:

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)

我是一个计算生物学家，当我必须计算随机过程的时间序列之间的自动/交叉相关时，我意识到 np.correlate没有做我需要的工作。

事实上，np.correlate似乎缺少的是距离 τ 的 对所有可能的时间点进行平均。

下面是我如何定义一个函数来完成我需要的任务:

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

在我看来，以前的答案似乎都没有涵盖这个自动/互相关的实例，希望这个答案可能对像我这样从事随机过程的人有用。

我认为有两件事情会使这个话题更加混乱:

1. 统计与信号处理定义: 正如其他人指出的，在统计学中，我们将自相关归一化为[ -1,1]。
2. 偏和非偏均值/方差: 当时间序列在滞后 > 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, no padding, non-partial(红线)的结果是不同的，因为它在做 FFT 之前没有用0填充时间序列，所以它是循环 FFT。我不能详细解释为什么，这是我从其他地方学到的。

一个没有熊猫的简单解决方案:

import numpy as np


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

用傅里叶变换和卷积定理

时间复杂度为 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]

利维提供的方法是可行的，但我在自己的电脑上测试了一下，它的时间复杂度似乎是 N * N

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

绘制给定熊猫数据的统计自相关性图表返回值序列:

import matplotlib.pyplot as plt


def plot_autocorr(returns, lags):
autocorrelation = []
for lag in range(lags+1):
corr_lag = returns.corr(returns.shift(-lag))
autocorrelation.append(corr_lag)
plt.plot(range(lags+1), autocorrelation, '--o')
plt.xticks(range(lags+1))
return np.array(autocorrelation)

在 Statmodels.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 nlags，但这可以为您的特定应用程序的 nlag=选项进行调整。在页面的底部有一个关于 函数背后的统计信息的引用。

要使用 numpy 实现 IDL 的 相关人员函数，需要使用 mode="full"运行 np.correlate并将 n-1添加到 lag数组。

def a_correlate(y, lag):
y = np.asarray(y)
lag = np.asarray(lag)
n = len(y)
yunbiased = y - np.mean(y)
ynorm = np.sum(yunbiased**2)
r = np.correlate(yunbiased, yunbiased, "full") / ynorm
return r[lag + (n - 1)]

示例(基于上面链接的 IDL 文档页面中的示例) :

# Define an n-element sample population:
X = np.array([3.73, 3.67, 3.77, 3.83, 4.67, 5.87, 6.70, 6.97, 6.40, 5.57])
# Compute the autocorrelation of X for LAG = -3, 0, 1, 3, 4, 8:
lag = [-3, 0, 1, 3, 4, 8]
result = a_correlate(X, lag)
print(result)
# prints: [ 0.01461851  1.          0.81087925  0.01461851 -0.32527914 -0.15168379]