subroutine mblik1(logL,pij,beta,lpsi,npar,x,y,theta,work,n) implicit double precision (a-h,o-z) dimension x(n,npar-1),beta(npar-1),work(n),theta(n),pij(n), *tpr(2),y(n),P0(2,2), P1(2,2), P2(2,2) double precision logL,lpsi,pij integer y,i0,i1,npar,n,n0,k psi= dexp(lpsi) psi1 = psi ps1 = psi1-1 call matp(x,beta,work,n,npar-1,1) do 10 i=1,n theta(i) = 1/(1+dexp(-work(i))) 10 continue i0 = 1 20 if (y(i0).eq.(-1)) then i0=i0+1 go to 20 end if n0 = n 30 if (y(n0).eq.(-1)) then n0=n0-1 go to 30 end if logL = 0 p = theta(i0) pij(i0)=p logL = y(i0)*dlog(p/(1-p))+dlog(1-p) if (i0.eq.n0) return i = i0+1 40 if (i.le.n0) then i1 = i 50 if (y(i1).eq.(-1)) then i1=i1+1 go to 50 end if C i0 is the most recent (past) observation time C i1 is the next observation time if (i1.eq.i) then th1 = theta(i) th2 = theta(i-1) call mcpj(th1,th2,psi1,tpr) else C (one or more intermediate missing datum) call mat2 (0.0D0,1.0D0,P0) do 60 k=(i0+1),i1 th1 = theta(k) th2 = theta(k-1) call mcpj (th1,th2,psi1,tpr) call mat2 (tpr(1),tpr(2),P1) call matp(P0,P1,P2,2,2,2) call matc(P2,P0,2,2) 60 continue tpr(1)= P0(1,2) tpr(2)= P0(2,2) end if p=tpr(y(i0)+1) pij(i1)=p logL = logL+y(i1)*dlog(p/(1-p))+dlog(1-p) i0=i1 i=i0+1 go to 40 end if return end