！！！一元线性回归，系数a画出来太平滑，不知道哪里出问题了！！！

星雨

!******************************************************************!
!                       本程序用来求解一元线性回归                 !  
!******************************************************************!
Program main   
Parameter my=33,mo=12,nx=144,ny=72
Real x(mo,my),y(nx,ny,mo,my),avex(4,my),avey(nx,ny,4,my)
        integer i,j,it,im
        !!!!!!!!!!it=year,im=month!!!!!!!!!!
        real a(nx,ny,mo),b(nx,ny,mo),r(nx,ny,mo),qx(mo),xx(mo),yy(nx,ny,mo),dx(nx,ny,mo),dy(nx,ny,mo),dxy(nx,ny,mo),t(nx,ny,mo)
        !!!!!!!!!xx=ave_x,yy=ave_y,dx,dy,dxy=zhongjianliang!!!!!!!!!!!!!
    open(11,file='c:\2013hg\precip.grd',form='binary',access='direct',recl=nx*ny)
    Open(22,file='c:\2013hg\aoi.grd',form='binary')
Print*,"开始读入数据："
     Do j=1,my
  do i=1,mo
            Read(22),x(i,j)
     End do
  end do
  print*,'aoi is ok'
     print*,x(1,1)
  print*,'precip'
irec=1
       do  it=1,my
    do  im=1,mo
       read(11,rec=irec)((y(i,j,im,it),i=1,nx),j=1,ny)
       irec=irec+1
      enddo
   enddo
print*,'read chang-var is ok!'
     Close(11)
     Close(22)
     Print*,"数据读入完毕！"
   do j=1,my  
  avex(1,j)=(x(3,j)+x(4,j)+x(5,j))/3
  avex(2,j)=(x(6,j)+x(7,j)+x(8,j))/3
  avex(3,j)=(x(9,j)+x(10,j)+x(11,j))/3
     avex(4,j)=(x(12,j)+x(1,j)+x(2,j))/3
   enddo

do j=1,ny
do i=1,nx
do it=1,my
     avey(i,j,1,it)=(y(i,j,3,it)+y(i,j,4,it)+y(i,j,5,it))/3
     avey(i,j,2,it)=(y(i,j,6,it)+y(i,j,7,it)+y(i,j,8,it))/3
     avey(i,j,3,it)=(y(i,j,9,it)+y(i,j,10,it)+y(i,j,11,it))/3
     avey(i,j,4,it)=(y(i,j,12,it)+y(i,j,1,it)+y(i,j,2,it))/3
enddo
enddo
enddo
!????????????????????????????????????????????????????????????????????

print*,'sub'
do im=1,4
do it=1,my
xx(im)=xx(im)+avex(im,it)
enddo
enddo

do j=1,ny
do i=1,nx
do im=1,4
do it=1,my
      yy(i,j,im)=yy(i,j,im)+avey(i,j,im,it)
      dy(i,j,im)=dy(i,j,im)+avey(i,j,im,it)*avey(i,j,im,it)
enddo
enddo
enddo
enddo
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

do j=1,ny
do i=1,nx
do im=1,4
   dy(i,j,im)=dy(i,j,im)-yy(i,j,im)*yy(i,j,im)/my
enddo
enddo
enddo

do im=1,4
   xx(im)=xx(im)/my
enddo
print*,'XXs'
print*, xx(2)
do j=1,ny
do i=1,nx
do im=1,4
yy(i,j,im)=yy(i,j,im)/my
enddo
enddo
enddo
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
do j=1,ny
do i=1,nx
do im=1,4
do it=1,my
qx(im)=avex(im,it)-xx(im)
dx(i,j,im)=dx(i,j,im)+qx(im)*qx(im)
dxy(i,j,im)=dxy(i,j,im)+qx(im)*(avey(i,j,im,it)-yy(i,j,im))
enddo
enddo
enddo
enddo
Print*,'begin'
do im=1,4
do j=1,ny
do i=1,nx
   a(i,j,im)=dxy(i,j,im)/dx(i,j,im)
   b(i,j,im)=yy(i,j,im)-a(i,j,im)*xx(im)
   r(i,j,im)=sqrt(dxy(i,j,im)*dxy(i,j,im)/(dx(i,j,im)*dy(i,j,im)))
enddo
enddo
enddo
print*,'a,b,r has finished !'
print*,'begin to write the data !'
do im=1,4
do j=1,ny
do i=1,nx
if(r(i,j,2)<=-1) then
print*,'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
endif
enddo
enddo
enddo
open(100,file='c:/2013hg/sj/pa4.grd',form='binary')
open(200,file='c:/2013hg/sj/pa4.txt',form='formatted')
open(300,file='c:/2013hg/sj/pr4.grd',form='binary')
open(400,file='c:/2013hg/sj/pr4.txt',form='formatted')
do im=1,4
do j=1,ny
do i=1,nx
write(100) a(i,j,im)
write(300) r(i,j,im)
write(200,1000) a(i,j,im)
write(400,1000) r(i,j,im)
enddo
enddo
enddo
1000 format(f8.2)

close(100)
close(200)
close(300)
close(400)
   Print *,'对以上相关系数进行显著性检验!'
do im=1,4
do j=1,ny
do i=1,nx
   t(i,j,im)=r(i,j,im)*sqrt((my-2)/(1-r(i,j,im)**2))
enddo
enddo
enddo
open(3333,file='c:/2013hg/sj/t4.grd',form='binary')
open(4444,file='c:/2013hg/sj/t4.txt',form='formatted')
do im=1,4
do j=1,ny
do i=1,nx
write(3333) t(i,j,im)
write(4444,1000) t(i,j,im)
enddo
enddo
enddo
close(3333)
close(4444)
End program



不知道哪位高手能给我RUN一下，看看哪里出问题了？？？？？？？？附图

river

你的相关系数绝对值那么大，相关性未免太好了吧。计算过程中都没有赋初值啊，除了求距平，剩下的求平均、均方差、协方差都要赋初值的啊

kongfeng0824

你的这个程序的图，是直接用fortran调用grads话的吗？还是输出数据之后才画的呢？

星雨

river 发表于 2013-5-3 22:15

                               
登录/注册后可看大图

你的相关系数绝对值那么大，相关性未免太好了吧。计算过程中都没有赋初值啊，除了求距平，剩下的求平均、均 ...

这是回归指数，找人看了看，已经找出问题了，赋初值问题谢谢你提醒,好人一生平安！

星雨

kongfeng0824 发表于 2013-5-4 02:56

                               
登录/注册后可看大图

你的这个程序的图，是直接用fortran调用grads话的吗？还是输出数据之后才画的呢？

输出后才画的，有错，谢谢关心！

river

星雨 发表于 2013-5-7 20:40

                               
登录/注册后可看大图

这是回归指数，找人看了看，已经找出问题了，赋初值问题谢谢你提醒,好人一生平安！

刚开始编程的时候都这样。但是你这都做论文快毕业了，开始的稍微有点儿晚哈，不过时间足够了

星雨

river 发表于 2013-5-7 20:44

                               
登录/注册后可看大图

刚开始编程的时候都这样。但是你这都做论文快毕业了，开始的稍微有点儿晚哈，不过时间足够了

还以为自己学的很好，一做论文发现啥也不会，啥也不精。。。没学好啊。。。

river

星雨 发表于 2013-5-7 21:25

                               
登录/注册后可看大图

还以为自己学的很好，一做论文发现啥也不会，啥也不精。。。没学好啊。。。

加油吧骚年，一个论文仔细做下来，四年的很多东西就掌握了

风儿飘飘

正在做回归分析，学习了

atlaszhang

灌水是为了下载。。。。