线性回归去除enso等信号的影响（程序分享）

happyaliyun

线性回归拟合去除ENSO信号：指数为自变量，对其做一元线性回归拟合，拟合方程的残差即表示统计上与ENSO 无关。

变量与ENSO做回归，再用原来的变量减去回归后的Y，

得到去ENSO趋势的变量：Y-（ax+b）

a=Cov(nino,z*)/Var(nino),协方差/方差

b=Ymean-aXmean

下面是ncl程序去除enso信号，对于去除iod等信号是一样的方法

;去除enso

  ;ts求回归方程的系数a,b。

cov=escovc(nino,ts({LAT|latS:latN},{LON|lonL:lonR},time|:));lat*lon

  var=variance(nino)                                         ;一个常数

  a=cov/var                    ;a=协方差/方差               ;lat*lon

  b=dim_avg_n(ts,0)-a*avg(nino);b=y(平均)-a*x(平均);注意此处ts平均仅是时间维度的平均;lat*lon

ts_guji=new((/ntime,nlat,nlon/),"float")

  do i=0,ntime-1

  ts_guji(i,:,:)=nino(i)*a(:,:)+b ;y=a*x+b

  end do

  ts_free=ts-ts_guji

  copy_VarMeta(ts,ts_free)

  printVarSummary(ts_free)

  printMinMax(ts_free, True)

叩容

楼主这个方法，结合楼上提到的escorc函数，可以写为：

1. a= escorc_n(sst_DJF_ano(:,{latS:latN},{lonL:lonR}), ENSO, 0, 0)
   
2. copy_VarMeta(sst_DJF_ano(0,:,:),a)
   
3. b=dim_avg_n(sst_DJF_ano,0)-a*avg(ENSO)
   
4. ;b=y(平均)-a*x(平均);注意此处ts平均仅是时间维度的平均;lat*lon
   
5. dims = dimsizes(sst_DJF_ano)
   
6. ntime = dims(0)
   
7. ;nlat = dims(1)
   
8. ;nlon = dims(2)
   
9. ;sst_guji=new((/ntime,nlat,nlon/),"float")
   
10. sst_guji=sst_DJF_ano
    
11. do i=0,ntime-1
    
12.         sst_guji(i,:,:)=ENSO(i)*a(:,:)+b ;y=a*x+b
    
13. end do
    
14. ;  printVarSummary(sst_guji)
    
15. sst_free=sst_DJF_ano-sst_guji
    
16. copy_VarMeta(sst_DJF_ano,sst_free)

复制代码

不过，上面这个方法，去除ENSO信号会去不干净。建议使用回归函数来做：

1. a= regCoef_n(ENSO,sst_DJF_ano,0,0)
   
2. copy_VarMeta(sst_DJF_ano(0,:,:),a)
   
3. sst_guji=sst_DJF_ano
   
4. dims = dimsizes(sst_DJF_ano)
   
5. ntime = dims(0)
   
6. do i=0,ntime-1
   
7.         sst_guji(i,:,:)=ENSO(i)*a(:,:)
   
8. end do
   
9. ;  printVarSummary(sst_guji)
   
10. sst_free=sst_DJF_ano-sst_guji
    
11. copy_VarMeta(sst_DJF_ano,sst_free)

复制代码

zrongli

我试了一下用你的方法算的a和escorc函数求得不同，发现你少计算了ts的方差var2
应该是a=cov/sqrt(var1*var2)

白衣风月倾天下

哇 感谢分享

zrongli

感谢分享

人可山

感谢分享

帅蛋fly

叩容 发表于 2023-8-11 20:29
楼主这个方法，结合楼上提到的escorc函数，可以写为：

天哪这个楼主太宝藏了！！！！给你点赞