我需要做的自相关的一组数字,因为我明白这仅仅是相关的设置本身。

我已经试过了使用顽固的相互关联的功能,但是我不相信的结果,因为它几乎总是给一个向量在第一个号码是 最大的,因为它应该是。

因此,这个问题是真正的两个问题:

  1. 到底是什么 numpy.correlate 在做什么?
  2. 我如何可以使用(或者别的东西)以做到自动相关的?
有帮助吗?

解决方案

来回答你的第一个问题, numpy.correlate(a, v, mode) 执行该卷的 a 与反的 v 和给予的结果裁剪过指定的模式。的 定义的卷积, C(t)=∑ -∞ < 我 < ∞ 一个vt+i 在那里-∞ < t < ∞,允许对结果从-∞到∞,但你明显地不能储存的一个无限长的阵列。因此,它已经被剪辑,这就是模式。有3种不同的模式:全,同样,&有效的:

  • "全面"方式返回的结果为每一个 t 其中两个 av 有一些重叠。
  • "相同的"模式返回的结果一样长短的矢量(av).
  • "有效的"方式返回的结果只有当 av 完全相互重叠。的 文档 对于 numpy.convolve 提供更详细的模式。

为你的第二个问题,我认为 numpy.correlate 给你自相关,它就是给你多一点为好。自相关的是用来寻找如何类似的一个信号,或功能,是本身在一定的时间差。在一个时间差异为0,自动相关的应当是最高的,因为的信号完全相同本身,所以你期望的第一个元素自相关结果列将是最大的。然而,相关性并不是起始时间的差异为0。它开始在一个负面的时间差,关闭为0,然后进积极的。就是,你期待:

自相关(a)=∑ -∞ < 我 < ∞ 一个vt+i 其中0 <=t < ∞

但你得到了什么是:

自相关(a)=∑ -∞ < 我 < ∞ 一个vt+i 在那里-∞ < t < ∞

你需要做的是采取最后一半的相关结果,而这应该是自相关你正在寻找。一个简单的蟒蛇的功能做到这将是:

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

使用 numpy.corrcoef 功能来代替 numpy.correlate 计算出的统计相关的滞后的t:

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

我只是碰到了同样的问题,我想和大家分享一些代码行与你同在。事实上,有几个相当类似的员额有关的自相关,在计算器了。如果你定义自相关的作为 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功能和它同意我所看到的在答案2的问题 #12269834],则以下似乎得到正确结果:

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" 正如其他人所做的。

我想有2个东西,增加混乱到这个话题:

  1. 统计v.s.信号处理的定义:正如其他人已经指出的,在统计,我们恢复正常自动的相关成[-1,1].
  2. 部分诉s.非部分的意思是/差额:当时间序列的转变,在一个滞后>0,他们的重叠大小会总是 < 原来的长度。我们使用的平均和标准的原始(非部分),或者总是计算新的平均和标准采用的不断变化的重叠(部分)有差别。(有可能是一个正式的术语,但我会使用"部分"为现在)。

我已经创建了5职能,计算自相关的一个1d阵列,部分诉s.非部分的区别。一些使用公式的统计数据,某些使用相关信号处理意义上,这也可以通过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()

这里是输出的数字:

enter image description here

我们没有看到所有5行,因为3人重叠(在紫).重叠是所有非局部自动的相关性。这是因为计算从信号处理方法(np.correlate, FFT)不计算不同的意思是/std对每一重叠。

还注意到, fft, no padding, non-partial (红线)的结果是不同的,因为它没有衬垫的时间序列与0s前做FFT,所以它是圆形的FFT。我不能详细解释为什么,那是什么我学到从其他地方。

你的问题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)

我使用塔利.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)

使用傅立叶变换和卷的定理

时间complexicity是 N*日志(N)

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

这里是一个规范化和公正的版本,它也是 N*日志(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.征税工作,但是我的测试,它在我的电脑,其时间complexicity似乎是 N*N

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

我认为真正的答案运的问题是简明扼要地载于本摘录的顽固.相关文件:

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

这意味着,在使用时没有'模式'的定义,顽固.相关的功能将返回的一个标,当给予同样的矢量,为其两个输入参数(即-当用来执行自相关).

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

import numpy as np

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

绘制统计的自相关性给予了大熊猫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)

我是一个计算生物学家,并且当我必须计算的自动/交的相关性夫妻之间的时间系列的随机过程中我意识到 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

在我看来没有任何先前的答复盖这种实例的自动/交相关:希望这个答案可能是有益的人工作上的随机处理等等我。

替代顽固.相关的提供 statsmodels.tsa。"菜单".acf().这将产生一个不断下降的自相关功能,像一个说明通过任择议定书》。实现这是相当简单:

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]

默认的行为是停止在40nlags,但这可以调整 nlag= 选择你的特定应用程序。有一个引文在页底部的 统计数据背后的功能.

许可以下: CC-BY-SA归因
不隶属于 StackOverflow
scroll top