C loglik.odds.gradient(param,dataset,print.level=0)
subroutine bgd2m(gbeta,glpsi1,glpsi2,beta,lpsi,npar,x,y,
*theta,work,dbeta,dbeta1,dbeta2,n,der,db,db1,db2)
implicit double precision (a-h,o-z)
dimension x(n,npar-2),beta(npar-2),work(n),theta(n),tpr(2),
*tpr1(4),der1(5),der2(5),psi(2),gbeta(npar-2),dbeta(npar-2),
*dbeta1(npar-2),dbeta2(npar-2),y(n),lpsi(2),
*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),
*db(3,npar-2),db1(4,npar-2),db2(5,npar-2),der(npar-2),
*P7(4,4), P8(4,4),dth(3),dth1(4),dth2(5)
double precision lpsi
integer y,k,k1,k2,i,i0,i1,i2,j,j1,j2,npar,m,n0,iaux,n
psi(1) = dexp(lpsi(1))
psi(2) = dexp(lpsi(2))
psi1 = psi(1)
psi2 = psi(2)
ps1 = psi1-1
ps2 = psi2-1
call matp(x,beta,work,n,npar-2,1)
do 10 k=1,n
theta(k) = 1/(1+dexp(-work(k)))
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
p = theta(i0)
do 40 k2=1,(npar-2)
dbeta2(k2) = p*(1-p)*x(i0,k2)
gbeta(k2)=(y(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 (y(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 = y(i0)
th1 = theta(i1)
th2 = theta(i1-1)
call mcpj(th1,th2,psi1,tpr)
p=tpr(y(i0)+1)
call deriv1(theta,psi1,psi2,n,i,j,der1)
dpth = der1(1)
dpth1 = der1(2)
dpth2 = der1(3)
dppsi1 = der1(4)
dppsi2 = der1(5)
dp = y(i1)/(p*(1-p))-1/(1-p)
do 50 k2=1,(npar-2)
dbeta1(k2) = theta(i)*(1-theta(i))*x(i,k2)
dbeta2(k2) = theta(i-1)*(1-theta(i-1))*x(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)
th01 = theta(i0+1)
th0 = theta(i0)
th1 = theta(i1)
call mcpj (th01,th0,psi1,tpr)
call mat2 (tpr(1),tpr(2),P1)
call mcpj (th1,th01,psi1,tpr)
call mat2 (tpr(1),tpr(2),P2)
call matp(P1,P2,P0,2,2,2)
j = y(i0)
tpr(1)= P0(1,2)
tpr(2)= P0(2,2)
p=tpr(y(i0)+1)
dp = y(i1)/(p*(1-p))-1/(1-p)
do 60 k1=0,2
do 70 k2=1,(npar-2)
db(k1+1,k2)= theta(k1+i0)*(1-theta(k1+i0))*x(k1+i0,k2)
70 continue
60 continue
call deriv1(theta,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(theta,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(theta,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
i=i1+1
100 if (i.le.n0) then
i2 = i
110 if (y(i2).eq.(-1)) then
i2=i2+1
go to 110
end if
if (i2.eq.i.and.i1.eq.(i0+1)) then
j = y(i0)
j1 = y(i1)
th = theta(i2)
th1 = theta(i2-1)
th2 = theta(i2-2)
call mcpij(th,th1,th2,psi1,psi2,tpr1)
call deriv2(theta,psi1,psi2,n,i2,j,j1,der2)
dpth = der2(1)
dpth1 = der2(2)
dpth2 = der2(3)
dppsi1 = der2(4)
dppsi2 = der2(5)
p=tpr1(y(i0)+2*y(i1)+1)
dp = y(i2)/(p*(1-p))-1/(1-p)
do 120 k2=1,(npar-2)
dbeta (k2) = theta(i2)*(1-theta(i2))*x(i2,k2)
dbeta1(k2) = theta(i2-1)*(1-theta(i2-1))*x(i2-1,k2)
dbeta2(k2) = theta(i2-2)*(1-theta(i2-2))*x(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 = theta(i0+1)
th0 = theta(i0)
th = theta(i2)
th1 = theta(i2-1)
th2 = theta(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*y(i1)+1)
Pr(2,1)=tpr1(2*y(i1)+2)
call matp(P1,Pr,Pr1,2,2,1)
p=Pr1(y(i0)+1,1)
dp = y(i2)/(p*(1-p))-1/(1-p)
j = y(i0)
j1 = y(i1)
do 130 k1=0,3
do 140 k2=1,(npar-2)
db1(k1+1,k2)= theta(k1+i0)*(1-theta(k1+i0))*x(k1+i0,k2)
140 continue
130 continue
call deriv1(theta,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(theta,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(theta,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 = theta(i0+2)
th01 = theta(i0+1)
th0 = theta(i0)
th12 = theta(i1+2)
th11 = theta(i1+1)
th22 = theta(i1)
call mcpij (th02,th01,th0,psi1,psi2,tpr1)
call mat4 (tpr1(1),tpr1(3),tpr1(2),tpr1(4),P5)
call mcpij (th12,th11,th22,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 = theta(i3)
th1 = theta(i3-1)
th2 = theta(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*y(i0)+y(i1)+1,2*y(i2)+y(i3)+1)
dp = 1/pstar
j = y(i0)
j1 = y(i1)
j2 = y(i2)
do 230 k1=0,4
do 240 k2=1,(npar-2)
db2(k1+1,k2)= theta(k1+i0)*(1-theta(k1+i0))*x(k1+i0,k2)
240 continue
230 continue
call deriv2(theta,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)*y(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)*y(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)*y(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)*y(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)*y(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)*y(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)*y(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)*y(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)*y(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)*y(i3)))
call deriv2(theta,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)*y(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)*y(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)*y(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)*y(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)*y(i3)))
call deriv2(theta,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)*y(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)*y(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)*y(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)*y(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)*y(i3)))
call deriv2(theta,psi1,psi2,n,i2+1,0,1,der2)
dth2(3)=dth2(3)+der2(3)*((1-P5(2*j+j1+1,2*j1+2))*(2*y(i3)-1)*
*(1-P6(2*j1+1,2)+(2*P6(2*j1+1,2)-1)*y(i2)))
dth2(4)=dth2(4)+der2(2)*((1-P5(2*j+j1+1,2*j1+2))*(2*y(i3)-1)*
*(1-P6(2*j1+1,2)+(2*P6(2*j1+1,2)-1)*y(i2)))
dth2(5)=der2(1)*((1-P5(2*j+j1+1,2*j1+2))*(2*y(i3)-1)*
*(1-P6(2*j1+1,2)+(2*P6(2*j1+1,2)-1)*y(i2)))
dpsi1=dpsi1+der2(4)*((1-P5(2*j+j1+1,2*j1+2))*(2*y(i3)-1)*
*(1-P6(2*j1+1,2)+(2*P6(2*j1+1,2)-1)*y(i2)))
dpsi2=dpsi2+der2(5)*((1-P5(2*j+j1+1,2*j1+2))*(2*y(i3)-1)*
*(1-P6(2*j1+1,2)+(2*P6(2*j1+1,2)-1)*y(i2)))
call deriv2(theta,psi1,psi2,n,i2+1,1,1,der2)
dth2(3)=dth2(3)+der2(3)*(P5(2*j+j1+1,2*j1+2)*(2*y(i3)-1)*
*(1-P6(2*j1+2,4)+(2*P6(2*j1+2,4)-1)*y(i2)))
dth2(4)=dth2(4)+der2(2)*(P5(2*j+j1+1,2*j1+2)*(2*y(i3)-1)*
*(1-P6(2*j1+2,4)+(2*P6(2*j1+2,4)-1)*y(i2)))
dth2(5)=dth2(5)+der2(1)*(P5(2*j+j1+1,2*j1+2)*(2*y(i3)-1)*
*(1-P6(2*j1+2,4)+(2*P6(2*j1+2,4)-1)*y(i2)))
dpsi1=dpsi1+der2(4)*(P5(2*j+j1+1,2*j1+2)*(2*y(i3)-1)*
*(1-P6(2*j1+2,4)+(2*P6(2*j1+2,4)-1)*y(i2)))
dpsi2=dpsi2+der2(5)*(P5(2*j+j1+1,2*j1+2)*(2*y(i3)-1)*
*(1-P6(2*j1+2,4)+(2*P6(2*j1+2,4)-1)*y(i2)))
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 = theta(i0+2)
th01= theta(i0+1)
th0 = theta(i0)
th12 = theta(i1+2)
th11 = theta(i1+1)
th1 = theta(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)
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*y(i0)+y(i1)+1,1)
dp = y(i2)/(p*(1-p))-1/(1-p)
j = y(i0)
j1 = y(i1)
do 330 k1=0,3
do 340 k2=1,(npar-2)
db1(k1+1,k2)= theta(k1+i0)*(1-theta(k1+i0))*x(k1+i0,k2)
340 continue
330 continue
call deriv2(theta,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(theta,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(theta,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
