C loglik.odds.gradient(param,dataset,print.level=0) subroutine mbgd2i(gbeta,glpsi1,glpsi2,npar,n) implicit double precision (a-h,o-z) DIMENSION x1(4500,10),theta1(4500), *work1(4500),y1(4500),lpsi1(2),beta1(10),bt1(10), *dbeta(10),dbeta1(10),dbeta2(10), *der(10),db(3,10),db1(4,10),db2(5,10),gbeta(10), *psi(2),tpr(2),tpr1(4),der1(5),der2(5), *P0(2,2), P1(2,2), P2(2,2),P4(4,4), P5(4,4), P6(4,4), *Pc(4,1),Pr(2,1),Pr0(4,1),Pr1(2,1),Paux(2,2),Pres(2,2), *dth(3),dth1(4),dth2(5),P7(4,4), P8(4,4) double precision lpsi1 integer y1,k,k1,k2,i,i0,i1,i2,j,j1,j2,npar,n0,iaux,n, *m,mpar COMMON/param/x1,theta1,work1, *y1,lpsi1,beta1,bt1,m,mpar,omega1 COMMON/grad/ dbeta,dbeta1,dbeta2, *der,db,db1,db2 psi(1) = dexp(lpsi1(1)) psi(2) = dexp(lpsi1(2)) psi1 = psi(1) psi2 = psi(2) ps1 = psi1-1 ps2 = psi2-1 call mati(x1,beta1,work1,4500,10,1,n,npar) do 10 k=1,n theta1(k) = 1/(1+dexp(-work1(k))) 10 continue i0 = 1 20 if (y1(i0).eq.(-1)) then i0=i0+1 go to 20 end if n0 = n 30 if (y1(n0).eq.(-1)) then n0=n0-1 go to 30 end if p = theta1(i0) do 40 k2=1,(npar-2) dbeta2(k2) = p*(1-p)*x1(i0,k2) gbeta(k2)=(y1(i0)/(p*(1-p))-1/(1-p))*dbeta2(k2) 40 continue glpsi1 = 0 glpsi2 = 0 if (i0.eq.n0) return gpsi1 = 0 gpsi2 = 0 i = i0+1 i1 = i 45 if (y1(i1).eq.(-1)) then i1=i1+1 go to 45 end if C i0 is the most recent (past) observation time C i1 is the next observation time if (i1.eq.i) then j = y1(i0) th1 = theta1(i1) th2 = theta1(i1-1) call mcpj(th1,th2,psi1,tpr) p=tpr(y1(i0)+1) call deriv1(theta1,psi1,psi2,n,i,j,der1) dpth = der1(1) dpth1 = der1(2) dpth2 = der1(3) dppsi1 = der1(4) dppsi2 = der1(5) dp = y1(i1)/(p*(1-p))-1/(1-p) do 50 k2=1,(npar-2) dbeta1(k2) = theta1(i)*(1-theta1(i))*x1(i,k2) dbeta2(k2) = theta1(i-1)*(1-theta1(i-1))*x1(i-1,k2) gbeta(k2) = gbeta(k2) + * dp*(dpth1*dbeta1(k2)+dpth2*dbeta2(k2)) 50 continue gpsi1 = gpsi1+dp*dppsi1 gpsi2 = gpsi2+dp*dppsi2 else C (exactly one intermediate missing datum between i0 and i1) call mcpj (theta1(i0+1),theta1(i0),psi1,tpr) call mat2 (tpr(1),tpr(2),P1) call mcpj (theta1(i1),theta1(i0+1),psi1,tpr) call mat2 (tpr(1),tpr(2),P2) call matp(P1,P2,P0,2,2,2) j = y1(i0) tpr(1)= P0(1,2) tpr(2)= P0(2,2) p=tpr(y1(i0)+1) dp = y1(i1)/(p*(1-p))-1/(1-p) do 60 k1=0,2 do 70 k2=1,(npar-2) db(k1+1,k2)= theta1(k1+i0)*(1-theta1(k1+i0))*x1(k1+i0,k2) 70 continue 60 continue call deriv1(theta1,psi1,psi2,n,i0+1,j,der1) dth(1)=der1(3)*(-P2(1,2)+P2(2,2)) dth(2)=der1(2)*(-P2(1,2)+P2(2,2)) dpsi1=der1(4)*(-P2(1,2)+P2(2,2)) dpsi2=der1(5)*(-P2(1,2)+P2(2,2)) call deriv1(theta1,psi1,psi2,n,i1,0,der1) dth(2)=dth(2)+der1(3)*(1-P1(j+1,2)) dth(3)=der1(2)*(1-P1(j+1,2)) dpsi1=dpsi1+der1(4)*(1-P1(j+1,2)) dpsi2=dpsi2+der1(5)*(1-P1(j+1,2)) call deriv1(theta1,psi1,psi2,n,i1,1,der1) dth(2)=dth(2)+der1(3)*P1(j+1,2) dth(3)=dth(3)+der1(2)*P1(j+1,2) dpsi1=dpsi1+der1(4)*P1(j+1,2) dpsi2=dpsi2+der1(5)*P1(j+1,2) do 80 k2=1,(npar-2) der(k2)=0.0D0 do 90 k1=1,3 der(k2)= der(k2)+dth(k1)*db(k1,k2) 90 continue 80 continue do 95 k2=1,(npar-2) gbeta(k2)= gbeta(k2)+dp*der(k2) 95 continue gpsi1 = gpsi1+dp*dpsi1 gpsi2 = gpsi2+dp*dpsi2 end if C if (i1.eq.n0) return i=i1+1 100 if (i.le.n0) then i2 = i 110 if (y1(i2).eq.(-1)) then i2=i2+1 go to 110 end if if (i2.eq.i.and.i1.eq.(i0+1)) then j = y1(i0) j1 = y1(i1) th = theta1(i2) th1 = theta1(i2-1) th2 = theta1(i2-2) call mcpij(th,th1,th2,psi1,psi2,tpr1) call deriv2(theta1,psi1,psi2,n,i2,j,j1,der2) dpth = der2(1) dpth1 = der2(2) dpth2 = der2(3) dppsi1 = der2(4) dppsi2 = der2(5) p=tpr1(y1(i0)+2*y1(i1)+1) dp = y1(i2)/(p*(1-p))-1/(1-p) do 120 k2=1,(npar-2) dbeta (k2) = theta1(i2)*(1-theta1(i2))*x1(i2,k2) dbeta1(k2) = theta1(i2-1)*(1-theta1(i2-1))*x1(i2-1,k2) dbeta2(k2) = theta1(i2-2)*(1-theta1(i2-2))*x1(i2-2,k2) gbeta(k2) = gbeta(k2) + * dp*(dpth*dbeta(k2)+dpth1*dbeta1(k2)+dpth2*dbeta2(k2)) 120 continue gpsi1 = gpsi1+dp*dppsi1 gpsi2 = gpsi2+dp*dppsi2 i0=i1 i1=i2 i=i1+1 else if (i1.ne.(i0+1).and.i2.eq.(i1+1)) then C (exactly one intermediate missing datum between i0 and i1) th01 = theta1(i0+1) th0 = theta1(i0) th = theta1(i2) th1 = theta1(i2-1) th2 = theta1(i2-2) call mcpj (th01,th0,psi1,tpr) call mat2 (tpr(1),tpr(2),P1) call mcpij (th,th1,th2,psi1,psi2,tpr1) call mat4 (tpr1(1),tpr1(3),tpr1(2),tpr1(4),P6) Pr(1,1)=tpr1(2*y1(i1)+1) Pr(2,1)=tpr1(2*y1(i1)+2) call matp(P1,Pr,Pr1,2,2,1) p=Pr1(y1(i0)+1,1) dp = y1(i2)/(p*(1-p))-1/(1-p) j = y1(i0) j1 = y1(i1) do 130 k1=0,3 do 140 k2=1,(npar-2) db1(k1+1,k2)= theta1(k1+i0)*(1-theta1(k1+i0))*x1(k1+i0,k2) 140 continue 130 continue call deriv1(theta1,psi1,psi2,n,i0+1,j,der1) dth1(1)=der1(3)*(-P6(j1+1,2*j1+2)+P6(j1+3,2*j1+2)) dth1(2)=der1(2)*(-P6(j1+1,2*j1+2)+P6(j1+3,2*j1+2)) dpsi1=der1(4)*(-P6(j1+1,2*j1+2)+P6(j1+3,2*j1+2)) call deriv2(theta1,psi1,psi2,n,i2,0,j1,der2) dth1(2)=dth1(2)+der2(3)*(1-P1(j+1,2)) dth1(3)=der2(2)*(1-P1(j+1,2)) dth1(4)=der2(1)*(1-P1(j+1,2)) dpsi1=dpsi1+der2(4)*(1-P1(j+1,2)) dpsi2=der2(5)*(1-P1(j+1,2)) call deriv2(theta1,psi1,psi2,n,i2,1,j1,der2) dth1(2)=dth1(2)+der2(3)*P1(j+1,2) dth1(3)=dth1(3)+der2(2)*P1(j+1,2) dth1(4)=dth1(4)+der2(1)*P1(j+1,2) dpsi1=dpsi1+der2(4)*P1(j+1,2) dpsi2=dpsi2+der2(5)*P1(j+1,2) do 150 k2=1,(npar-2) der(k2)=0.0D0 do 160 k1=1,4 der(k2)= der(k2)+dth1(k1)*db1(k1,k2) 160 continue 150 continue do 170 k2=1,(npar-2) gbeta(k2)= gbeta(k2)+dp*der(k2) 170 continue gpsi1 = gpsi1+dp*dpsi1 gpsi2 = gpsi2+dp*dpsi2 i0=i1 i1=i2 i=i1+1 else if (i2.ne.(i1+1).and.i2.ne.n0) then C (exactly one intermediate missing datum between i1 and i2) th02 = theta1(i0+2) th01 = theta1(i0+1) th0 = theta1(i0) th12 = theta1(i1+2) th11 = theta1(i1+1) th1 = theta1(i1) call mcpij (th02,th01,th0,psi1,psi2,tpr1) call mat4 (tpr1(1),tpr1(3),tpr1(2),tpr1(4),P5) call mcpij (th12,th11,th1,psi1,psi2,tpr1) call mat4 (tpr1(1),tpr1(3),tpr1(2),tpr1(4),P6) call matp(P5,P6,P4,4,4,4) i3=i2+1 th = theta1(i3) th1 = theta1(i3-1) th2 = theta1(i3-2) call mcpij (th,th1,th2,psi1,psi2,tpr1) call mat4 (tpr1(1),tpr1(3),tpr1(2),tpr1(4),P7) call matp(P4,P7,P8,4,4,4) pstar=P8(2*y1(i0)+y1(i1)+1,2*y1(i2)+y1(i3)+1) dp = 1/pstar j = y1(i0) j1 = y1(i1) j2 = y1(i2) do 230 k1=0,4 do 240 k2=1,(npar-2) db2(k1+1,k2)= theta1(k1+i0)*(1-theta1(k1+i0))*x1(k1+i0,k2) 240 continue 230 continue call deriv2(theta1,psi1,psi2,n,i0+2,j,j1,der2) dth2(1)=der2(3)*((P6(2*j1+1,2)-1+(1-2*P6(2*j1+1,2))*j2)* *(1-P7(j2+1,2*j2+2)+(2*P7(j2+1,2*j2+2)-1)*y1(i3))+ *(1-P6(2*j1+2,4)+(2*P6(2*j1+2,4)-1)*j2)* *(1-P7(j2+3,2*j2+2)+(2*P7(j2+3,2*j2+2)-1)*y1(i3))) dth2(2)=der2(2)*((P6(2*j1+1,2)-1+(1-2*P6(2*j1+1,2))*j2)* *(1-P7(j2+1,2*j2+2)+(2*P7(j2+1,2*j2+2)-1)*y1(i3))+ *(1-P6(2*j1+2,4)+(2*P6(2*j1+2,4)-1)*j2)* *(1-P7(j2+3,2*j2+2)+(2*P7(j2+3,2*j2+2)-1)*y1(i3))) dth2(3)=der2(1)*((P6(2*j1+1,2)-1+(1-2*P6(2*j1+1,2))*j2)* *(1-P7(j2+1,2*j2+2)+(2*P7(j2+1,2*j2+2)-1)*y1(i3))+ *(1-P6(2*j1+2,4)+(2*P6(2*j1+2,4)-1)*j2)* *(1-P7(j2+3,2*j2+2)+(2*P7(j2+3,2*j2+2)-1)*y1(i3))) dpsi1=der2(4)*((P6(2*j1+1,2)-1+(1-2*P6(2*j1+1,2))*j2)* *(1-P7(j2+1,2*j2+2)+(2*P7(j2+1,2*j2+2)-1)*y1(i3))+ *(1-P6(2*j1+2,4)+(2*P6(2*j1+2,4)-1)*j2)* *(1-P7(j2+3,2*j2+2)+(2*P7(j2+3,2*j2+2)-1)*y1(i3))) dpsi2=der2(5)*((P6(2*j1+1,2)-1+(1-2*P6(2*j1+1,2))*j2)* *(1-P7(j2+1,2*j2+2)+(2*P7(j2+1,2*j2+2)-1)*y1(i3))+ *(1-P6(2*j1+2,4)+(2*P6(2*j1+2,4)-1)*j2)* *(1-P7(j2+3,2*j2+2)+(2*P7(j2+3,2*j2+2)-1)*y1(i3))) call deriv2(theta1,psi1,psi2,n,i2,j1,0,der2) dth2(2)=dth2(2)+der2(3)*((1-P5(2*j+j1+1,2*j1+2))*(2*j2-1)* *(1-P7(j2+1,2*j2+2)+(2*P7(j2+1,2*j2+2)-1)*y1(i3))) dth2(3)=dth2(3)+der2(2)*((1-P5(2*j+j1+1,2*j1+2))*(2*j2-1)* *(1-P7(j2+1,2*j2+2)+(2*P7(j2+1,2*j2+2)-1)*y1(i3))) dth2(4)=der2(1)*((1-P5(2*j+j1+1,2*j1+2))*(2*j2-1)* *(1-P7(j2+1,2*j2+2)+(2*P7(j2+1,2*j2+2)-1)*y1(i3))) dpsi1=dpsi1+der2(4)*((1-P5(2*j+j1+1,2*j1+2))*(2*j2-1)* *(1-P7(j2+1,2*j2+2)+(2*P7(j2+1,2*j2+2)-1)*y1(i3))) dpsi2=dpsi2+der2(5)*((1-P5(2*j+j1+1,2*j1+2))*(2*j2-1)* *(1-P7(j2+1,2*j2+2)+(2*P7(j2+1,2*j2+2)-1)*y1(i3))) call deriv2(theta1,psi1,psi2,n,i2,j1,1,der2) dth2(2)=dth2(2)+der2(3)*(P5(2*j+j1+1,2*j1+2)*(2*j2-1)* *(1-P7(j2+3,2*j2+2)+(2*P7(j2+3,2*j2+2)-1)*y1(i3))) dth2(3)=dth2(3)+der2(2)*(P5(2*j+j1+1,2*j1+2)*(2*j2-1)* *(1-P7(j2+3,2*j2+2)+(2*P7(j2+3,2*j2+2)-1)*y1(i3))) dth2(4)=dth2(4)+der2(1)*(P5(2*j+j1+1,2*j1+2)*(2*j2-1)* *(1-P7(j2+3,2*j2+2)+(2*P7(j2+3,2*j2+2)-1)*y1(i3))) dpsi1=dpsi1+der2(4)*(P5(2*j+j1+1,2*j1+2)*(2*j2-1)* *(1-P7(j2+3,2*j2+2)+(2*P7(j2+3,2*j2+2)-1)*y1(i3))) dpsi2=dpsi2+der2(5)*(P5(2*j+j1+1,2*j1+2)*(2*j2-1)* *(1-P7(j2+3,2*j2+2)+(2*P7(j2+3,2*j2+2)-1)*y1(i3))) call deriv2(theta1,psi1,psi2,n,i2+1,0,1,der2) dth2(3)=dth2(3)+der2(3)*((1-P5(2*j+j1+1,2*j1+2))*(2*y1(i3)-1)* *(1-P6(2*j1+1,2)+(2*P6(2*j1+1,2)-1)*j2)) dth2(4)=dth2(4)+der2(2)*((1-P5(2*j+j1+1,2*j1+2))*(2*y1(i3)-1)* *(1-P6(2*j1+1,2)+(2*P6(2*j1+1,2)-1)*j2)) dth2(5)=der2(1)*((1-P5(2*j+j1+1,2*j1+2))*(2*y1(i3)-1)* *(1-P6(2*j1+1,2)+(2*P6(2*j1+1,2)-1)*j2)) dpsi1=dpsi1+der2(4)*((1-P5(2*j+j1+1,2*j1+2))*(2*y1(i3)-1)* *(1-P6(2*j1+1,2)+(2*P6(2*j1+1,2)-1)*j2)) dpsi2=dpsi2+der2(5)*((1-P5(2*j+j1+1,2*j1+2))*(2*y1(i3)-1)* *(1-P6(2*j1+1,2)+(2*P6(2*j1+1,2)-1)*j2)) call deriv2(theta1,psi1,psi2,n,i2+1,1,1,der2) dth2(3)=dth2(3)+der2(3)*(P5(2*j+j1+1,2*j1+2)*(2*y1(i3)-1)* *(1-P6(2*j1+2,4)+(2*P6(2*j1+2,4)-1)*j2)) dth2(4)=dth2(4)+der2(2)*(P5(2*j+j1+1,2*j1+2)*(2*y1(i3)-1)* *(1-P6(2*j1+2,4)+(2*P6(2*j1+2,4)-1)*j2)) dth2(5)=dth2(5)+der2(1)*(P5(2*j+j1+1,2*j1+2)*(2*y1(i3)-1)* *(1-P6(2*j1+2,4)+(2*P6(2*j1+2,4)-1)*j2)) dpsi1=dpsi1+der2(4)*(P5(2*j+j1+1,2*j1+2)*(2*y1(i3)-1)* *(1-P6(2*j1+2,4)+(2*P6(2*j1+2,4)-1)*j2)) dpsi2=dpsi2+der2(5)*(P5(2*j+j1+1,2*j1+2)*(2*y1(i3)-1)* *(1-P6(2*j1+2,4)+(2*P6(2*j1+2,4)-1)*j2)) do 250 k2=1,(npar-2) der(k2)=0.0D0 do 260 k1=1,5 der(k2)= der(k2)+dth2(k1)*db2(k1,k2) 260 continue 250 continue do 270 k2=1,(npar-2) gbeta(k2)= gbeta(k2)+dp*der(k2) 270 continue gpsi1 = gpsi1+dp*dpsi1 gpsi2 = gpsi2+dp*dpsi2 i0=i2 i1=i3 i=i1+1 else if (i2.ne.(i1+1).and.i2.eq.n0) then C (exactly one intermediate missing datum between i1 and i2(last value)) th02 = theta1(i0+2) th01 = theta1(i0+1) th0 = theta1(i0) th = theta1(i1+2) th1 = theta1(i1+1) th2 = theta1(i1) call mcpij (th02,th01,th0,psi1,psi2,tpr1) call mat4 (tpr1(1),tpr1(3),tpr1(2),tpr1(4),P5) call mcpij (th,th1,th2,psi1,psi2,tpr1) call mat4 (tpr1(1),tpr1(3),tpr1(2),tpr1(4),P6) call matp(P5,P6,P4,4,4,4) Pc(1,1)=0 Pc(2,1)=1 Pc(3,1)=0 Pc(4,1)=1 call matp(P4,Pc,Pr0,4,4,1) p=Pr0(2*y1(i0)+y1(i1)+1,1) dp = y1(i2)/(p*(1-p))-1/(1-p) j = y1(i0) j1 = y1(i1) do 330 k1=0,3 do 340 k2=1,(npar-2) db1(k1+1,k2)= theta1(k1+i0)*(1-theta1(k1+i0))*x1(k1+i0,k2) 340 continue 330 continue call deriv2(theta1,psi1,psi2,n,i0+2,j,j1,der2) dth1(1)=der2(3)*(-P6(2*j1+1,2)+P6(2*j1+2,4)) dth1(2)=der2(2)*(-P6(2*j1+1,2)+P6(2*j1+2,4)) dth1(3)=der2(1)*(-P6(2*j1+1,2)+P6(2*j1+2,4)) dpsi1=der2(4)*(-P6(2*j1+1,2)+P6(2*j1+2,4)) dpsi2=der2(5)*(-P6(2*j1+1,2)+P6(2*j1+2,4)) call deriv2(theta1,psi1,psi2,n,i2,j1,0,der2) dth1(2)=dth1(2)+der2(3)*(1-P5(2*j+j1+1,2*j1+2)) dth1(3)=dth1(3)+der2(2)*(1-P5(2*j+j1+1,2*j1+2)) dth1(4)=der2(1)*(1-P5(2*j+j1+1,2*j1+2)) dpsi1=dpsi1+der2(4)*(1-P5(2*j+j1+1,2*j1+2)) dpsi2=dpsi2+der2(5)*(1-P5(2*j+j1+1,2*j1+2)) call deriv2(theta1,psi1,psi2,n,i2,j1,1,der2) dth1(2)=dth1(2)+der2(3)*P5(2*j+j1+1,2*j1+2) dth1(3)=dth1(3)+der2(2)*P5(2*j+j1+1,2*j1+2) dth1(4)=dth1(4)+der2(1)*P5(2*j+j1+1,2*j1+2) dpsi1=dpsi1+der2(4)*P5(2*j+j1+1,2*j1+2) dpsi2=dpsi2+der2(5)*P5(2*j+j1+1,2*j1+2) do 350 k2=1,(npar-2) der(k2)=0.0D0 do 360 k1=1,4 der(k2)= der(k2)+dth1(k1)*db1(k1,k2) 360 continue 350 continue do 370 k2=1,(npar-2) gbeta(k2)= gbeta(k2)+dp*der(k2) 370 continue gpsi1 = gpsi1+dp*dpsi1 gpsi2 = gpsi2+dp*dpsi2 i0=i1 i1=i2 i=i1+1 end if go to 100 end if glpsi1=gpsi1*psi1 glpsi2=gpsi2*psi2 return end