N.I.M.R.O.D.  

marquardt.f90

Go to the documentation of this file.
00001 
00002 !------------------------------------------------------------------------------
00003 ! N.I.M.R.O.D. - Normal approximation Inference in Models with Random
00004 !                effects based on Ordinary Differential equations
00005 !------------------------------------------------------------------------------
00006 !
00007 ! VERSION : 1.0
00008 !
00009 ! MODULE: marquardt.f90 function marquardt
00039 subroutine marquardt(b,m,ni,v,rl,istop,maxiter,dd)
00040 
00041     use ToleranceThreshold
00042     use VectorLength
00043     USE m_vpropre
00044     use mpimod
00045 
00046     IMPLICIT NONE
00047     !   variables globales
00048     integer,intent(in) :: m,maxiter
00049     integer,intent(inout)::ni,istop
00050     double precision,dimension(m*(m+3)/2),intent(out)::v
00051     double precision,intent(out)::rl
00052     double precision,dimension(m),intent(inout)::b
00053     double precision,intent(out)::dd
00054     double precision::ca,cb
00055 
00056     !   variables locales
00057     integer::nql,ii,nfmax,idpos,ncount,id,jd,j,i,ij,integer,ier
00058     double precision,dimension(m*(m+3)/2)::fu
00059     double precision,dimension(m)::delta,b1,bh
00060     double precision::da,dm,ga,tr
00061     double precision::GHG,det,step,eps,vw,fi,maxt, 
00062         z,rl1,th,ep,epsa,epsb,epsd
00063     double precision,external::funcpaRandomEffect
00064 
00065 
00066     epsa=marquardt_ca
00067     epsb=marquardt_cb
00068     epsd=marquardt_rdm
00069     eps=gonfle_eps
00070     da=gonfle_da
00071     dm=gonfle_dm
00072     ep=dsinv_ep
00073     th=searpas_th
00074 
00075     id=0
00076     jd=0
00077     z=0.d0
00078 
00079 
00080     nfmax=m*(m+1)/2
00081     ca=epsa+1.d0
00082     cb=epsb+1.d0
00083     rl1=-1.d+10
00084     ni=0
00085     istop=0
00086 
00087     nql=1
00088 
00089 
00090     Main:Do
00091         if(debugMARC)write(0,*)
00092         if(debugMARC)write(0,*)'***     Iteration',ni,'      ***'
00093         if(debugMARC) write(0,*)'B=',(b(j),j=1,m)
00094         call derivRANDOMEFFECT(b,m,v,rl)
00095         if(debugMARC)write(0,*)'Log-likelihood=',rl
00096 
00097 
00098         rl1=rl
00099 
00100         dd = 0.d0
00101         fu=0.D0
00102         do i=1,m
00103             do j=i,m
00104                 ij=(j-1)*j/2+i
00105                 fu(ij)=v(ij)
00106             end do
00107         end do
00108 
00109         if(debugMARC)write(0,*)"mat. hessian",fu
00110 
00111         ier=0
00112         call dsinv(fu,m,ep,ier,det)
00113         if (ier.eq.-1) then
00114             dd=epsd+1.d0
00115         else
00116             GHG = 0.d0
00117             do i=1,m
00118                 do j=1,m
00119                     if(j.ge.i) then
00120                         ij=(j-1)*j/2+i
00121                     else
00122                         ij=(i-1)*i/2+j
00123                     end if
00124                     GHG = GHG + v(nfmax+i)*fu(ij)*V(nfmax+j)
00125                 end do
00126             end do
00127             dd=GHG/dble(m)
00128         end if
00129 
00130         if(ca.lt.epsa.and.cb.lt.epsb.and.dd.lt.epsd) exit main
00131 
00132         tr=0.d0
00133         do i=1,m
00134             ii=i*(i+1)/2
00135             tr=tr+dabs(v(ii))
00136         end do
00137         tr=tr/dble(m)
00138 
00139         ncount=0
00140         ga=gonfle_ga
00141 400 continue
00142     do i=1,nfmax+m
00143         fu(i)=v(i)
00144     end do
00145     do i=1,m
00146         ii=i*(i+1)/2
00147         if (v(ii).ne.0) then
00148             fu(ii)=v(ii)+da*((1.d0-ga)*dabs(v(ii))+ga*tr)
00149         else
00150             fu(ii)=da*ga*tr
00151         endif
00152     end do
00153     call dchole(fu,m,nql,idpos)
00154     if (idpos.ne.0) then
00155         ncount=ncount+1
00156         if (ncount.le.3.or.ga.ge.1.d0) then
00157             da=da*dm
00158         else
00159             ga=ga*dm
00160             if (ga.gt.1.d0) ga=1.d0
00161         endif
00162         if (ncount .LE. 100000)goto 400
00163     else
00164         vw=1.0D0
00165         do i=1,m
00166             delta(i)=fu(nfmax+i)
00167             b1(i)=b(i)+delta(i)
00168         end do
00169         rl=funcpaRandomEffect(b1,m,id,z,jd,z)
00170         if (rl1.lt.rl) then
00171             if(da.lt.eps) then
00172                 da=eps
00173             else
00174                 da=da/(dm+2.d0)
00175             endif
00176             goto 800
00177         endif
00178     endif
00179     call dmaxt(maxt,delta,m)
00180     if(maxt.eq.0.D0) then
00181         vw=th
00182     else
00183         vw=th/maxt
00184     endif
00185     step=dlog(1.5d0)
00186     call searpas(vw,step,b,bh,m,delta,fi,funcpaRandomEffect)
00187     rl=-fi
00188     if(rl.eq.-1.D9) then
00189         istop=4
00190         goto 110
00191     end if
00192 
00193     do i=1,m
00194         delta(i)=vw*delta(i)
00195     end do
00196     da=(dm-3.d0)*da
00197 
00198 800 cb=dabs(rl1-rl)
00199     ca=0.d0
00200     do i=1,m
00201         ca=ca+delta(i)*delta(i)
00202     end do
00203     if(debugMARC)write(0,*) 'ca =',ca,' cb =',cb,' dd =',dd
00204     if(debugMARC)write(0,*)"delta= ",delta,rl
00205     if(debugMARC)write(0,*)"vw= ",vw
00206     if(debugMARC.and.stepbystep)pause
00207     do i=1,m
00208         b(i)=b(i)+delta(i)
00209     end do
00210 
00211     ni=ni+1
00212     if (ni.ge.maxiter) then
00213         istop=2
00214         if(debugMARC)write(0,*) 'maximum number of iteration reached'
00215         goto 110
00216     end if
00217 End do Main
00218 v=0.D0
00219 v(1:m*(m+1)/2)=fu(1:m*(m+1)/2)
00220 istop=1
00221 
00222 110 continue
00223     return
00224 end subroutine marquardt
00225 
00226