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

likelihood.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: likelihood.f90 function vraisobs
00032 subroutine vraisobs(b,num_pat1,m)
00033     use WorkingSharedValues
00034     use mpimod
00035     implicit none
00036 
00037     integer, intent(in) :: m,num_pat1
00038     double precision, dimension(m),intent(in)::b
00039 
00040     external ::vraistot
00041     INTEGER NDIM2, NF2, MINPTS, MAXPTS, RESTAR, NEVAL, IFAIL
00042     DOUBLE PRECISION EPSABS, EPSREL, RESULT, ABSERR2
00043     DOUBLE PRECISION, dimension(10000) :: WORK
00044 
00045     ndim2=ndim
00046     nf2=nf
00047     b1=b
00048 
00049     MINPTS=likelihood_MINPTS
00050     MAXPTS=likelihood_MAXPTS
00051     EPSABS=likelihood_EPSABS
00052     epsrel=abserrfuncpa !1d-3! *** depend on the iteration and the individual likelihood value
00053 
00054     RESTAR=0
00055     adaptive2=.TRUE.
00056     numpat1=num_pat1
00057     call  hrmsym( NDIM2, NF2, MINPTS, MAXPTS, vraistot, EPSABS, &
00058         EPSREL, RESTAR, RESULT, ABSERR2, NEVAL, IFAIL, WORK)
00059     if (result.NE.result) then
00060         vrais_obs(numpat1)=1.d-300
00061     else
00062         vrais_obs(numpat1)=result
00063     end if
00064     likelihoodPRECISION(3,numpat1)=ABSERR2
00065     IF(IFAIL.NE.0)likelihoodERROR(3,numpat1)=1
00066     return
00067 end subroutine VRAISOBS
00068 
00069 
00070 
00071 !------------------------------------------------------------------------------
00072 ! N.I.M.R.O.D. - Normal approximation Inference in Models with Random
00073 !                effects based on Ordinary Differential equations
00074 !------------------------------------------------------------------------------
00075 !
00076 ! VERSION : 1.0
00077 !
00078 ! MODULE: likelihood.f90 function vraistot
00104 subroutine vraistot(NDIM2,X,NF2,FUNVLS)
00105     use WorkingSharedValues
00106     use Constante
00107     use mpimod
00108     implicit none
00109 
00110     integer,intent(in)::ndim2,nf2
00111     double precision,dimension(ndim2),intent(in)::X
00112     double precision,intent(out) :: funvls
00113 
00114     double precision ::vrais_tot,funvls2,funvls1,phid
00115     double precision, dimension(tdef2,npmcomp)::simul,deriv
00116     double precision,dimension(ndim2)::Xprime,Xprime2
00117     integer :: i,j,k
00118 
00119 
00120     funvls1=0.d0
00121     funvls2=0.d0
00122     systeme=1
00123     xprime=0.0D0
00124     if(debugVRAISTOT)print*,"avant solveur",adaptive2
00125     if (adaptive2) then
00126         Xprime(1:ndim)=startsauv(1:ndim,numpat1) + &
00127             matmul(scaleinvsauv(1:ndim,1:ndim,numpat1),X(1:ndim))
00128         if(debugVRAISTOT)print*,"startsauv(1:ndim,numpat1)",startsauv(1:ndim,numpat1)
00129         if(debugVRAISTOT)print*,"scaleinvsauv(1:ndim,1:ndim,numpat1)",scaleinvsauv(1:ndim,1:ndim,numpat1)
00130         if(debugVRAISTOT)print*,"X(1:ndim)",X(1:ndim)
00131 
00132         call solution(xprime,simul,deriv)
00133     else
00134         call solution(x,simul,deriv)
00135     end if
00136     if(debugVRAISTOT)print*,"apres solveur",x,xprime
00137 
00138     deriv=0.d0
00139     vrais_tot=1.d0
00140     do k=1,npmcomp
00141         do j=1,nbobs(numpat1,k)
00142             if (censor(numpat1,j,k).NE.1) then
00143                 funvls1=0.d0
00144                 funvls2=0.d0
00145                 funvls1=donnees(numpat1,j,k)-simul(schedule(numpat1,j,k)+1,k)
00146                 funvls1 = (funvls1 /b1(npm-npmcomp+k))**2
00147                 if(debugVRAISTOT) then
00148                     print*,"vraitot",idpat(numpat1),schedule(numpat1,j,k),donnees(numpat1,j,k),simul(schedule(numpat1,j,k)+1,k)
00149                 end if
00150 
00151                 funvls2=-funvls1
00152                 funvls2=funvls2/2
00153                 funvls2=exp(funvls2)/abs(b1(npm-npmcomp+k))
00154                 funvls2=funvls2/sqrt(2.d0*pigrec)
00155             else
00156                 if(debugVRAISTOT) then
00157                     if((numproc.EQ.0).and.(LOG_PRINT))write(0,*)"vraitot SS  seuil detec",donnees(numpat1,j,k), &
00158                         "detec ",donnees(numpat1,j,k)
00159                 end if
00160                 funvls2=phid((donnees(numpat1,j,k)-simul(schedule(numpat1,j,k)+1,k))/DABS(b1(npm-npmcomp+k)))
00161             end if
00162             vrais_tot=vrais_tot*funvls2
00163         end do
00164     end do
00165     funvls=vrais_tot
00166 
00167     if(debugVRAISTOT)print*,"avant adaptive"
00168     if (adaptive2) then
00169         Xprime2(1:ndim)=matmul(scaleinvsauv(1:ndim,1:ndim,numpat1),X(1:ndim))
00170         Xprime(1:ndim)=startsauv(1:ndim,numpat1) + Xprime2(1:ndim)
00171         funvls=funvls*exp(-(sum(startsauv(1:ndim,numpat1)**2) )/2)/ &
00172             detersauv(numpat1)
00173         funvls=funvls*exp(-(dot_product(Xprime2(1:ndim),startsauv(1:ndim,numpat1))))
00174         funvls=funvls*exp((5.d-1)*dot_product(X(1:ndim), &
00175             matmul(scaleinv2sauv(1:ndim,1:ndim,numpat1),X(1:ndim))))
00176     end if
00177     if(debugVRAISTOT)print*,"apres adaptive => sortir de vraistot"
00178     if (abs(funvls).LE.1.d-300) funvls=1.d-300
00179     if (funvls.NE.funvls) funvls=1.d-300
00180     if(debugVRAISTOT) then
00181         print*,"apres calcul vraistot = ",funvls
00182         pause
00183     end if
00184     return
00185 end subroutine vraistot
00186 
00187 !------------------------------------------------------------------------------
00188 ! N.I.M.R.O.D. - Normal approximation Inference in Models with Random
00189 !                effects based on Ordinary Differential equations
00190 !------------------------------------------------------------------------------
00191 !
00192 ! VERSION : 1.0
00193 !
00194 ! MODULE: likelihood.f90 function vraistotEXP
00222 subroutine vraistotEXP(NDIM2,X,NF2,FUNVLS)
00223     use WorkingSharedValues
00224     use Constante
00225     use mpimod
00226     implicit none
00227 
00228     integer,intent(in)::ndim2,nf2
00229     double precision,dimension(ndim2),intent(in)::X
00230     double precision,intent(out) :: funvls
00231 
00232     double precision ::vrais_tot,funvls2,funvls1,phid
00233     double precision, dimension(tdef2,npmcomp)::simul,deriv
00234     double precision,dimension(ndim2)::Xprime,Xprime2
00235     integer :: i,j,k
00236 
00237 
00238     funvls1=0.d0
00239     funvls2=0.d0
00240     systeme=1
00241     xprime=0.0D0
00242 
00243     if(debugVRAISTOT)print*,"avant solveur",adaptive2
00244     if (adaptive2) then
00245         Xprime(1:ndim)=startsauv(1:ndim,numpat1) + &
00246             matmul(scaleinvsauv(1:ndim,1:ndim,numpat1),X(1:ndim))
00247         if(debugVRAISTOT)print*,"startsauv(1:ndim,numpat1)",startsauv(1:ndim,numpat1)
00248         if(debugVRAISTOT)print*,"scaleinvsauv(1:ndim,1:ndim,numpat1)",scaleinvsauv(1:ndim,1:ndim,numpat1)
00249         if(debugVRAISTOT)print*,"X(1:ndim)",X(1:ndim)
00250 
00251         call solution(xprime,simul,deriv)
00252     else
00253         call solution(x,simul,deriv)
00254     end if
00255     if(debugVRAISTOT)print*,"apres solveur",x,xprime
00256     deriv=0.d0
00257     vrais_tot=0.d0
00258     do k=1,npmcomp
00259         do j=1,nbobs(numpat1,k)
00260             if (censor(numpat1,j,k).NE.1) then
00261                 if(debugVRAISTOT) then
00262                     print*,"vraitotexp",idpat(numpat1),schedule(numpat1,j,k),donnees(numpat1,j,k),simul(schedule(numpat1,j,k)+1,k)
00263                 end if
00264 
00265                 funvls1=0.d0
00266                 funvls2=0.d0
00267                 funvls1=donnees(numpat1,j,k)-simul(schedule(numpat1,j,k)+1,k)
00268                 funvls1 = (funvls1 /b1(npm-npmcomp+k))**2
00269 
00270                 funvls2=-funvls1
00271                 funvls2=funvls2/2
00272 
00273                 funvls2=dexp(funvls2)/dabs(b1(npm-npmcomp+k))
00274                 funvls2=funvls2/dsqrt(2.d0*pigrec)
00275                 funvls2=dlog(funvls2)
00276             else
00277                 funvls2=phid((donnees(numpat1,j,k)-simul(schedule(numpat1,j,k)+1,k))/DABS(b1(npm-npmcomp+k)))
00278                 if(debugVRAISTOT) then
00279                     if((numproc.EQ.0).and.(LOG_PRINT))write(0,*)"vraistotexp SS  seuil detec: pat",numpat1,"time ",j,"comp ", &
00280                         k,donnees(numpat1,j,k),"donnee",simul(schedule(numpat1,j,k)+1,k), &
00281                         "detec ",donnees(numpat1,j,k),"funvls2",funvls2
00282                 end if
00283 
00284                 funvls2=dlog(funvls2)
00285             end if
00286             vrais_tot=vrais_tot+funvls2
00287         end do
00288     end do
00289     funvls=vrais_tot
00290 
00291 
00292     if(debugVRAISTOT)print*,"avant adaptive"
00293     if (adaptive2) then
00294         Xprime2(1:ndim)=matmul(scaleinvsauv(1:ndim,1:ndim,numpat1),X(1:ndim))
00295         Xprime(1:ndim)=startsauv(1:ndim,numpat1) + Xprime2(1:ndim)
00296         funvls=funvls*exp(-(sum(startsauv(1:ndim,numpat1)**2) )/2)/ &
00297             detersauv(numpat1)
00298         funvls=funvls*exp(-(dot_product(Xprime2(1:ndim),startsauv(1:ndim,numpat1))))
00299         funvls=funvls*exp((5.d-1)*dot_product(X(1:ndim), &
00300             matmul(scaleinv2sauv(1:ndim,1:ndim,numpat1),X(1:ndim))))
00301     end if
00302     if(debugVRAISTOT)print*,"apres adaptive => sortir de vraistot"
00303     if (funvls.LE.-1.d9) funvls=-1.d9
00304     if (funvls.NE.funvls) funvls=-1.d9
00305     if(debugVRAISTOT) then
00306         print*,"apres calcul vraistot = ",funvls
00307         pause
00308     end if
00309     return
00310 end subroutine vraistotEXP
00311 
00312 
00313 
00314 !------------------------------------------------------------------------------
00315 ! N.I.M.R.O.D. - Normal approximation Inference in Models with Random
00316 !                effects based on Ordinary Differential equations
00317 !------------------------------------------------------------------------------
00318 !
00319 ! VERSION : 1.0
00320 !
00321 ! MODULE: likelihood.f90 function ODEschedule
00340 subroutine  ODEschedule
00341     use WorkingSharedValues
00342     implicit none
00343     integer :: i,j,k,iii,jjj,kkk,compteur
00344     listtime=0.0D0
00345     tailletime=0
00346     DO iii=1,nbpatient
00347         compteur=1
00348         listtime(iii,compteur)=schedule(iii,1,1)
00349         compteur=compteur+1
00350         DO kkk=1,npmcomp
00351             DO jjj=1,nbobs(iii,kkk)
00352                 if(.NOT.(ANY(listtime(iii,:).eq.schedule(iii,jjj,kkk)))) then
00353                     listtime(iii,compteur)=schedule(iii,jjj,kkk)
00354                     compteur=compteur+1
00355                 end if
00356             END DO
00357         END DO
00358         tailletime(iii)=compteur-1
00359     END DO
00360     return
00361 end subroutine  ODEschedule
00362