nlinfit.f90

非线性最小二乘法,L-M算法

PROGRAM D9R8A
!Driver for routine MRQMIN
! EXTERNAL FGAUSS
PARAMETER(NPT=100,MA=3,SPREAD=0.001)
DIMENSION X(NPT),Y(NPT),SIG(NPT),A(MA),LISTA(MA),&
          COVAR(MA,MA),ALPHA(MA,MA),GUES(MA)
DATA A/5.0,2.0,3.0/
DATA GUES/4.5,2.2,2.8/
IDUM=-911
!First try a sum of two Gaussians
DO I=1,100
    X(I)=0.1*I
    Y(I)=0.0
!    DO J=1,4,3
        J=1
        Y(I)=Y(I)+A(J)*EXP(-((X(I)-A(J+1))/A(J+2))**2)
!    END DO
    Y(I)=Y(I)*(1.0+SPREAD*GASDEV(IDUM))
    SIG(I)=SPREAD*Y(I)
END DO
MFIT=MA
DO I=1,MFIT
    LISTA(I)=I
END DO
ALAMDA=-1
DO I=1,MA
    A(I)=GUES(I)
END DO
CALL MRQMIN(X,Y,SIG,NPT,A,MA,LISTA,MFIT,COVAR,ALPHA,&
                 MA,CHISQ,FGAUSS,ALAMDA)
K=1
ITST=0
DO
    WRITE(*,'(/1X,A,I2,T18,A,F10.4,T43,A,E9.2)')&
    'Iteration #',K,'Chi-squared:',CHISQ,'ALAMDA:',ALAMDA
    WRITE(*,'(1X,T5,A,T13,A,T21,A,T29,A,T37,A,T45,A)')&
            'A(1)','A(2)','A(3)','A(4)','A(5)','A(6)'
    WRITE(*,'(1X,6F8.4)') (A(I),I=1,6)
    K=K+1
    OCHISQ=CHISQ
    CALL MRQMIN(X,Y,SIG,NPT,A,MA,LISTA,MFIT,COVAR,ALPHA,&
                 MA,CHISQ,FGAUSS,ALAMDA)
    IF (CHISQ>OCHISQ) THEN 
        ITST=0
    ELSE IF (ABS(OCHISQ-CHISQ)<0.1) THEN
        ITST=ITST+1
    ENDIF
    IF (ITST>=2)  EXIT
END DO
ALAMDA=0.0
CALL MRQMIN(X,Y,SIG,NPT,A,MA,LISTA,MFIT,COVAR,ALPHA,&
                 MA,CHISQ,FGAUSS,ALAMDA)
WRITE(*,*) 'Uncertainties:'
WRITE(*,'(1X,6F8.4/)') (SQRT(COVAR(I,I)),I=1,6)
WRITE(*,'(1X,A)') 'Expected results:'
WRITE(*,'(1X,F7.2,5F8.2/)') 5.0,2.0,3.0,2.0,5.0,3.0
END




!!!!!!!!!!!    GASDEV     !!!!!!!!!!!!!!!!!!!!!!!
FUNCTION gasdev(idum)
INTEGER idum
REAL gasdev
! USES ran1
INTEGER iset
REAL fac,gset,rsq,v1,v2,ran1
SAVE iset,gset
DATA iset/0/
if (iset==0) then
  do
    v1=2.*ran1(idum)-1.
    v2=2.*ran1(idum)-1.
    rsq=v1**2+v2**2
    if(.NOT.(rsq>=1..or.rsq==0.)) EXIT
  end do
  fac=sqrt(-2.*log(rsq)/rsq)
  gset=v1*fac
  gasdev=v2*fac
  iset=1
else
  gasdev=gset
  iset=0
endif
END FUNCTION gasdev

!!!!!!!!!!!!!!!!!!!!!!!!!!! 
FUNCTION ran1(idum)
INTEGER idum,ia,im,iq,ir,ntab,ndiv
REAL ran1,am,eps,rnmx
PARAMETER (ia=16807,im=2147483647,am=1./im,&
     iq=127773,ir=2836,ntab=32,ndiv=1+(im-1)/ntab,&
	 eps=1.2e-7,rnmx=1.-eps)
INTEGER j,k,iv(ntab),iy
SAVE iv,iy
DATA iv /ntab*0/, iy /0/
if (idum<=0.or.iy==0) then
  idum=max(-idum,1)
  do j=ntab+8,1,-1
    k=idum/iq
    idum=ia*(idum-k*iq)-ir*k
    if (idum<0) idum=idum+im
    if (j<=ntab) iv(j)=idum
  end do
  iy=iv(1)
endif
k=idum/iq
idum=ia*(idum-k*iq)-ir*k
if (idum<0) idum=idum+im
j=1+iy/ndiv
iy=iv(j)
iv(j)=idum
ran1=min(am*iy,rnmx)
END FUNCTION ran1
!!!!!!!!!!!!!!!!!!!!!!!!!!!


!!!!!!!!!!!!!!!!MRQMIN!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE mrqmin(X,Y,SIG,NDATA,A,MA,LISTA,MFIT,&
           COVAR,ALPHA,NCA,CHISQ,FUNCS,ALAMDA)
INTEGER ma,nca,ndata,lista(ma),MMAX
REAL alamda,chisq,funcs,a(ma),alpha(nca,nca),&
     covar(nca,nca),sig(ndata),x(ndata),y(ndata)
PARAMETER (MMAX=20)
!USES covsrt,gaussj,mrqcof
INTEGER j,k,l,mfit
REAL ochisq,atry(MMAX),beta(MMAX),da(MMAX)
if(alamda<0.) then
  kk=mfit+1
  do j=1,ma
    ihit=0
    do k=1,mfit
      if(lista(k)==j) ihit=ihit+1
    end do
    if(ihit==0) then
      lista(kk)=j
      kk=kk+1
    else if(ihit>1) then
      pause 'improper permutation in lista'
    endif
  end do
  if(kk/=(ma+1)) pause 'improper permutation in lista'
  alamda=0.001
  call mrqcof(x,y,sig,ndata,a,ma,lista,mfit,alpha,&
                     beta,nca,chisq,funcs)
  ochisq=chisq
  do j=1,ma
    atry(j)=a(j)
  end do
endif
do j=1,mfit
  do k=1,mfit
    covar(j,k)=alpha(j,k)
  end do
  covar(j,j)=alpha(j,j)*(1.+alamda)
  da(j)=beta(j)
end do
call gaussj(covar,mfit,da)
if(alamda==0.) then
  call covsrt(covar,nca,ma,lista,mfit)
  return
endif
do j=1,mfit
  atry(lista(j))=a(lista(j))+da(j)
end do
call mrqcof(x,y,sig,ndata,atry,ma,lista,mfit,covar,da,&
                nca,chisq,funcs)
if(chisq<ochisq) then
  alamda=0.1*alamda
  ochisq=chisq
  do j=1,mfit
    do k=1,mfit
      alpha(j,k)=covar(j,k)
    end do
    beta(j)=da(j)
    a(lista(j))=atry(lista(j))
  end do
else
  alamda=10.*alamda
  chisq=ochisq
endif
END SUBROUTINE mrqmin

!!!!!!!!!!!!!!!!!MRQCOF!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE mrqcof(x,y,sig,ndata,a,ma,lista,mfit,&
             alpha,beta,nalp,chisq,funcs)
PARAMETER (mmax=20)
DIMENSION x(ndata),y(ndata),alpha(nalp,nalp),&
   beta(ma),dyda(mmax),lista(mfit),sig(ndata),a(ma)
do j=1,mfit
  do k=1,j
    alpha(j,k)=0.
  end do
  beta(j)=0.
end do
chisq=0.
do i=1,ndata
  call funcs(x(i),a,ymod,dyda,ma)
  sig2i=1./(sig(i)*sig(i))
  dy=y(i)-ymod
  do j=1,mfit
    wt=dyda(lista(j))*sig2i
    do k=1,j
      alpha(j,k)=alpha(j,k)+wt*dyda(lista(k))
    end do
    beta(j)=beta(j)+dy*wt
  end do
  chisq=chisq+dy*dy*sig2i
end do
do j=2,mfit
  do k=1,j-1
    alpha(k,j)=alpha(j,k)
  end do
end do
END SUBROUTINE mrqcof

!!!!!!!!!!!!!  funcs  !!!!!!!!!!!!!!!!!!!
SUBROUTINE funcs(x,a,ymod,dyda,ma)
INTEGER ma
REAL x,ymod,a(ma),dyda(ma)
INTEGER i
!REAL arg,ex,fac
ymod=A(1)*EXP(-((X(I)-A(2))/A(3))**2)
do 	i=1,ma-2
  dyda(i)=EXP(-((X(I)-A(2))/A(3))**2)
  dyda(i+1)=EXP(-((X(I)-A(2))/A(3))**2)*2.*(X(I)-A(2))/A(3)
  dyda(i+2)=EXP(-((X(I)-A(2))/A(3))**2)*2.*(X(I)-a(2))**2./a(3)**3
enddo
end SUBROUTINE funcs




!!!!!!!!!!!!!  COVSRT  !!!!!!!!!!!!!!!!!!!

SUBROUTINE covsrt(covar,ncvm,ma,lista,mfit)
INTEGER ma,mfit,ncvm,lista(mfit)
REAL covar(ncvm,ncvm)
INTEGER i,j,k
REAL swap
do j=1,ma-1
  do i=j+1,ma
    covar(i,j)=0.
  end do
end do
do i=1,mfit-1
  do j=i+1,mfit
    if(lista(j)>lista(i)) then
      covar(lista(j),lista(i))=covar(i,j)
    else
      covar(lista(i),lista(j))=covar(i,j)
    endif
  end do
end do
swap=covar(1,1)
do j=1,ma
  covar(1,j)=covar(j,j)
  covar(j,j)=0.
end do
covar(lista(1),lista(1))=swap
do j=2,mfit
  covar(lista(j),lista(j))=covar(1,j)
end do
do j=2,ma
  do i=1,j-1
    covar(i,j)=covar(j,i)
  end do
end do
END SUBROUTINE covsrt
!!!!!!!!!!!!!  GAUSSJ  !!!!!!!!!!!!!!!!!!!!!!!!!!

SUBROUTINE gaussj(a,n,b)
INTEGER n,NMAX
REAL a(n,n),b(n)
PARAMETER (NMAX=50)
INTEGER i,icol,irow,j,k,l,ll,indxc(NMAX),&
        indxr(NMAX),ipiv(NMAX)
REAL big,dum,pivinv
do j=1,n
  ipiv(j)=0
end do
do i=1,n
  big=0.
  do j=1,n
    if(ipiv(j)/=1) then
      do k=1,n
        if (ipiv(k)==0) then
          if (abs(a(j,k))>=big) then
            big=abs(a(j,k))
            irow=j
            icol=k
          endif
        else if (ipiv(k)>1) then
          pause 'singular matrix in gaussj'
        endif
      end do
    endif
  end do
  ipiv(icol)=ipiv(icol)+1
  if (irow/=icol) then
    do l=1,n
      dum=a(irow,l)
      a(irow,l)=a(icol,l)
      a(icol,l)=dum
    end do
    dum=b(irow)
	b(irow)=b(icol)
	b(icol)=dum
  endif
  indxr(i)=irow
  indxc(i)=icol
  if (a(icol,icol)==0.) pause 'singular matrix in gaussj'
  pivinv=1./a(icol,icol)
  a(icol,icol)=1.
  do l=1,n
    a(icol,l)=a(icol,l)*pivinv
  end do
  b(icol)=b(icol)*pivinv
  do ll=1,n
    if(ll/=icol) then
      dum=a(ll,icol)
      a(ll,icol)=0.
      do l=1,n
        a(ll,l)=a(ll,l)-a(icol,l)*dum
      end do
      b(ll)=b(ll)-b(icol)*dum
    endif
  end do
end do
do l=n,1,-1
  if(indxr(l)/=indxc(l)) then
    do k=1,n
      dum=a(k,indxr(l))
      a(k,indxr(l))=a(k,indxc(l))
      a(k,indxc(l))=dum
    end do
  endif
end do
END SUBROUTINE gaussj