![956d21b6cdfb7e46e5d01b1268eedf85.png](https://i-blog.csdnimg.cn/blog_migrate/ad4654b15f9d3c66fa5529e0c9d43f36.jpeg)
最近再做一个小作业,是关于寻找太平洋Nino3.4区和热带印度洋(TI)海温(SST)的最大超前滞后关系,不可避免需要求解互相关函数。其中numpy,scipy,statsmodels均有求解互相关的函数,最后选择numpy.correlate。但是我发现其中关于互相关的定义好像有点问题(或者函数其实并不完善),其中问题如下:
先来看一下numpy.correlate计算定义:c_{av}[k] = sum_n a[n+k] * conj(v[n]) (其中conj为求解共轭,对于实数时间序列,不必考虑),简单来说numpy.correlate做的事情就是计算出两个序列之间的"错位点积"
a
那么当时滞为1(即a落后b一位)时,"错位点积"为2×2+3×4=16,当时滞为-1(即a超前b一位),“错位点积”为1×4+2×5=14,其他情况也是以此类推。
那么问题就来了,函数中并没提供'unbiased'和'normalize'的选项,仅仅是求解了"错位点积",这样就可以代表互相关函数吗?
在我查阅资料过程中,我翻到了知乎上的这篇文章:
明天:卷积、互相关与自相关zhuanlan.zhihu.com![742d4991db6ec2e3f14783a91beda7e4.png](https://i-blog.csdnimg.cn/blog_migrate/f18c7dff0d201c928d8b0861cbcec412.png)
似乎在信号处理中,定义的确是这样。但是类似功能的MATLAB函数xcorr提供了'unbiased'和'normalize'的选项,同时还有'maxlags'规定最大时滞的参数。statsmodels.tsa.stattools.ccf中有提供'unbiased'选项。
那么我在numpy.correlate中至少存在三点疑问:
①是否应该在函数中加入'unbiased'选项,具体的无偏化如下:
a=np.array([1,2,3,4,5])
b=np.array([2,3,4,5,6])
a_d=a-a.mean()
b_d=b-b.mean()
p_unbiased=np.correlate(a_d,b_d,'full')/ (len(a)-np.arange(-4,5,1))
即协方差分母为n-k,n为a的长度,k为时滞的绝对值,这一部分是由statsmodels.tsa.stattools.ccf源码为我提供的思路。
②是否应该加入'normalize'选项,具体归一化过程如下:
a=np.array([1,2,3,4,5])
b=np.array([2,3,4,5,6])
a_d=a-a.mean()
b_d=b-b.mean()
p_unbiased=np.correlate(a_d,b_d,'full')/ (len(a)-np.arange(-4,5,1))
p_norm=p_unbiased/(np.std(a_d)*np.std(b_d))
这正好是MATLAB函数xcorr使用的归一化方法。
③是否应该加入maxlags选项,有时候我们并不需要太长互相关函数值序列,只需要某一最大时滞的互相关函数序列,比如逐月SST值,一般取maxlags=12即可,再长就没有实际意义了。同时加入maxlags选项后,用户可以很方便的截取,否则需要自行处理。这一部分其实在GitHub中四年前就有提出疑问。
Add maxlag kw for limited cross correlation window size to correlate/convolve functions · Issue #4940 · scipy/scipygithub.com![1dc527113705f56e93975cb495ba4182.png](https://i-blog.csdnimg.cn/blog_migrate/974ee8a65c74db44aea77cbdad1f13bb.jpeg)
上述问题中,唯一可以满足我的需求的竟然是matplotlib.pyplot中的plt.xcorr(x, y, normed=True, detrend=<function detrend_none at 0x7f1ed58c08c0>, usevlines=True, maxlags=10, *, data=None, **kwargs)