回归系数计算开始对，但是后来加上求滤去线性时候，输出的是0kb

dzxconan

        parameter (nx=180,ny=89,m=55,undef=32767)
    real sst(nx,ny,m),in(m),rxy(nx,ny),sx(nx,ny),sy,avex(nx,ny),avey,b(nx,ny),f(nx,ny)
        real a(nx,ny,m),rr(nx,ny),t(nx,ny),sxa(nx,ny),sya,avexa(nx,ny),aveya
!        real sst1(nx,ny,n1),in1(n1)
       
!-----input------
open(1,file='ssts.grd',form='binary')
        do k=1,m
        read(1)((sst(i,j,k),i=1,nx),j=1,ny)
        enddo
        open(2,file='index.grd',form='binary')
        read(2)(in(i),i=1,m)

!------计算系数和F检验-----

call xiangguan(sst,in,rxy,sx,sy,avex,avey)
do j=1,ny
do i=1,nx
if(sst(i,j,1)/=undef)then
b(i,j)=rxy(i,j)*sqrt(sx(i,j))/sqrt(sy)
f(i,j)=rxy(i,j)**2/(1-rxy(i,j)**2)*(m-2)
else
b(i,j)=undef
f(i,j)=undef
endif
enddo
enddo

do k=1,m
do j=1,ny
do i=1,nx
if(sst(i,j,k)/=undef)then
a(i,j,k)=sst(i,j,k)-b(i,j)*avey
else
a(i,j,k)=undef
endif
enddo
enddo
enddo


open(6,file='d:\paper\160\asummer.txt')

write(6,*)((f(i,j),i=1,nx),j=i,ny)

!call xiangguan(a,in,rr,sxa,sya,avexa,aveya)
!-------------t检验---------------------------


!------output-------------------
open(3,file='d:\paper\160\bsummer.grd',form='binary')
open(4,file='d:\paper\160\fsummer.grd',form='binary')
open(5,file='d:\paper\160\atsummer.grd',form='binary')
!open(6,file='d:\paper\160\asummer.grd',form='binary')
write(3)((b(i,j),i=1,nx),j=1,ny)
write(4)((f(i,j),i=1,nx),j=1,ny)
write(5)((b(i,j),i=1,nx),j=i,ny)
!write(5)((t(i,j),i=1,nx),j=i,ny)本来想写输出a，发现是0kb后来改成b或者f还是0kb好神奇啊~~~~~~~~~~~
!do k=1,m
!write(6)((a(i,j,k),i=1,nx),j=i,ny)
!enddo
end

!---------子程序：计算相关系数-----
subroutine xiangguan(a,b,r,sa,sb,avea,aveb)
        parameter (nx=180,ny=89,m=55,undef=32767)
real a(nx,ny,m),b(m),sa(nx,ny),sb,sab(nx,ny),r(nx,ny),avea(nx,ny),aveb

do k=1,m
do j=1,ny
do i=1,nx
if(a(i,j,1)/=undef)then
avea(i,j)=avea(i,j)+a(i,j,k)/m
endif
enddo
enddo
aveb=aveb+b(k)/m
enddo

do k=1,m
do j=1,ny
do i=1,nx
if(a(i,j,1)/=undef)then
sa(i,j)=sa(i,j)+(a(i,j,k)-avea(i,j))**2/m
sab(i,j)=sab(i,j)+(a(i,j,k)-avea(i,j))*(b(k)-aveb)/m
endif
enddo
enddo
sb=sb+(b(k)-aveb)**2/m
enddo

do j=1,ny
do i=1,nx
if(a(i,j,1)/=undef)then
r(i,j)=sab(i,j)/sqrt(sa(i,j))/sqrt(sb)
endif
!if(abs(r(i,j))>=1)then
!write(*,*)r(i,j)
!endif
enddo
enddo


end

言深深

程序太长，不容易看出来的···建议你逐步检测数据吧，看看哪一步出错了···

dzxconan

言深深 发表于 2013-4-29 17:46

                               
登录/注册后可看大图

程序太长，不容易看出来的···建议你逐步检测数据吧，看看哪一步出错了···

print是对的 可是写不进文件是怎么回事啊？？

天籁之音

程序太长，不容易看出来的···建议你逐步检测数据吧，看看哪一步出错了···多谢

言深深

dzxconan 发表于 2013-4-29 21:07

                               
登录/注册后可看大图

print是对的 可是写不进文件是怎么回事啊？？

你写到txt里面看看是否出错···类似于6的写法，给一个输出格式···如果数据不出错的话再write称binary格式的

dzxconan

言深深 发表于 2013-4-30 09:16

                               
登录/注册后可看大图

你写到txt里面看看是否出错···类似于6的写法，给一个输出格式···如果数据不出错的话再write称binar ...

open(3,file='d:\paper\160\bsummer.grd',form='binary')
open(4,file='d:\paper\160\fsummer.grd',form='binary')
open(50,file='d:\paper\160\arsummer.txt')
!open(6,file='d:\paper\160\asummer.grd',form='binary')
write(3)((b(i,j),i=1,nx),j=1,ny)
write(4)((f(i,j),i=1,nx),j=1,ny)
print*,rr
write(50,*)((rr (i,j),i=1,nx),j=i,ny)、
我这样写输出 可是print没有问题 txt只有1k 里面什么都没有

dzxconan

天籁之音 发表于 2013-4-29 22:57

                               
登录/注册后可看大图

程序太长，不容易看出来的···建议你逐步检测数据吧，看看哪一步出错了···多谢

但是后来同时输出 为什么开始没有问题的数组后来就有问题了呢？？

言深深

dzxconan 发表于 2013-4-30 13:41

                               
登录/注册后可看大图

open(3,file='d:\paper\160\bsummer.grd',form='binary')
open(4,file='d:\paper\160\fsummer.grd',form ...

用循环输出
do j=i,ny
do i=1,nx
print*,rr (i,j)
enddo
enddo

风儿飘飘

怎么不对回归系数进行显著性检验呢

风儿飘飘

学习到不少知识，谢谢