计算有效自由度，俩公式，蒙特卡洛，MATLAB

伽蓝鸟

经常被问到有效自由度的事情，正好写个帖子说下这事儿。算相关或者回归的时候，需要用到t检验。这里有个假定，数据是独立同分布的。
独立同分布这个假定，在年际数据上常常能得以满足。今年的降水跟去年的降水，一般来说是没啥关系的（TBO暂且不讨论）。所以自由度按照样本数-2来算就成了。
但是如果对数据做了低通滤波，那么独立性就没法得以满足了。举个栗子：

这张图里，左边是原始数据，基本可以认为是独立的；右边是低通滤波之后的数据，有很强的的自相关（滞相关），今年的数据跟去年的数据差不了太多。
肉眼来看的话，右边的图里只有6~7个波，也就是说大概只有六七个独立样本。这时候算这俩序列的相关系数，非常大。不能用原始序列的自由度（62-2=60）来检验，应该考虑有效自由度。
有效自由度的估算公式有很多，例如我之前用过这个：

舍友上次推荐了另一个公式，是这个：

公式是咋推导出来的就别问我了。。我的统计没学到那么好。。
根据这俩公式，算出有效自由度，再根据有效自由度查到相关系数临界值，就能知道相关是否显著了。


我之前一直用公式（1）嘛，它比较简单，只考虑了一阶滞相关。最近看到公式（2）之后，想看看哪个公式更靠谱。
于是就折腾了俩蒙特卡洛模拟，来对比下结果。
思路很简单，数据一共62年，我扔随机数，从62年里有放回地取62年出来。
对重新采样的数据，算个相关系数——这个是跟年际尺度的相关做对比的。
对重新采样的数据，低通滤波，再算相关——这个跟年代际的相关来做对比。
重复很多次（比如一万次），就能得到相关系数的分布了。

重采样之后直接算相关，5%和95%分位数是-0.2090和0.2104，而自由度为60时alpha=0.1的相关系数临界值是0.211，对应得很好。
重采样，滤波，再算相关，5%和95%分位数是-0.447和0.456，更接近于第二个公式算出来的相关系数临界值。所以考虑了更高阶滞相关的公式在这个例子里可能是更靠谱的。

当然，实际应用的时候，有两种做法。
一个是根据公式，算出有效自由度，查到相关系数临界值，跟算到的相关系数对比。
还可以直接用蒙特卡洛模拟。拿算到的相关系数，跟蒙特卡洛得到的相关系数临界值，对比下，就知道了。可以不必反着去算有效自由度的。



测试用的数据和代码（俩公式、蒙特卡洛）见附件。

pearal

之前自己算过,检验两个高度自相关的序列是否相关,第二个公式计算结果确实更靠谱.

vivianal

给力！！！！

croton

很有用！

JGD飘

感谢楼主分享！没接触过Matlab，可以问一下function "autocorr"自相关的函数需要自己写吗还是哪里可以下载的

JGD飘

提示“Undefined function or variable 'autocorr'.”

伽蓝鸟

JGD飘 发表于 2021-12-14 21:41
提示“Undefined function or variable 'autocorr'.”

这是内置函数啊，咋可能没有呢

JGD飘

可能是我在虚拟机上运行的原因吧，可能相关函数包没有，我没找着autocorr的源码，所以我用了xcorr，依据网上的公式，变换了一下

爽歪歪

请问高通滤波后，需要重新计算有效自由度吗？谢谢！

伽蓝鸟

爽歪歪 发表于 2021-12-18 17:11
请问高通滤波后，需要重新计算有效自由度吗？谢谢！

也需要啊，但是高通滤波的时候，r1可能小于0，自由度看起来会增加。实际上还是取序列数-2作为自由度，不会让它额外增加的。