GJRM
Generalised Joint Regression Modelling

Package index
Search the GJRM package
Vignettes
Functions
448
Source code
220
Man pages
93
Home
/
CRAN
/
GJRM
/
R/bprobgHsDiscr2ROY.r

R/bprobgHsDiscr2ROY.r
In GJRM: Generalised Joint Regression Modelling

Defines functions bprobgHsDiscr2ROY

Documented in bprobgHsDiscr2ROY 

bprobgHsDiscr2ROY <- function(params, respvec, VC, ps, AT = FALSE){

 
  p1 <- p2 <- pdf1 <- pdf2 <- c.copula.be2 <- c.copula.be1 <- l.par <- dl.dbe1 <- d2l.be1.be1 <- l.parL <- NA
  etad1 <- etad2 <- etas1 <- etas2 <- etan1 <- etan2 <- l.ln <- NULL 
 
  ########
 
    eta1 <-         VC$X1%*%params[1:VC$X1.d2]
    eta2 <- eta.tr( VC$X2%*%params[(VC$X1.d2+1):(VC$X1.d2+VC$X2.d2)],                   VC$margins[2])
    eta3 <- eta.tr( VC$X3%*%params[(VC$X1.d2+VC$X2.d2+1):(VC$X1.d2+VC$X2.d2+VC$X3.d2)], VC$margins[3]) 
 
    sigma1.st <- VC$X4%*%params[(VC$X1.d2+VC$X2.d2+VC$X3.d2+1):(VC$X1.d2+VC$X2.d2+VC$X3.d2+VC$X4.d2)]
    sigma2.st <- VC$X5%*%params[(VC$X1.d2+VC$X2.d2+VC$X3.d2+VC$X4.d2+1):(VC$X1.d2+VC$X2.d2+VC$X3.d2+VC$X4.d2+VC$X5.d2)]   
 
    teta.st1 <- VC$X6%*%params[(VC$X1.d2+VC$X2.d2+VC$X3.d2+VC$X4.d2+VC$X5.d2+1):(VC$X1.d2+VC$X2.d2+VC$X3.d2+VC$X4.d2+VC$X5.d2 + VC$X6.d2)]
    teta.st2 <- VC$X7%*%params[(VC$X1.d2+VC$X2.d2+VC$X3.d2+VC$X4.d2+VC$X5.d2 + VC$X6.d2+1):(VC$X1.d2+VC$X2.d2+VC$X3.d2+VC$X4.d2+VC$X5.d2 + VC$X6.d2+VC$X7.d2)]
 

    pd1 <- probm(eta1, VC$margins[1], bc = TRUE, min.dn = VC$min.dn, min.pr = VC$min.pr, max.pr = VC$max.pr) 
        
    p1 <- pd1$pr    #        
    p0 <- 1 - p1    #       
    
    derp1.dereta1      <- pd1$derp1.dereta1  # recall that these derivs are wrt to 1 - p1
    der2p1.dereta1eta1 <- pd1$der2p1.dereta1eta1
        
    s1tr      <- esp.tr(sigma1.st, VC$margins[2])  
    sigma1.st <- s1tr$vrb.st 
    sigma1    <- s1tr$vrb 
    
    s2tr      <- esp.tr(sigma2.st, VC$margins[3])  
    sigma2.st <- s2tr$vrb.st 
    sigma2    <- s2tr$vrb 
    
    VC1   <- list(BivD = VC$BivD1, BivD2 = NULL)  
    resT1 <- teta.tr(VC1, teta.st1)
    
    VC2   <- list(BivD = VC$BivD2, BivD2 = NULL)
    resT2 <- teta.tr(VC2, teta.st2) 
    
    teta.st1 <- resT1$teta.st
    teta.st2 <- resT2$teta.st
    teta1    <- resT1$teta
    teta2    <- resT2$teta 
     
    Cop1 <- VC$BivD1
    Cop2 <- VC$BivD2 
    nC1  <- VC$nC1
    nC2  <- VC$nC2 
      
#######################    

 dHs <- distrHsDiscr(respvec$y2, eta2, sigma1, sigma1.st, nu = 1, nu.st = 1, margin2 = VC$margins[2], naive = FALSE, y2m = VC$y2m, min.dn = VC$min.dn, min.pr = VC$min.pr, max.pr = VC$max.pr)
   
   
 pdf2.M2                         <- dHs$pdf2
 p2.M2                           <- dHs$p2 
 derpdf2.dereta2.M2              <- dHs$derpdf2.dereta2 
 derpdf2.dersigma2.st.M2         <- dHs$derpdf2.dersigma2.st 
 derp2.dersigma.st.M2            <- dHs$derp2.dersigma.st
 derp2.dereta2.M2                <- dHs$derp2.dereta2
 der2p2.dereta2eta2.M2           <- dHs$der2p2.dereta2eta2 
 der2pdf2.dereta2.M2             <- dHs$der2pdf2.dereta2
 der2p2.dersigma2.st2.M2         <- dHs$der2p2.dersigma2.st2
 der2pdf2.dersigma2.st2.M2       <- dHs$der2pdf2.dersigma2.st2
 der2p2.dereta2dersigma2.st.M2   <- dHs$der2p2.dereta2dersigma2.st            
 der2pdf2.dereta2dersigma2.st.M2 <- dHs$der2pdf2.dereta2dersigma2.st
   
   
   
  dHs <- distrHsDiscr(respvec$y3, eta3, sigma2, sigma2.st, nu = 1, nu.st = 1, margin2 = VC$margins[3], naive = FALSE, y2m = VC$y3m, min.dn = VC$min.dn, min.pr = VC$min.pr, max.pr = VC$max.pr)
   
   
 pdf2.M3                         <- dHs$pdf2
 p2.M3                           <- dHs$p2 
 derpdf2.dereta2.M3              <- dHs$derpdf2.dereta2 
 derpdf2.dersigma2.st.M3         <- dHs$derpdf2.dersigma2.st 
 derp2.dersigma.st.M3            <- dHs$derp2.dersigma.st
 derp2.dereta2.M3                <- dHs$derp2.dereta2
 der2p2.dereta2eta2.M3           <- dHs$der2p2.dereta2eta2 
 der2pdf2.dereta2.M3             <- dHs$der2pdf2.dereta2
 der2p2.dersigma2.st2.M3         <- dHs$der2p2.dersigma2.st2
 der2pdf2.dersigma2.st2.M3       <- dHs$der2pdf2.dersigma2.st2
 der2p2.dereta2dersigma2.st.M3   <- dHs$der2p2.dereta2dersigma2.st            
 der2pdf2.dereta2dersigma2.st.M3 <- dHs$der2pdf2.dereta2dersigma2.st
    
    

    
  ########################################################################################################
  
    C1.M2 <- mm(BiCDF(p0[VC$inde0], p2.M2,                                                       nC1, teta1, VC$dof1), min.pr = VC$min.pr, max.pr = VC$max.pr  )
    C2.M2 <- mm(BiCDF(p0[VC$inde0], mm(p2.M2 - pdf2.M2, min.pr = VC$min.pr, max.pr = VC$max.pr), nC1, teta1, VC$dof1), min.pr = VC$min.pr, max.pr = VC$max.pr  )
    
    A.M2  <- mm(C1.M2 - C2.M2, min.pr = VC$min.pr, max.pr = VC$max.pr)
    
    C1.M3 <- mm(BiCDF(p0[VC$inde1], p2.M3,                                                       nC2, teta2, VC$dof2), min.pr = VC$min.pr, max.pr = VC$max.pr  )
    C2.M3 <- mm(BiCDF(p0[VC$inde1], mm(p2.M3 - pdf2.M3, min.pr = VC$min.pr, max.pr = VC$max.pr), nC2, teta2, VC$dof2), min.pr = VC$min.pr, max.pr = VC$max.pr  )
    
    A.M3  <- mm(C1.M3 - C2.M3,  min.pr = VC$min.pr, max.pr = VC$max.pr)    
    B.M3  <- mm(pdf2.M3 - A.M3, min.pr = VC$min.pr, max.pr = VC$max.pr)
  
    l.par[VC$inde0] <- log( A.M2 )
    l.par[VC$inde1] <- log( B.M3 )  
    
    l.par <- VC$weights*l.par 
     
  ########################################################################################################
     

    dH1F.M2 <- copgHs(p0[VC$inde0], p2.M2,                                                       eta1=NULL, eta2=NULL, teta1, teta.st1, Cop1, VC$dof1, min.dn = VC$min.dn, min.pr = VC$min.pr, max.pr = VC$max.pr)  
    dH2F.M2 <- copgHs(p0[VC$inde0], mm(p2.M2 - pdf2.M2, min.pr = VC$min.pr, max.pr = VC$max.pr), eta1=NULL, eta2=NULL, teta1, teta.st1, Cop1, VC$dof1, min.dn = VC$min.dn, min.pr = VC$min.pr, max.pr = VC$max.pr) 
    
    c.copula.be1.C1.M2   <- dH1F.M2$c.copula.be1 
    c.copula.be1.C2.M2   <- dH2F.M2$c.copula.be1 
    
    c.copula.be2.C1.M2   <- dH1F.M2$c.copula.be2 
    c.copula.be2.C2.M2   <- dH2F.M2$c.copula.be2 
    
    c.copula.theta.C1.M2 <- dH1F.M2$c.copula.theta 
    c.copula.theta.C2.M2 <- dH2F.M2$c.copula.theta
    
    derp2m1.dereta2.M2   <- derp2.dereta2.M2 - derpdf2.dereta2.M2
    derp2m1.dersig2.M2   <- derp2.dersigma.st.M2 - derpdf2.dersigma2.st.M2

    
    Cc.M2    <- c.copula.be1.C1.M2 - c.copula.be1.C2.M2  
    C.M2     <- Cc.M2*derp1.dereta1[VC$inde0] # contains already -
    
    Cs.M2    <- c.copula.theta.C1.M2  - c.copula.theta.C2.M2

    Ceta2 <- (c.copula.be2.C1.M2 - c.copula.be2.C2.M2)*derp2.dereta2.M2     + c.copula.be2.C2.M2*derpdf2.dereta2.M2       
    Csig1 <- (c.copula.be2.C1.M2 - c.copula.be2.C2.M2)*derp2.dersigma.st.M2 + c.copula.be2.C2.M2*derpdf2.dersigma2.st.M2  




    dH1F.M3 <- copgHs(p0[VC$inde1], p2.M3,                                                       eta1=NULL, eta2=NULL, teta2, teta.st2, Cop2, VC$dof2, min.dn = VC$min.dn, min.pr = VC$min.pr, max.pr = VC$max.pr)  
    dH2F.M3 <- copgHs(p0[VC$inde1], mm(p2.M3 - pdf2.M3, min.pr = VC$min.pr, max.pr = VC$max.pr), eta1=NULL, eta2=NULL, teta2, teta.st2, Cop2, VC$dof2, min.dn = VC$min.dn, min.pr = VC$min.pr, max.pr = VC$max.pr) 
    
    c.copula.be1.C1.M3   <- dH1F.M3$c.copula.be1 
    c.copula.be1.C2.M3   <- dH2F.M3$c.copula.be1 
    
    c.copula.be2.C1.M3   <- dH1F.M3$c.copula.be2 
    c.copula.be2.C2.M3   <- dH2F.M3$c.copula.be2 
    
    c.copula.theta.C1.M3 <- dH1F.M3$c.copula.theta 
    c.copula.theta.C2.M3 <- dH2F.M3$c.copula.theta
    
    derp2m1.dereta2.M3   <- derp2.dereta2.M3 - derpdf2.dereta2.M3
    derp2m1.dersig2.M3   <- derp2.dersigma.st.M3 - derpdf2.dersigma2.st.M3

    
    Cc.M3    <- c.copula.be1.C1.M3 - c.copula.be1.C2.M3   
    C.M3     <- Cc.M3*derp1.dereta1[VC$inde1] # contains already -
    
    Cs.M3    <- c.copula.theta.C1.M3 - c.copula.theta.C2.M3
    
    
    Ceta3 <- derpdf2.dereta2.M3      - ((c.copula.be2.C1.M3 -  c.copula.be2.C2.M3)*derp2.dereta2.M3     + c.copula.be2.C2.M3*derpdf2.dereta2.M3)
    Csig2 <- derpdf2.dersigma2.st.M3 - ((c.copula.be2.C1.M3 -  c.copula.be2.C2.M3)*derp2.dersigma.st.M3 + c.copula.be2.C2.M3*derpdf2.dersigma2.st.M3)
 
    # ************ check from above what is really needed and not
    

      dl.dbe1[VC$inde0] <-  C.M2/A.M2
      dl.dbe1[VC$inde1] <- -C.M3/B.M3
      dl.dbe1           <-  VC$weights*dl.dbe1

      dl.dbe2       <- VC$weights[VC$inde0]*( Ceta2/A.M2)    
      dl.dsigma1.st <- VC$weights[VC$inde0]*( Csig1/A.M2)
      
      dl.dbe3       <- VC$weights[VC$inde1]*( Ceta3/B.M3)
      dl.dsigma2.st <- VC$weights[VC$inde1]*( Csig2/B.M3)
      
      dl.dteta1.st <- VC$weights[VC$inde0]*( Cs.M2/A.M2)  
      dl.dteta2.st <- VC$weights[VC$inde1]*(-Cs.M3/B.M3)                    

    ######################################################################################################## 
     
      
      c.copula2.be1.C1.M2           <- dH1F.M2$c.copula2.be1
      c.copula2.be1.C2.M2           <- dH2F.M2$c.copula2.be1
      
      c.copula2.be2.C1.M2           <- dH1F.M2$c.copula2.be2
      c.copula2.be2.C2.M2           <- dH2F.M2$c.copula2.be2
      
      c.copula2.be1be2.C1.M2        <- dH1F.M2$c.copula2.be1be2
      c.copula2.be1be2.C2.M2        <- dH2F.M2$c.copula2.be1be2
      
      c.copula2.be2th.C1.M2         <- dH1F.M2$c.copula2.be2th # with star
      c.copula2.be2th.C2.M2         <- dH2F.M2$c.copula2.be2th  
      
      c.copula2.theta.C1.M2         <- dH1F.M2$bit1.th2ATE # no start
      c.copula2.theta.C2.M2         <- dH2F.M2$bit1.th2ATE 
      
      c.copula.thet.C1.M2           <- dH1F.M2$c.copula.thet # NO star
      c.copula.thet.C2.M2           <- dH2F.M2$c.copula.thet    
     
      derteta.derteta.st.M2         <- dH1F.M2$derteta.derteta.st         # does not matter dH1 or dH2 
      der2teta.derteta.stteta.st.M2 <- dH1F.M2$der2teta.derteta.stteta.st   
     
      c.copula2.be1th.C1.M2         <- dH1F.M2$c.copula2.be1th 
      c.copula2.be1th.C2.M2         <- dH2F.M2$c.copula2.be1th    
     
      c.copula2.be1.C1.M3           <- dH1F.M3$c.copula2.be1
      c.copula2.be1.C2.M3           <- dH2F.M3$c.copula2.be1
      
      c.copula2.be2.C1.M3           <- dH1F.M3$c.copula2.be2
      c.copula2.be2.C2.M3           <- dH2F.M3$c.copula2.be2
      
      c.copula2.be1be2.C1.M3        <- dH1F.M3$c.copula2.be1be2
      c.copula2.be1be2.C2.M3        <- dH2F.M3$c.copula2.be1be2
      
      c.copula2.be2th.C1.M3         <- dH1F.M3$c.copula2.be2th
      c.copula2.be2th.C2.M3         <- dH2F.M3$c.copula2.be2th  
      
      c.copula2.theta.C1.M3         <- dH1F.M3$bit1.th2ATE 
      c.copula2.theta.C2.M3         <- dH2F.M3$bit1.th2ATE 
      
      c.copula.thet.C1.M3           <- dH1F.M3$c.copula.thet # NO star
      c.copula.thet.C2.M3           <- dH2F.M3$c.copula.thet    
     
      derteta.derteta.st.M3         <- dH1F.M3$derteta.derteta.st         # does not matter dH1 or dH2 
      der2teta.derteta.stteta.st.M3 <- dH1F.M3$der2teta.derteta.stteta.st   
     
      c.copula2.be1th.C1.M3         <- dH1F.M3$c.copula2.be1th 
      c.copula2.be1th.C2.M3         <- dH2F.M3$c.copula2.be1th   
  
##########  
  
      der2p2m1.dereta2eta2.M2 <- der2p2.dereta2eta2.M2 - der2pdf2.dereta2.M2
      der2p2m1.dereta2eta2.M3 <- der2p2.dereta2eta2.M3 - der2pdf2.dereta2.M3
      
      der2p2m1.dersig2.M2 <- der2p2.dersigma2.st2.M2 - der2pdf2.dersigma2.st2.M2
      der2p2m1.dersig2.M3 <- der2p2.dersigma2.st2.M3 - der2pdf2.dersigma2.st2.M3

##########################

      b1b1CY <- (c.copula2.be1.C1.M2 - c.copula2.be1.C2.M2)*derp1.dereta1[VC$inde0]^2 + Cc.M2*der2p1.dereta1eta1[VC$inde0]
      b1b1Y  <- (c.copula2.be1.C1.M3 - c.copula2.be1.C2.M3)*derp1.dereta1[VC$inde1]^2 + Cc.M3*der2p1.dereta1eta1[VC$inde1]  
      d2l.be1.be1[VC$inde0] <-  b1b1CY/A.M2 - C.M2^2/A.M2^2
      d2l.be1.be1[VC$inde1] <- -b1b1Y/B.M3  - C.M3^2/B.M3^2
      d2l.be1.be1           <- -VC$weights*d2l.be1.be1         
      

           
      b3b3Y  <- der2pdf2.dereta2.M3 - ( c.copula2.be2.C1.M3*derp2.dereta2.M3^2 + c.copula.be2.C1.M3*der2p2.dereta2eta2.M3 - (c.copula2.be2.C2.M3*derp2m1.dereta2.M3^2 + c.copula.be2.C2.M3*der2p2m1.dereta2eta2.M3)  ) 
      d2l.be3.be3 <- -VC$weights[VC$inde1]*( b3b3Y/B.M3  - Ceta3^2/B.M3^2    ) 

             
      b2b2CY <-   c.copula2.be2.C1.M2*derp2.dereta2.M2^2 + c.copula.be2.C1.M2*der2p2.dereta2eta2.M2 - (c.copula2.be2.C2.M2*derp2m1.dereta2.M2^2 + c.copula.be2.C2.M2*der2p2m1.dereta2eta2.M2)   
      d2l.be2.be2 <- -VC$weights[VC$inde0]*( b2b2CY/A.M2 - Ceta2^2/A.M2^2    ) 
 
 
 
 
      b1b3 <-     -(c.copula2.be1be2.C1.M3*derp2.dereta2.M3 - c.copula2.be1be2.C2.M3*derp2m1.dereta2.M3)*derp1.dereta1[VC$inde1]
      d2l.be1.be3 <- -VC$weights[VC$inde1]*(b1b3/B.M3 - -C.M3*Ceta3/B.M3^2  )  
      
      
      
      
      b1t2 <- -(c.copula2.be1th.C1.M3 - c.copula2.be1th.C2.M3)*derp1.dereta1[VC$inde1] 
      d2l.be1.th2 <- -VC$weights[VC$inde1]*( b1t2/B.M3 - -C.M3*-Cs.M3/B.M3^2  ) 
     
     
     
     
     
      t1t1CY <- c.copula2.theta.C1.M2*derteta.derteta.st.M2^2 + c.copula.thet.C1.M2*der2teta.derteta.stteta.st.M2 - (c.copula2.theta.C2.M2*derteta.derteta.st.M2^2 + c.copula.thet.C2.M2*der2teta.derteta.stteta.st.M2 ) 
      d2l.th1.th1 <- -VC$weights[VC$inde0]*( t1t1CY/A.M2 - Cs.M2^2/A.M2^2    )  
     
     
     
     

     
      t2t2Y  <- -(c.copula2.theta.C1.M3*derteta.derteta.st.M3^2 + c.copula.thet.C1.M3*der2teta.derteta.stteta.st.M3 - (c.copula2.theta.C2.M3*derteta.derteta.st.M3^2 + c.copula.thet.C2.M3*der2teta.derteta.stteta.st.M3)  ) 
      d2l.th2.th2 <- -VC$weights[VC$inde1]*(t2t2Y/B.M3 - Cs.M3^2/B.M3^2    ) 





      b1t1 <-  (c.copula2.be1th.C1.M2 - c.copula2.be1th.C2.M2)*derp1.dereta1[VC$inde0] 
      d2l.be1.th1 <- -VC$weights[VC$inde0]*( b1t1/A.M2 - C.M2*Cs.M2/A.M2^2  ) 





      b1b2 <-      (c.copula2.be1be2.C1.M2*derp2.dereta2.M2 - c.copula2.be1be2.C2.M2*derp2m1.dereta2.M2)*derp1.dereta1[VC$inde0]
      d2l.be1.be2 <- -VC$weights[VC$inde0]*( b1b2/A.M2 - C.M2*Ceta2/A.M2^2 ) 
 

 
 
 
      
      

      b2t1 <-   c.copula2.be2th.C1.M2*derp2.dereta2.M2 - c.copula2.be2th.C2.M2*derp2m1.dereta2.M2
      d2l.be2.th1 <- -VC$weights[VC$inde0]*( b2t1/A.M2 -  Ceta2*Cs.M2/A.M2^2 ) 
 
 

      b3t2 <- -(c.copula2.be2th.C1.M3*derp2.dereta2.M3 - c.copula2.be2th.C2.M3*derp2m1.dereta2.M3)
      d2l.be3.th2 <- -VC$weights[VC$inde1]*( b3t2/B.M3 - -Ceta3*Cs.M3/B.M3^2 ) 

















      be1.si1T <-  (c.copula2.be1be2.C1.M2*derp2.dersigma.st.M2 - c.copula2.be1be2.C2.M2*derp2m1.dersig2.M2)*derp1.dereta1[VC$inde0]
      d2l.be1.si1 <- -VC$weights[VC$inde0]*( be1.si1T/A.M2 - C.M2*Csig1/A.M2^2 ) # looks ok



      be1.si2T <- -(c.copula2.be1be2.C1.M3*derp2.dersigma.st.M3 - c.copula2.be1be2.C2.M3*derp2m1.dersig2.M3)*derp1.dereta1[VC$inde1]
      d2l.be1.si2 <- -VC$weights[VC$inde1]*(be1.si2T/B.M3 - -C.M3*Csig2/B.M3^2  ) # looks ok 
      



      be2.si1T <- c.copula2.be2.C1.M2*derp2.dersigma.st.M2*derp2.dereta2.M2 + c.copula.be2.C1.M2*der2p2.dereta2dersigma2.st.M2 - (c.copula2.be2.C2.M2*derp2m1.dereta2.M2*derp2m1.dersig2.M2 + c.copula.be2.C2.M2*( der2p2.dereta2dersigma2.st.M2-der2pdf2.dereta2dersigma2.st.M2))   
      d2l.be2.si1 <- -VC$weights[VC$inde0]*( be2.si1T/A.M2 - Csig1*Ceta2/A.M2^2    ) # DOUBT HERE, checked and fixed
      
      
 

      be3.si2T  <- der2pdf2.dereta2dersigma2.st.M3 - ( c.copula2.be2.C1.M3*derp2.dereta2.M3*derp2.dersigma.st.M3 + c.copula.be2.C1.M3*der2p2.dereta2dersigma2.st.M3 - (c.copula2.be2.C2.M3*derp2m1.dereta2.M3*derp2m1.dersig2.M3 + c.copula.be2.C2.M3*( der2p2.dereta2dersigma2.st.M3-der2pdf2.dereta2dersigma2.st.M3))  ) 
      d2l.be3.si2 <- -VC$weights[VC$inde1]*(be3.si2T/B.M3  - Csig2*Ceta3/B.M3^2    ) # DOUBT HERE, checked and fixed
 
 
 
 
      si1.si1T    <- c.copula2.be2.C1.M2*derp2.dersigma.st.M2^2 + c.copula.be2.C1.M2*der2p2.dersigma2.st2.M2 - (c.copula2.be2.C2.M2*derp2m1.dersig2.M2^2 + c.copula.be2.C2.M2*der2p2m1.dersig2.M2)   
      d2l.si1.si1 <- -VC$weights[VC$inde0]*( si1.si1T/A.M2 - Csig1^2/A.M2^2    ) # looks ok
 
      si2.si2T  <- der2pdf2.dersigma2.st2.M3 - ( c.copula2.be2.C1.M3*derp2.dersigma.st.M3^2 + c.copula.be2.C1.M3*der2p2.dersigma2.st2.M3 - (c.copula2.be2.C2.M3*derp2m1.dersig2.M3^2 + c.copula.be2.C2.M3*der2p2m1.dersig2.M3)  ) 
      d2l.si2.si2 <- -VC$weights[VC$inde1]*( si2.si2T/B.M3  - Csig2^2/B.M3^2    )  # looks ok
 
 
      si1t1T <-   c.copula2.be2th.C1.M2*derp2.dersigma.st.M2 - c.copula2.be2th.C2.M2*derp2m1.dersig2.M2
      d2l.si1.th1 <- -VC$weights[VC$inde0]*( si1t1T/A.M2 -  Csig1*Cs.M2/A.M2^2 )  # looks fine 
  
   
      si2t2T <- -(c.copula2.be2th.C1.M3*derp2.dersigma.st.M3 - c.copula2.be2th.C2.M3*derp2m1.dersig2.M3)
      d2l.si2.th2 <- -VC$weights[VC$inde1]*( si2t2T/B.M3 - -Csig2*Cs.M3/B.M3^2 ) # looks fine 
 
 

 
 
 
 
 
 

  be1.be1 <- crossprod(VC$X1*c(d2l.be1.be1),           VC$X1)
  be1.be2 <- crossprod(VC$X1[VC$inde0,]*c(d2l.be1.be2),VC$X2)
  be1.be3 <- crossprod(VC$X1[VC$inde1,]*c(d2l.be1.be3),VC$X3)
  be1.si1 <- crossprod(VC$X1[VC$inde0,]*c(d2l.be1.si1),VC$X4)
  be1.si2 <- crossprod(VC$X1[VC$inde1,]*c(d2l.be1.si2),VC$X5)
  be1.th1 <- crossprod(VC$X1[VC$inde0,]*c(d2l.be1.th1),VC$X6)
  be1.th2 <- crossprod(VC$X1[VC$inde1,]*c(d2l.be1.th2),VC$X7) 
 
 
  be2.be2 <-    crossprod(VC$X2*c(d2l.be2.be2),VC$X2)
  be2.be3 <- matrix(0,dim(VC$X2)[2],       dim(VC$X3)[2])
  be2.si1 <-    crossprod(VC$X2*c(d2l.be2.si1),VC$X4)
  be2.si2 <- matrix(0,dim(VC$X2)[2],       dim(VC$X5)[2])  
  be2.th1 <-    crossprod(VC$X2*c(d2l.be2.th1),VC$X6)
  be2.th2 <- matrix(0,dim(VC$X2)[2],       dim(VC$X7)[2])  
  
  be3.be3 <-    crossprod(VC$X3*c(d2l.be3.be3),VC$X3)
  be3.si1 <- matrix(0,dim(VC$X3)[2],       dim(VC$X4)[2])
  be3.si2 <-    crossprod(VC$X3*c(d2l.be3.si2),VC$X5) 
  be3.th1 <- matrix(0,dim(VC$X3)[2],       dim(VC$X6)[2])   
  be3.th2 <-    crossprod(VC$X3*c(d2l.be3.th2),VC$X7)    


  
  si1.si1 <-    crossprod(VC$X4*c(d2l.si1.si1),VC$X4)
  si1.si2 <- matrix(0,dim(VC$X4)[2],       dim(VC$X5)[2])
  si1.th1 <-    crossprod(VC$X4*c(d2l.si1.th1),VC$X6)
  si1.th2 <- matrix(0,dim(VC$X4)[2],       dim(VC$X7)[2])    
  
  si2.si2 <-    crossprod(VC$X5*c(d2l.si2.si2),VC$X5)
  si2.th1 <- matrix(0,dim(VC$X5)[2],       dim(VC$X6)[2])
  si2.th2 <-    crossprod(VC$X5*c(d2l.si2.th2),VC$X7)     
  
  th1.th1 <-    crossprod(VC$X6*c(d2l.th1.th1),VC$X6)
  th1.th2 <- matrix(0,dim(VC$X6)[2],       dim(VC$X7)[2])   
  
  th2.th2 <-    crossprod(VC$X7*c(d2l.th2.th2),VC$X7)

    
 



  
  H <- rbind( cbind(   be1.be1 ,   be1.be2 ,   be1.be3 ,   be1.si1 ,   be1.si2 ,   be1.th1 ,  be1.th2 ), 
              cbind( t(be1.be2),   be2.be2 ,   be2.be3 ,   be2.si1 ,   be2.si2 ,   be2.th1 ,  be2.th2 ), 
              cbind( t(be1.be3), t(be2.be3),   be3.be3 ,   be3.si1 ,   be3.si2 ,   be3.th1 ,  be3.th2 ),
              cbind( t(be1.si1), t(be2.si1), t(be3.si1),   si1.si1 ,   si1.si2 ,   si1.th1 ,  si1.th2 ),
              cbind( t(be1.si2), t(be2.si2), t(be3.si2), t(si1.si2),   si2.si2 ,   si2.th1 ,  si2.th2 ),
              cbind( t(be1.th1), t(be2.th1), t(be3.th1), t(si1.th1), t(si2.th1),   th1.th1 ,  th1.th2 ),
              cbind( t(be1.th2), t(be2.th2), t(be3.th2), t(si1.th2), t(si2.th2), t(th1.th2),  th2.th2 )             
            ) 
            
           
  G   <- -c( colSums(c(dl.dbe1      )*VC$X1),
             colSums(c(dl.dbe2      )*VC$X2),
             colSums(c(dl.dbe3      )*VC$X3),
             colSums(c(dl.dsigma1.st)*VC$X4),
             colSums(c(dl.dsigma2.st)*VC$X5),
             colSums(c(dl.dteta1.st )*VC$X6), 
             colSums(c(dl.dteta2.st )*VC$X7) 
             
             ) 


  res <- -sum(l.par)

if(VC$extra.regI == "pC" && VC$hess==FALSE) H <- regH(H, type = 1)
  
  S.h <- ps$S.h
  
  if( length(S.h) != 1){
  
    S.h1 <- 0.5*crossprod(params,S.h)%*%params
    S.h2 <- S.h%*%params
  
  } else S.h <- S.h1 <- S.h2 <- 0   
  
  S.res <- res
  res   <- S.res + S.h1
  G     <- G + S.h2
  H     <- H + S.h 
      
if(VC$extra.regI == "sED") H <- regH(H, type = 2) 
  
  
  
         list(value = res, gradient = G, hessian = H, S.h = S.h, S.h1 = S.h1, S.h2 = S.h2, l = S.res, l.par = l.par, ps = ps, l.ln = NULL,
              eta1 = eta1, eta2 = eta2, eta3 = eta3, etad1 = etad1, 
              etad2 = etad2, etas1 = etas1, etas2 = etas2, etan1 = etan1, etan2 = etan2,
              dl.dbe1 = dl.dbe1, dl.dbe2 = dl.dbe2, dl.dbe3 = dl.dbe3, dl.dsigma1.st = dl.dsigma1.st, dl.dsigma2.st = dl.dsigma2.st,
              dl.dteta1.st = dl.dteta1.st, dl.dteta2.st = dl.dteta2.st,
              BivD1 = VC$BivD1, BivD2 = VC$BivD2,                              
              p1 = p1, p0 = p0, pdf1 = pdf2.M2, pdf2 = pdf2.M3, c.copula2.be1be2 = c(c.copula2.be1be2.C1.M2, c.copula2.be1be2.C2.M2, c.copula2.be1be2.C1.M3, c.copula2.be1be2.C2.M3),         
              teta.st1 = teta.st1, teta.st2 = teta.st2,
              Cop1 = Cop1, Cop2 = Cop2, teta1 = teta1, teta2 = teta2)      

}

Try the GJRM package in your browser

Any scripts or data that you put into this service are public.

GJRM documentation built on July 9, 2023, 7:15 p.m.

R Package Documentation

rdrr.io home R language documentation Run R code online

Browse R Packages

CRAN packages Bioconductor packages R-Forge packages GitHub packages

We want your feedback!

Note that we can't provide technical support on individual packages. You should contact the package authors for that.
GitHub issue tracker
ian@mutexlabs.com
Personal blog

 
Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.

Close