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