t_test_cha.f90

t检验的fortran程序 看是否对大家有帮助

!---------------- t检验 -------------------------－－--!
!气候稳定性检验，均值检验法之一：                      !
!  X1(n1)、X1(n2):需作检验的2个序列;                   !
!  输出t，确定显著性水平a=0.05等等,自由度v＝n1+n2-2    !
!  查t分布表ta＝？，t>ta则认为在此显著性水平下样本     !
!     X1(n1)、X1(n2)有显著差异                         !
!  a=0.1, v=11, ta＝1.80                               !
!  a=0.05, v=11, ta＝2.20                              !
!  a=0.01, v=11, ta＝3.11                              !
!                                     李晓峰           !
!                                   2005年4月1日       !
!------------------------------------------------------!

parameter (n1=8,n2=5,mm=144,nn=73,nz=1)
real xu(mm,nn,nz,n1),yu(mm,nn,nz,n2),tu(mm,nn,nz),x1(n1),x2(n2)
real xv(mm,nn,nz,n1),yv(mm,nn,nz,n2),tv(mm,nn,nz)
real xh(mm,nn,nz,n1),yh(mm,nn,nz,n2),th(mm,nn,nz)
real tuv(mm,nn,nz)
character*50 infile1,infile2,outfile
infile1='i:\t-test\h1.grd'
infile2='i:\t-test\h2.grd'
outfile='i:\t-test\T.grd'

open(10,file=infile1,form='binary')
do it=1,n1
read(10)(((xu(i,j,iz,it),i=1,mm),j=1,nn),iz=1,nz)
read(10)(((xv(i,j,iz,it),i=1,mm),j=1,nn),iz=1,nz)
read(10)(((xh(i,j,iz,it),i=1,mm),j=1,nn),iz=1,nz)
enddo
close(10)

open(11,file=infile2,form='binary')
do it=1,n2
read(11)(((yu(i,j,iz,it),i=1,mm),j=1,nn),iz=1,nz)
read(11)(((yv(i,j,iz,it),i=1,mm),j=1,nn),iz=1,nz)
read(11)(((yh(i,j,iz,it),i=1,mm),j=1,nn),iz=1,nz)
enddo
close(11)

write(*,*)'ok'

do i=1,mm
  do j=1,nn
    do iz=1,nz

      do it=1,n1
 	     x1(it)=xu(i,j,iz,it)
	  enddo

	  do it=1,n2
	     x2(it)=yu(i,j,iz,it)
      enddo

      call T_test_cha(n1,n2,x1,x2,t)
    
	  tu(i,j,iz)=t
	enddo
  enddo
enddo

do i=1,mm
  do j=1,nn
    do iz=1,nz

      do it=1,n1
 	     x1(it)=xv(i,j,iz,it)
	  enddo

	  do it=1,n2
	     x2(it)=yv(i,j,iz,it)
      enddo

      call T_test_cha(n1,n2,x1,x2,t)
    
	  tv(i,j,iz)=t
	enddo
  enddo
enddo

do i=1,mm
  do j=1,nn
    do iz=1,nz

      do it=1,n1
 	     x1(it)=xh(i,j,iz,it)
	  enddo

	  do it=1,n2
	     x2(it)=yh(i,j,iz,it)
      enddo

      call T_test_cha(n1,n2,x1,x2,t)
    
	  th(i,j,iz)=t
	enddo
  enddo
enddo

do i=1,mm
  do j=1,nn
    do iz=1,nz
      tuv(i,j,iz)=sqrt(tu(i,j,iz)**2+tv(i,j,iz)**2)
	enddo
  enddo
enddo

open(12,file=outfile,form='binary')
write(12)(((tuv(i,j,iz),i=1,mm),j=1,nn),iz=1,nz)
write(12)(((th(i,j,iz),i=1,mm),j=1,nn),iz=1,nz)
close(12)

write(*,*)(((tu(i,j,iz),i=1,mm),j=1,nn),iz=1,nz)
end

!***************************************************************88

!----------------------------------------
subroutine T_test_cha(n1,n2,X1,X2,t)
implicit none
integer,intent(in)::n1,n2
real,dimension(n1),intent(in)::X1,X2
real,intent(out)::t
real::ave1,ave2,var1,var2,Sp

call avevar(n1,X1,ave1,var1)
call avevar(n2,X2,ave2,var2)
Sp=((N1-1)*VAR1+(N2-1)*VAR2)/(N1+N2-2)
Sp=SQRT(Sp)
t=(ave2-ave1)/(Sp*(sqrt(1.0/n1+1.0/n2)))

end

!--------------------------------
subroutine AveVar(n,X,ave,var)
implicit none
integer,intent(in)::n
real,dimension(n)::X(n)
real,intent(out)::ave,var
integer::i

ave=0.0
do i=1,n
   ave=ave+X(i)/float(n)                          !作用相当于 ave=0.0
enddo                                             !           do i=1,n;ave=ave+X(i)/float(n);enddo
var=0.0                                           !           ave=ave/float(n) 
do i=1,n
   var=var+(X(i)-ave)*(X(i)-ave)/float(n)
enddo
end