一元回归分析的代码问题

chieniou

需要做云南（1960年到2011年52个春季）平均温度（做过标准化处理）与necp全球500hPa  omega（也是春季平均，且做过距平处理）的回归分析并做显著性检验
以下是求回归系数b的代码
real varx(52),sumvarx,avevarx,dvarx(52),vary(144,73,52),sumvary(144,73),
avevary(144,73),dvary(144,73,52),sumdvary(144,73),re(144,73)
open(1,file='e:\730\dat\omega_d_500_spring.dat',form='binary',status='old')
open(2,file='e:\730\dat\tema_spring_standardization.dat',form='binary',status='old')
open(3,file='e:\730\dat\re_omega_d_500_spring.dat',form='binary',status='replace')
read(1)(((vary(x,y,t),x=1,144),y=1,73),t=1,52)
read(2)(varx(t),t=1,52)

!以下是求预报因子(云南春季平均温度)距平平方和。（因为坐过标准化，因此距平等于本身）
sumvarx=0
do t=1,52
  sumvarx=sumdvarx+varx(t)**2
enddo

!以下是对预报量距平与预报因子距平之积求和
sumvary(1,1)=0
do y=1,73
  do x=1,144
    do t=1,52
       sumvary(x,y)=sumvary(x,y)+vary(x,y,t)
    enddo
  enddo
enddo

sumdvary(1,1)=0
do y=1,73
  do x=1,144
    do t=1,52
       avevary(x,y)=sumvary(x,y)/52
       dvary(x,y,t)=vary(x,y,t)-avevary(x,y)
       sumdvary(x,y)=sumdvary(x,y)+dvary(x,y,t)*dvarx(t)
    enddo
  enddo
enddo

!以下是求回归系数
do y=1,73
  do x=1,144
     write(3)sumdvary(x,y)/sumdvarx
  enddo
enddo
close(1)
close(2)
close(3)
end

以下是做t检验的代码（感觉就是这部分不对）
real varx(52),sumvarx,c,vary(144,73,52),sumvary(144,73),avevary(144,73),
dvary(144,73,52),Q(144,73),re(144,73),tt(144,73),estimator(144,73,52)
open(1,file='e:\730\dat\omega_d_300_spring.dat',form='binary',status='old')
open(2,file='e:\730\dat\tema_spring_standardization.dat',form='binary',status='old')
open(3,file='e:\730\dat\re_omega_d_300_spring.dat',form='binary',status='old')
open(4,file='e:\730\dat\test_omega_d_300_spring.dat',form='binary',status='replace')
open(5,file='e:\730\dat\omega_d_300_spring_estimator.dat',form='binary',status='replace')
read(1)(((vary(x,y,t),x=1,144),y=1,73),t=1,52)
read(2)(varx(t),t=1,52)
read(3)((re(x,y),x=1,144),y=1,73)
!print*,vary(34,54,33)
!print*,varx(44)
!print*,re(55,33)

!以下是求预报因子距平平方和
c=0
do t=1,52
  c=c+varx(t)**2
enddo
!以下是求残差平方和Q
do x=1,144
  do y=1,73
    do t=1,52
      dvary(x,y,t)=vary(x,y,t)-re(x,y)*varx(t)
    enddo
  enddo
enddo
Q(1,1)=0
do x=1,144
  do y=1,73
    do t=1,52
      Q(x,y)=Q(x,y)+dvary(x,y,t)**2
    enddo
  enddo
enddo

!以下是做t检验，写出估计量vary
do t=1,52
  do y=1,73
    do x=1,144
      tt(x,y)=sqrt(50*(re(x,y)**2)*c/Q(x,y))
      write(4)tt(x,y)
      estimator(x,y,t)=re(x,y)*varx(t)
      write(5)estimator(x,y,t)
    enddo
  enddo
enddo

close(1)
close(2)
close(3)
close(4)
close(5)
end

我做出来的结果是，相关显著的地方预报量vary与估计量vary相差很远，这不对吧？
我自己检查了几次程序逻辑，“没发现”问题，
所以请各位大神帮看程序有没有逻辑错误，或者说算法错误什么的

chieniou

zhouling 发表于 2015-5-3 15:13
貌似你写的t-检验公式有问题，楼主可参考http://bbs.06climate.com/forum.php?mod=viewthread&tid=18528&ex ...

我之前弄错了   我写的是对回归系数是否为0的F检验  只是统计量用的是t遵从自由度为（n-2）的t分布

chieniou

有个地方打错了，就是求回归系数那，sumvarx少打了个d  应该是sumdvarx，还有定义变量的地方也是掉了个d  不过这都不是问题所在

zhouling

貌似你写的t-检验公式有问题，楼主可参考http://bbs.06climate.com/forum.p ... p;extra=&page=1