我需要对一组数字进行自相关，据我所知，这只是该组与自身的相关。
我已经尝试使用 numpy 的相关函数，但我不相信结果，因为它几乎总是给出第一个数字不是最大的向量，因为它应该是。
所以，这个问题实际上是两个问题：
- 究竟在
numpy.correlate做什么？ - 我如何使用它（或其他东西）进行自相关？
Ben
问问题
224010 次
13 回答
134
要回答您的第一个问题，numpy.correlate(a, v, mode)是执行a与反向的卷积v并给出由指定模式裁剪的结果。卷积的定义，C(t)=∑ -∞ < i < ∞ a i v t+i where -∞ < t < ∞，允许从 -∞ 到 ∞ 的结果，但你显然不能存储无限长的大批。所以它必须被剪裁，这就是模式的用武之地。有 3 种不同的模式：完整、相同和有效：
- “完整”模式返回每个
t位置的结果，a并且v有一些重叠。 a“相同”模式返回与最短向量 (或)长度相同的结果v。- “有效”模式仅在完全重叠
a时返回结果。v的文档提供numpy.convolve了有关模式的更多详细信息。
对于你的第二个问题，我认为numpy.correlate 是给你自相关，它也只是给你更多。自相关用于找出信号或函数在某个时间差与其自身的相似程度。在时间差为 0 时，自相关应该是最高的，因为信号与其自身相同，因此您预计自相关结果数组中的第一个元素将是最大的。但是，相关性并不是从 0 的时间差开始。它从负的时间差开始，接近 0，然后变为正数。也就是说，您期望：
自相关 (a) = ∑ -∞ < i < ∞ a i v t+i其中 0 <= t < ∞
但你得到的是：
自相关(a) = ∑ -∞ < i < ∞ a i v t+i其中 -∞ < t < ∞</p>
您需要做的是获取相关结果的后半部分，这应该是您正在寻找的自相关。一个简单的python函数可以做到这一点：
def autocorr(x):
result = numpy.correlate(x, x, mode='full')
return result[result.size/2:]
当然，您需要进行错误检查以确保它x实际上是一维数组。此外，这种解释可能不是最严格的数学解释。我一直在抛出无穷大，因为卷积的定义使用它们，但这并不一定适用于自相关。因此，这种解释的理论部分可能有点不可靠，但希望实际结果会有所帮助。这些关于自相关的 页面非常有帮助，如果您不介意涉足符号和繁重的概念，它们可以为您提供更好的理论背景。
于 2009-03-24T06:09:02.647 回答
33
自相关有两个版本：统计和卷积。它们都做同样的事情，除了一点细节：统计版本被标准化为区间 [-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)])
于 2011-11-02T13:27:37.917 回答
27
我认为有两件事给这个话题增加了混乱：
- 统计与信号处理的定义：正如其他人指出的那样，在统计中，我们将自相关归一化为 [-1,1]。
- 部分与非部分均值/方差：当时间序列以滞后 > 0 移动时，它们的重叠大小将始终 < 原始长度。我们是使用原始（非部分）的均值和标准差，还是总是使用不断变化的重叠（部分）计算新的均值和标准差。（这可能有一个正式的术语，但我现在要使用“部分”）。
我创建了 5 个函数来计算一维数组的自相关，具有部分与非部分的区别。有些使用统计公式，有些使用信号处理意义上的相关，也可以通过 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。我无法详细解释为什么，这是我从其他地方学到的。
于 2018-07-04T07:32:40.810 回答
24
使用numpy.corrcoef函数而不是numpy.correlate计算滞后 t 的统计相关性：
def autocorr(x, t=1):
return numpy.corrcoef(numpy.array([x[:-t], x[t:]]))
于 2014-02-18T19:10:37.863 回答
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()
如您所见，我已经使用 sin 曲线和均匀随机分布对此进行了测试，这两个结果看起来都符合我的预期。请注意，我使用mode="same"而不是mode="full"像其他人那样使用。
于 2013-06-13T14:52:10.863 回答
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)
于 2016-10-20T12:46:02.253 回答
8
numpy.correlate 的替代方法在statsmodels.tsa.stattools.acf()中可用。这会产生一个不断减小的自相关函数，就像 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 处停止，但这可以通过nlag=针对您的特定应用程序的选项进行调整。该功能背后的统计数据在页面底部有一个引文。
于 2019-04-24T18:38:24.547 回答
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
在我看来，以前的答案都没有涵盖这个自/互相关的例子：希望这个答案对像我这样从事随机过程的人有用。
于 2018-04-12T15:25:30.047 回答
1
我使用 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)
于 2015-12-22T05:06:48.107 回答
1
没有熊猫的简单解决方案：
import numpy as np
def auto_corrcoef(x):
return np.corrcoef(x[1:-1], x[2:])[0,1]
于 2018-07-23T13:22:38.553 回答
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]
A. Levy 提供的方法有效，但我在我的 PC 上测试过，它的时间复杂度似乎是N*N
def autocorr(x):
result = numpy.correlate(x, x, mode='full')
return result[result.size/2:]
于 2018-09-17T06:20:55.287 回答
0
我认为 OP 问题的真正答案简洁地包含在 Numpy.correlate 文档的摘录中：
mode : {'valid', 'same', 'full'}, optional
Refer to the `convolve` docstring. Note that the default
is `valid`, unlike `convolve`, which uses `full`.
这意味着，当没有“模式”定义使用时，Numpy.correlate 函数将返回一个标量，当给定其两个输入参数的相同向量时（即 - 当用于执行自相关时）。
于 2015-02-01T21:17:42.833 回答
0
绘制给定 pandas datatime 返回系列的统计自相关：
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)
于 2018-10-08T12:42:32.750 回答