请选择 进入手机版 | 继续访问电脑版
爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 20540|回复: 28

[源程序] 计算有效自由度,俩公式,蒙特卡洛,MATLAB

[复制链接]

新浪微博达人勋

发表于 2021-12-13 10:48:48 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x
本帖最后由 伽蓝鸟 于 2021-12-13 10:52 编辑

经常被问到有效自由度的事情,正好写个帖子说下这事儿。算相关或者回归的时候,需要用到t检验。这里有个假定,数据是独立同分布的。
独立同分布这个假定,在年际数据上常常能得以满足。今年的降水跟去年的降水,一般来说是没啥关系的(TBO暂且不讨论)。所以自由度按照样本数-2来算就成了。
但是如果对数据做了低通滤波,那么独立性就没法得以满足了。举个栗子:
1.png
这张图里,左边是原始数据,基本可以认为是独立的;右边是低通滤波之后的数据,有很强的的自相关(滞相关),今年的数据跟去年的数据差不了太多。
肉眼来看的话,右边的图里只有6~7个波,也就是说大概只有六七个独立样本。这时候算这俩序列的相关系数,非常大。不能用原始序列的自由度(62-2=60)来检验,应该考虑有效自由度。
有效自由度的估算公式有很多,例如我之前用过这个:
2.png
舍友上次推荐了另一个公式,是这个:
3.png
公式是咋推导出来的就别问我了。。我的统计没学到那么好。。
根据这俩公式,算出有效自由度,再根据有效自由度查到相关系数临界值,就能知道相关是否显著了。


我之前一直用公式(1)嘛,它比较简单,只考虑了一阶滞相关。最近看到公式(2)之后,想看看哪个公式更靠谱。
于是就折腾了俩蒙特卡洛模拟,来对比下结果。
思路很简单,数据一共62年,我扔随机数,从62年里有放回地取62年出来。
对重新采样的数据,算个相关系数——这个是跟年际尺度的相关做对比的。
对重新采样的数据,低通滤波,再算相关——这个跟年代际的相关来做对比。
重复很多次(比如一万次),就能得到相关系数的分布了。
4.png 5.png
重采样之后直接算相关,5%和95%分位数是-0.2090和0.2104,而自由度为60时alpha=0.1的相关系数临界值是0.211,对应得很好。
重采样,滤波,再算相关,5%和95%分位数是-0.447和0.456,更接近于第二个公式算出来的相关系数临界值。所以考虑了更高阶滞相关的公式在这个例子里可能是更靠谱的。

当然,实际应用的时候,有两种做法。

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



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

compare_Ndof.m

1.68 KB, 下载次数: 161, 下载积分: 金钱 -5

sampledata.mat

1.24 KB, 下载次数: 100, 下载积分: 金钱 -5

评分

参与人数 2金钱 +11 贡献 +1 收起 理由
huangbicheng + 6 + 1 很给力!
vivianal + 5 很给力!

查看全部评分

密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2021-12-13 11:35:31 | 显示全部楼层
之前自己算过,检验两个高度自相关的序列是否相关,第二个公式计算结果确实更靠谱.
密码修改失败请联系微信:mofangbao
回复 支持 1 反对 0

使用道具 举报

新浪微博达人勋

发表于 2021-12-13 10:56:55 | 显示全部楼层
给力!!!!
密码修改失败请联系微信:mofangbao
回复

使用道具 举报

新浪微博达人勋

发表于 2021-12-14 11:50:44 | 显示全部楼层
很有用!
密码修改失败请联系微信:mofangbao
回复

使用道具 举报

新浪微博达人勋

发表于 2021-12-14 21:05:49 | 显示全部楼层
感谢楼主分享!没接触过Matlab,可以问一下function "autocorr"自相关的函数需要自己写吗还是哪里可以下载的
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2021-12-14 21:41:30 | 显示全部楼层
提示“Undefined function or variable 'autocorr'.”
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2021-12-15 08:58:30 | 显示全部楼层
JGD飘 发表于 2021-12-14 21:41
提示“Undefined function or variable 'autocorr'.”

这是内置函数啊,咋可能没有呢
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2021-12-15 09:56:37 | 显示全部楼层
可能是我在虚拟机上运行的原因吧,可能相关函数包没有,我没找着autocorr的源码,所以我用了xcorr,依据网上的公式,变换了一下
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2021-12-18 17:11:24 | 显示全部楼层
请问高通滤波后,需要重新计算有效自由度吗?谢谢!
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2021-12-18 19:23:00 | 显示全部楼层
爽歪歪 发表于 2021-12-18 17:11
请问高通滤波后,需要重新计算有效自由度吗?谢谢!

也需要啊,但是高通滤波的时候,r1可能小于0,自由度看起来会增加。实际上还是取序列数-2作为自由度,不会让它额外增加的。
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表