zbhg.f90

逐步回归分析的通用程序(Fortran95版)

!  *********************************************************************
!  逐步回归程序
!  功能:求出逐步回归方程,其中a(m)为优选出来的回归方程系数,a0为优选的回归方程常数项
!  部分重要变量说明:
!  m----自变量个数,qn----原始数据行数,qq(qn)----已知观测数据的输出值(因变量)
!  xy(qn,m+1)----增广矩阵,前M列相应于m个变元的输入值,m+1列相应于因变量qq的输出值
!  s----最终回归方程产生的误差
!  ********************************************************************
program main
  implicit none
  integer,parameter :: m=4,qn=365
  integer i,j,k,i0,k0
  real a0,pmax,pmin,s
  real,parameter ::  f1=0.1,f2=0.1
  integer nmax,nmin
  real :: qq(qn),xy(qn,m+1),a(m),q(qn),xp(m+1),&
                    cgma(m+1),rxy(m+1,m+1),rx(m+1,m+1),qql(qn),qc(qn)
  open(1,file='data2.txt')
  do i=1,qn
  read(1,*)xy(i,1),xy(i,2),xy(i,3),xy(i,4),xy(i,5)
  end do
  close(1)
  qq=xy(:,m+1)
  call zupu(m,qn,f1,f2,qq,xy,a0,a,q,s,xp,cgma,rxy,rx,qql,qc)
  open(3,file='jieguo.txt')
  write(3,*)a0,a,s
  close(3)
  write(*,*)a0,a,s


contains
subroutine zupu(m,qn,f1,f2,qq,xy,a0,a,q,s,xp,cgma,rxy,rx,qql,qc)
implicit none
  integer m,qn
  integer i,j,k,i0,k0,k1,i1,i2,i3,i4
  real a0,pmax,pmin,s
  real f1,f2
  real r,f
  integer nmax,nmin
  real :: qq(qn),xy(qn,m+1),a(m),q(qn),xp(m+1),&
                    cgma(m+1),rxy(m+1,m+1),rx(m+1,m+1),qql(qn),qc(qn)
  do i=1,qn
     qql(i)=qq(i)
  end do
  do j=1,m+1
     xp(j)=0.0
     do i=1,qn
        xp(j)=xp(j)+xy(i,j)
     end do
	 xp(j)=xp(j)/qn
  end do
  do i=1,m+1
  do j=1,m+1
     rxy(i,j)=0.0
     do k=1,qn
        rxy(i,j)=rxy(i,j)+(xy(k,i)-xp(i))*(xy(k,j)-xp(j))
     end do
     if(i==j)cgma(i)=sqrt(rxy(i,j))
  end do
  end do
  do i=1,m+1
  do j=1,m+1
     rxy(i,j)=rxy(i,j)/(cgma(i)*cgma(j))
 end do
  end do

  k0=qn-1
  
do i1=1,m
     a(i1)=0.0
 end do
 do i=1,m
     a0=0.0
     pmax=0.0
     nmax=0
     nmin=0
     pmin=1e9
     do i0=1,m
        if(rxy(i0,i0)>1e-2)then
          r=(rxy(i0,m+1)*rxy(m+1,i0))/rxy(i0,i0)
        if(r==0)exit
        if(r<0)then
          a(i0)=cgma(m+1)*rxy(i0,m+1)/cgma(i0)
          if(abs(r)<abs(pmin))then
          pmin=r
          nmin=i0
        end if
        else
        if(r>0)then
          if(r>pmax)then
          pmax=r
          nmax=i0
          end if
        end if
        end if
        end if
     end do
     a0=xp(m+1)
       do i2=1,m
          a0=a0-a(i2)*xp(i2)
	   end do
     f=abs(pmin)*k0/rxy(m+1,m+1)
       if(f<=f1)then
         k1=nmin
         k0=k0+1
       else
     f=pmax*(k0-1)/(rxy(m+1,m+1)-pmax)
       if(f<=(f1+f2))then
        exit
      else
        k1=nmax
        k0=k0-1
      end if
      end if
      do i3=1,m+1
      do j=1,m+1
         rx(i3,j)=rxy(i3,j)
      end do
      end do
      do i4=1,m+1
      do j=1,m+1
         if(i4==k1 .and. j==k1)then
           rxy(i4,j)=1/rx(k1,k1)
         else if(i4==k1)then
           rxy(i4,j)=rx(k1,j)/rx(k1,k1)
         else if(j==k1)then
           rxy(i4,j)=-rx(i4,k1)/rx(k1,k1)
         else
           rxy(i4,j)=rx(i4,j)-rx(i4,k1)*rx(k1,j)/rx(k1,k1)
         end if
      end do
      end do
   end do
  
    s=0
         do k=1,qn
            q(k)=a0
            do j=1,m
               q(k)=a(j)*xy(k,j)+q(k)
            end do
            qc(k)=(qql(k)-q(k))/qql(k)*100
            s=s+qc(k)*qc(k)
        end do
   
   return
   end subroutine zupu
end program main
!*******************************************************