matlab std函数_关于numpy互相关函数np.correlate的一点疑问

956d21b6cdfb7e46e5d01b1268eedf85.png

最近再做一个小作业,是关于寻找太平洋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

似乎在信号处理中,定义的确是这样。但是类似功能的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/scipy​github.com
1dc527113705f56e93975cb495ba4182.png

上述问题中,唯一可以满足我的需求的竟然是matplotlib.pyplot中的plt.xcorr(x, y, normed=True, detrend=<function detrend_none at 0x7f1ed58c08c0>, usevlines=True, maxlags=10, *, data=None, **kwargs)

  • 0
    点赞
  • 2
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值