H.triESS <- function(params, respvec, VC, TIn, LgTRI){
dst.1 <- dnorm( (TIn$mar2 - TIn$theta12 * TIn$mar1[VC$inde] )/sqrt(1 - TIn$theta12^2) )
pst.1 <- pnorm( ( ((TIn$mar3 - TIn$theta13 * TIn$mar1[VC$inde])/sqrt(1 - TIn$theta13^2)) - ((TIn$theta23 - TIn$theta12 * TIn$theta13)/sqrt((1 - TIn$theta12^2) * (1 - TIn$theta13^2))) * ((TIn$mar2 - TIn$theta12 * TIn$mar1[VC$inde])/sqrt(1 - TIn$theta12^2)) )/sqrt(1 - ((TIn$theta23 - TIn$theta12 * TIn$theta13)/sqrt((1 - TIn$theta12^2) * (1 - TIn$theta13^2)))^2))
pst.1 <- mm(pst.1, min.pr = VC$min.pr, max.pr = VC$max.pr)
st.1 <- -TIn$theta12/sqrt(1 - TIn$theta12^2)
dst.2 <- dnorm((TIn$mar3 - TIn$theta13 * TIn$mar1[VC$inde])/sqrt(1 - TIn$theta13^2))
pst.2 <- pnorm( ( ((TIn$mar2 - TIn$theta12 * TIn$mar1[VC$inde])/sqrt(1 - TIn$theta12^2)) - ((TIn$theta23 - TIn$theta12 * TIn$theta13)/sqrt((1 - TIn$theta12^2) * (1 - TIn$theta13^2))) * ((TIn$mar3 - TIn$theta13 * TIn$mar1[VC$inde])/sqrt(1 - TIn$theta13^2)) )/sqrt(1 - ((TIn$theta23 - TIn$theta12 * TIn$theta13)/sqrt((1 - TIn$theta12^2) * (1 - TIn$theta13^2)))^2))
pst.2 <- mm(pst.2, min.pr = VC$min.pr, max.pr = VC$max.pr)
st.2 <- -TIn$theta13/sqrt(1 - TIn$theta13^2)
dst.3 <- dnorm((TIn$mar1[VC$inde] - TIn$theta12 * TIn$mar2)/sqrt(1 - TIn$theta12^2))
pst.3 <- pnorm( ( ((TIn$mar3 - TIn$theta23 * TIn$mar2)/sqrt(1 - TIn$theta23^2)) - ((TIn$theta13 - TIn$theta12 * TIn$theta23)/sqrt((1 - TIn$theta12^2) * (1 - TIn$theta23^2))) * ((TIn$mar1[VC$inde] - TIn$theta12 * TIn$mar2)/sqrt(1 - TIn$theta12^2)) )/sqrt(1 - ((TIn$theta13 - TIn$theta12 * TIn$theta23)/sqrt((1 - TIn$theta12^2) * (1 - TIn$theta23^2)))^2))
pst.3 <- mm(pst.3, min.pr = VC$min.pr, max.pr = VC$max.pr)
st.3 <- -TIn$theta12/sqrt(1 - TIn$theta12^2)
dst.4 <- dnorm((TIn$mar3 - TIn$theta23 * TIn$mar2)/sqrt(1 - TIn$theta23^2))
pst.4 <- pnorm( ( ((TIn$mar1[VC$inde] - TIn$theta12 * TIn$mar2)/sqrt(1 - TIn$theta12^2)) - ((TIn$theta13 - TIn$theta12 * TIn$theta23)/sqrt((1 - TIn$theta12^2) * (1 - TIn$theta23^2))) * ((TIn$mar3 - TIn$theta23 * TIn$mar2)/sqrt(1 - TIn$theta23^2)) )/sqrt(1 - ((TIn$theta13 - TIn$theta12 * TIn$theta23)/sqrt((1 - TIn$theta12^2) * (1 - TIn$theta23^2)))^2))
pst.4 <- mm(pst.4, min.pr = VC$min.pr, max.pr = VC$max.pr)
st.4 <- -TIn$theta23/sqrt(1 - TIn$theta23^2)
dst.5 <- dnorm((TIn$mar1[VC$inde] - TIn$theta13 * TIn$mar3)/sqrt(1 - TIn$theta13^2))
pst.5 <- pnorm( ( ((TIn$mar2 - TIn$theta23 * TIn$mar3)/sqrt(1 - TIn$theta23^2)) - ((TIn$theta12 - TIn$theta13 * TIn$theta23)/sqrt((1 - TIn$theta13^2) * (1 - TIn$theta23^2))) * ((TIn$mar1[VC$inde] - TIn$theta13 * TIn$mar3)/sqrt(1 - TIn$theta13^2)) )/sqrt(1 - ((TIn$theta12 - TIn$theta13 * TIn$theta23)/sqrt((1 - TIn$theta13^2) * (1 - TIn$theta23^2)))^2))
pst.5 <- mm(pst.5, min.pr = VC$min.pr, max.pr = VC$max.pr)
st.5 <- -TIn$theta13/sqrt(1 - TIn$theta13^2)
dst.6 <- dnorm((TIn$mar2 - TIn$theta23 * TIn$mar3)/sqrt(1 - TIn$theta23^2))
pst.6 <- pnorm( ( ((TIn$mar1[VC$inde] - TIn$theta13 * TIn$mar3)/sqrt(1 - TIn$theta13^2)) - ((TIn$theta12 - TIn$theta13 * TIn$theta23)/sqrt((1 - TIn$theta13^2) * (1 - TIn$theta23^2))) * ((TIn$mar2 - TIn$theta23 * TIn$mar3)/sqrt(1 - TIn$theta23^2)) )/sqrt(1 - ((TIn$theta12 - TIn$theta13 * TIn$theta23)/sqrt((1 - TIn$theta13^2) * (1 - TIn$theta23^2)))^2))
pst.6 <- mm(pst.6, min.pr = VC$min.pr, max.pr = VC$max.pr)
st.6 <- -TIn$theta23/sqrt(1 - TIn$theta23^2)
dp.1.11.de1 <- dst.1 * pst.1 * st.1 + dst.2 * pst.2 * st.2
dp.1.10.de1 <- dst.1 * (1 - pst.1) * st.1 + dst.2 * pst.2 * ( -st.2 )
dp.1.00.de1 <- dst.1 * (1 - pst.1) * ( -st.1 ) + dst.2 * (1 - pst.2) * ( -st.2 )
dp.1.01.de1 <- dst.1 * pst.1 * ( -st.1 ) + dst.2 * (1 - pst.2) * st.2
dp.2.11.de2 <- dst.3 * pst.3 * st.3 + dst.4 * pst.4 * st.4
dp.2.10.de2 <- dst.3 * (1 - pst.3) * st.3 + dst.4 * pst.4 * ( -st.4 )
dp.3.11.de3 <- dst.5 * pst.5 * st.5 + dst.6 * pst.6 * st.6
dp.3.10.de3 <- dst.5 * (1 - pst.5) * st.5 + dst.6 * pst.6 * ( -st.6 )
der2p.dereta1 <- probm(TIn$eta1, VC$margins[1], only.pr = FALSE, min.dn = VC$min.dn, min.pr = VC$min.pr, max.pr = VC$max.pr)$der2p.dereta
der2p.dereta2 <- probm(TIn$eta2, VC$margins[2], only.pr = FALSE, min.dn = VC$min.dn, min.pr = VC$min.pr, max.pr = VC$max.pr)$der2p.dereta
der2p.dereta3 <- probm(TIn$eta3, VC$margins[3], only.pr = FALSE, min.dn = VC$min.dn, min.pr = VC$min.pr, max.pr = VC$max.pr)$der2p.dereta
d2F1.de1 <- (TIn$mar1 * LgTRI$dmar1^2)/LgTRI$d.1^2 + der2p.dereta1/LgTRI$d.1
d2F2.de2 <- (TIn$mar2 * LgTRI$dmar2^2)/LgTRI$d.2^2 + der2p.dereta2/LgTRI$d.2
d2F3.de3 <- (TIn$mar3 * LgTRI$dmar3^2)/LgTRI$d.3^2 + der2p.dereta3/LgTRI$d.3
d2l.dF1.F1.1 <- respvec$cy1 * ( -1/TIn$p0^2 * (LgTRI$d.1)^2 + 1/TIn$p0 * (TIn$mar1 * LgTRI$d.1) )
d2l.dF1.F1.1[VC$inde] <- respvec$y1.y2.y3 * ( -1/TIn$p111^2 * (LgTRI$d.1[VC$inde] * LgTRI$p.1.11)^2 + 1/TIn$p111 * ( -TIn$mar1[VC$inde] * LgTRI$d.1[VC$inde] * LgTRI$p.1.11 + LgTRI$d.1[VC$inde] * dp.1.11.de1) ) +
respvec$y1.y2.cy3 * ( -1/TIn$p110^2 * (LgTRI$d.1[VC$inde] * LgTRI$p.1.10)^2 + 1/TIn$p110 * ( -TIn$mar1[VC$inde] * LgTRI$d.1[VC$inde] * LgTRI$p.1.10 + LgTRI$d.1[VC$inde] * dp.1.10.de1) ) +
respvec$y1.cy2.cy3 * ( -1/TIn$p100^2 * (LgTRI$d.1[VC$inde] * LgTRI$p.1.00)^2 + 1/TIn$p100 * ( -TIn$mar1[VC$inde] * LgTRI$d.1[VC$inde] * LgTRI$p.1.00 + LgTRI$d.1[VC$inde] * dp.1.00.de1) ) +
respvec$y1.cy2.y3 * ( -1/TIn$p101^2 * (LgTRI$d.1[VC$inde] * LgTRI$p.1.01)^2 + 1/TIn$p101 * ( -TIn$mar1[VC$inde] * LgTRI$d.1[VC$inde] * LgTRI$p.1.01 + LgTRI$d.1[VC$inde] * dp.1.01.de1) )
d2l.dF1.F1 <- d2l.dF1.F1.1
d2l.de1.e1 <- d2l.dF1.F1 * LgTRI$dF1.de1^2 + LgTRI$dl.dF1 * d2F1.de1
d2l.dF2.F2 <- respvec$y1.y2.y3 * ( -1/TIn$p111^2 * (LgTRI$d.2 * LgTRI$p.2.11)^2 + 1/TIn$p111 * ( -TIn$mar2 * LgTRI$d.2 * LgTRI$p.2.11 + LgTRI$d.2 * dp.2.11.de2) ) +
respvec$y1.y2.cy3 * ( -1/TIn$p110^2 * (LgTRI$d.2 * LgTRI$p.2.10)^2 + 1/TIn$p110 * ( -TIn$mar2 * LgTRI$d.2 * LgTRI$p.2.10 + LgTRI$d.2 * dp.2.10.de2) ) -
respvec$y1.cy2.cy3 * ( 1/TIn$p100^2 * (LgTRI$d.2 * LgTRI$p.2.10)^2 + 1/TIn$p100 * ( -TIn$mar2 * LgTRI$d.2 * LgTRI$p.2.10 + LgTRI$d.2 * dp.2.10.de2) ) -
respvec$y1.cy2.y3 * ( 1/TIn$p101^2 * (LgTRI$d.2 * LgTRI$p.2.11)^2 + 1/TIn$p101 * ( -TIn$mar2 * LgTRI$d.2 * LgTRI$p.2.11 + LgTRI$d.2 * dp.2.11.de2) )
d2l.de2.e2 <- d2l.dF2.F2 * LgTRI$dF2.de2^2 + LgTRI$dl.dF2 * d2F2.de2
d2l.dF3.F3 <- respvec$y1.y2.y3 * ( -1/TIn$p111^2 * (LgTRI$d.3 * LgTRI$p.3.11)^2 + 1/TIn$p111 * ( -TIn$mar3 * LgTRI$d.3 * LgTRI$p.3.11 + LgTRI$d.3 * dp.3.11.de3) ) -
respvec$y1.y2.cy3 * ( 1/TIn$p110^2 * (LgTRI$d.3 * LgTRI$p.3.11)^2 + 1/TIn$p110 * ( -TIn$mar3 * LgTRI$d.3 * LgTRI$p.3.11 + LgTRI$d.3 * dp.3.11.de3) ) -
respvec$y1.cy2.cy3 * ( 1/TIn$p100^2 * (LgTRI$d.3 * LgTRI$p.3.10)^2 + 1/TIn$p100 * ( -TIn$mar3 * LgTRI$d.3 * LgTRI$p.3.10 + LgTRI$d.3 * dp.3.10.de3) ) +
respvec$y1.cy2.y3 * ( -1/TIn$p101^2 * (LgTRI$d.3 * LgTRI$p.3.10)^2 + 1/TIn$p101 * ( -TIn$mar3 * LgTRI$d.3 * LgTRI$p.3.10 + LgTRI$d.3 * dp.3.10.de3) )
d2l.de3.e3 <- d2l.dF3.F3 * LgTRI$dF3.de3^2 + LgTRI$dl.dF3 * d2F3.de3
mean11 <- TIn$theta23 * TIn$mar2 + ((TIn$theta13 - TIn$theta12 * TIn$theta23) * (TIn$mar1[VC$inde] - TIn$theta12 * TIn$mar2))/(1 - TIn$theta12^2)
mean22 <- TIn$theta23 * TIn$mar3 + ((TIn$theta12 - TIn$theta13 * TIn$theta23) * (TIn$mar1[VC$inde] - TIn$theta13 * TIn$mar3))/(1 - TIn$theta13^2)
mean33 <- TIn$theta13 * TIn$mar3 + ((TIn$theta12 - TIn$theta13 * TIn$theta23) * (TIn$mar2 - TIn$theta23 * TIn$mar3))/(1 - TIn$theta23^2)
deno <- 1 - TIn$theta12^2 - TIn$theta13^2 - TIn$theta23^2 + 2 * TIn$theta12 * TIn$theta13 * TIn$theta23
sd11 <- sqrt( deno / ( 1 - TIn$theta12^2 ) )
sd22 <- sqrt( deno / ( 1 - TIn$theta13^2 ) )
sd33 <- sqrt( deno / ( 1 - TIn$theta23^2 ) )
p1.1 <- mm( pnorm((TIn$mar3 - mean11)/sd11) , min.pr = VC$min.pr, max.pr = VC$max.pr)
p2.2 <- mm( pnorm((TIn$mar2 - mean22)/sd22) , min.pr = VC$min.pr, max.pr = VC$max.pr)
p3.3 <- mm( pnorm((TIn$mar1[VC$inde] - mean33)/sd33) , min.pr = VC$min.pr, max.pr = VC$max.pr)
p1.1.c <- mm(1 - p1.1, min.pr = VC$min.pr, max.pr = VC$max.pr)
p2.2.c <- mm(1 - p2.2, min.pr = VC$min.pr, max.pr = VC$max.pr)
d.1.1 <- dnorm(TIn$mar1[VC$inde], mean = TIn$theta12 * TIn$mar2, sd = sqrt(1 - TIn$theta12^2))
d.1.2 <- dnorm(TIn$mar1[VC$inde], mean = TIn$theta13 * TIn$mar3, sd = sqrt(1 - TIn$theta13^2))
d.1.3 <- dnorm(TIn$mar2 , mean = TIn$theta23 * TIn$mar3, sd = sqrt(1 - TIn$theta23^2))
d2l.dF1.F2 <- respvec$y1.y2.y3 * ( -1/TIn$p111^2 * LgTRI$d.2 * LgTRI$p.2.11 * LgTRI$d.1[VC$inde] * LgTRI$p.1.11 + 1/TIn$p111 * LgTRI$d.2 * d.1.1 * p1.1 ) +
respvec$y1.y2.cy3 * ( -1/TIn$p110^2 * LgTRI$d.2 * LgTRI$p.2.10 * LgTRI$d.1[VC$inde] * LgTRI$p.1.10 + 1/TIn$p110 * LgTRI$d.2 * d.1.1 * p1.1.c ) -
respvec$y1.cy2.cy3 * ( -1/TIn$p100^2 * LgTRI$d.2 * LgTRI$p.2.10 * LgTRI$d.1[VC$inde] * LgTRI$p.1.00 + 1/TIn$p100 * LgTRI$d.2 * d.1.1 * p1.1.c ) -
respvec$y1.cy2.y3 * ( -1/TIn$p101^2 * LgTRI$d.2 * LgTRI$p.2.11 * LgTRI$d.1[VC$inde] * LgTRI$p.1.01 + 1/TIn$p101 * LgTRI$d.2 * d.1.1 * p1.1 )
d2l.de1.e2 <- d2l.dF1.F2 * LgTRI$dF1.de1[VC$inde] * LgTRI$dF2.de2
d2l.dF1.F3 <- respvec$y1.y2.y3 * ( -1/TIn$p111^2 * LgTRI$d.3 * LgTRI$p.3.11 * LgTRI$d.1[VC$inde] * LgTRI$p.1.11 + 1/TIn$p111 * LgTRI$d.3 * d.1.2 * p2.2 ) -
respvec$y1.y2.cy3 * ( -1/TIn$p110^2 * LgTRI$d.3 * LgTRI$p.3.11 * LgTRI$d.1[VC$inde] * LgTRI$p.1.10 + 1/TIn$p110 * LgTRI$d.3 * d.1.2 * p2.2 ) -
respvec$y1.cy2.cy3 * ( -1/TIn$p100^2 * LgTRI$d.3 * LgTRI$p.3.10 * LgTRI$d.1[VC$inde] * LgTRI$p.1.00 + 1/TIn$p100 * LgTRI$d.3 * d.1.2 * p2.2.c ) +
respvec$y1.cy2.y3 * ( -1/TIn$p101^2 * LgTRI$d.3 * LgTRI$p.3.10 * LgTRI$d.1[VC$inde] * LgTRI$p.1.01 + 1/TIn$p101 * LgTRI$d.3 * d.1.2 * p2.2.c )
d2l.de1.e3 <- d2l.dF1.F3 * LgTRI$dF1.de1[VC$inde] * LgTRI$dF3.de3
d2l.dF2.F3 <- respvec$y1.y2.y3 * ( -1/TIn$p111^2 * LgTRI$d.3 * LgTRI$p.3.11 * LgTRI$d.2 * LgTRI$p.2.11 + 1/TIn$p111 * LgTRI$d.3 * d.1.3 * p3.3 ) -
respvec$y1.y2.cy3 * ( -1/TIn$p110^2 * LgTRI$d.3 * LgTRI$p.3.11 * LgTRI$d.2 * LgTRI$p.2.10 + 1/TIn$p110 * LgTRI$d.3 * d.1.3 * p3.3 ) +
respvec$y1.cy2.cy3 * ( -1/TIn$p100^2 * LgTRI$d.3 * LgTRI$p.3.10 * LgTRI$d.2 * LgTRI$p.2.10 + 1/TIn$p100 * LgTRI$d.3 * d.1.3 * p3.3 ) -
respvec$y1.cy2.y3 * ( -1/TIn$p101^2 * LgTRI$d.3 * LgTRI$p.3.10 * LgTRI$d.2 * LgTRI$p.2.11 + 1/TIn$p101 * LgTRI$d.3 * d.1.3 * p3.3 )
d2l.de2.e3 <- d2l.dF2.F3 * LgTRI$dF2.de2 * LgTRI$dF3.de3
d12 <- dnorm( (TIn$mar3 - LgTRI$mean.12)/LgTRI$sd.12 )
d13 <- dnorm( (TIn$mar2 - LgTRI$mean.13)/LgTRI$sd.13 )
d23 <- dnorm( (TIn$mar1[VC$inde] - LgTRI$mean.23)/LgTRI$sd.23 )
d2l.dF1.theta12 <- respvec$y1.y2.y3 * ( -1/TIn$p111^2 * LgTRI$d.1[VC$inde] * LgTRI$p.1.11 * LgTRI$d11.12 * LgTRI$p12.g + 1/TIn$p111 * ( LgTRI$d11.12 * (TIn$theta12*TIn$mar2 - TIn$mar1[VC$inde])/(1 - TIn$theta12^2) * LgTRI$p12.g - LgTRI$d11.12 * d12/LgTRI$sd.12 * ((TIn$theta13 - TIn$theta12 * TIn$theta23)/(1 - TIn$theta12^2)) )) +
respvec$y1.y2.cy3 * ( -1/TIn$p110^2 * LgTRI$d.1[VC$inde] * LgTRI$p.1.10 * LgTRI$d11.12 * LgTRI$p12.g.c + 1/TIn$p110 * ( LgTRI$d11.12 * (TIn$theta12*TIn$mar2 - TIn$mar1[VC$inde])/(1 - TIn$theta12^2) * LgTRI$p12.g.c + LgTRI$d11.12 * d12/LgTRI$sd.12 * ((TIn$theta13 - TIn$theta12 * TIn$theta23)/(1 - TIn$theta12^2)) )) -
respvec$y1.cy2.cy3 * ( -1/TIn$p100^2 * LgTRI$d.1[VC$inde] * LgTRI$p.1.00 * LgTRI$d11.12 * LgTRI$p12.g.c + 1/TIn$p100 * ( LgTRI$d11.12 * (TIn$theta12*TIn$mar2 - TIn$mar1[VC$inde])/(1 - TIn$theta12^2) * LgTRI$p12.g.c + LgTRI$d11.12 * d12/LgTRI$sd.12 * ((TIn$theta13 - TIn$theta12 * TIn$theta23)/(1 - TIn$theta12^2)) )) -
respvec$y1.cy2.y3 * ( -1/TIn$p101^2 * LgTRI$d.1[VC$inde] * LgTRI$p.1.01 * LgTRI$d11.12 * LgTRI$p12.g + 1/TIn$p101 * ( LgTRI$d11.12 * (TIn$theta12*TIn$mar2 - TIn$mar1[VC$inde])/(1 - TIn$theta12^2) * LgTRI$p12.g - LgTRI$d11.12 * d12/LgTRI$sd.12 * ((TIn$theta13 - TIn$theta12 * TIn$theta23)/(1 - TIn$theta12^2)) ))
d2l.de1.theta12 <- d2l.dF1.theta12 * LgTRI$dF1.de1[VC$inde]
d2l.dF1.theta13 <- respvec$y1.y2.y3 * ( -1/TIn$p111^2 * LgTRI$d.1[VC$inde] * LgTRI$p.1.11 * LgTRI$d11.13 * LgTRI$p13.g + 1/TIn$p111 * ( LgTRI$d11.13 * (TIn$theta13 * TIn$mar3 - TIn$mar1[VC$inde])/(1 - TIn$theta13^2) * LgTRI$p13.g - LgTRI$d11.13 * d13/LgTRI$sd.13 * ((TIn$theta12 - TIn$theta13 * TIn$theta23)/(1 - TIn$theta13^2)) )) -
respvec$y1.y2.cy3 * ( -1/TIn$p110^2 * LgTRI$d.1[VC$inde] * LgTRI$p.1.10 * LgTRI$d11.13 * LgTRI$p13.g + 1/TIn$p110 * ( LgTRI$d11.13 * (TIn$theta13 * TIn$mar3 - TIn$mar1[VC$inde])/(1 - TIn$theta13^2) * LgTRI$p13.g - LgTRI$d11.13 * d13/LgTRI$sd.13 * ((TIn$theta12 - TIn$theta13 * TIn$theta23)/(1 - TIn$theta13^2)) )) -
respvec$y1.cy2.cy3 * ( -1/TIn$p100^2 * LgTRI$d.1[VC$inde] * LgTRI$p.1.00 * LgTRI$d11.13 * LgTRI$p13.g.c + 1/TIn$p100 * ( LgTRI$d11.13 * (TIn$theta13 * TIn$mar3 - TIn$mar1[VC$inde])/(1 - TIn$theta13^2) * LgTRI$p13.g.c + LgTRI$d11.13 * d13/LgTRI$sd.13 * ((TIn$theta12 - TIn$theta13 * TIn$theta23)/(1 - TIn$theta13^2)) )) +
respvec$y1.cy2.y3 * ( -1/TIn$p101^2 * LgTRI$d.1[VC$inde] * LgTRI$p.1.01 * LgTRI$d11.13 * LgTRI$p13.g.c + 1/TIn$p101 * ( LgTRI$d11.13 * (TIn$theta13 * TIn$mar3 - TIn$mar1[VC$inde])/(1 - TIn$theta13^2) * LgTRI$p13.g.c + LgTRI$d11.13 * d13/LgTRI$sd.13 * ((TIn$theta12 - TIn$theta13 * TIn$theta23)/(1 - TIn$theta13^2)) ))
d2l.de1.theta13 <- d2l.dF1.theta13 * LgTRI$dF1.de1[VC$inde]
d2l.dF1.theta23 <- respvec$y1.y2.y3 * ( -1/TIn$p111^2 * LgTRI$d.1[VC$inde] * LgTRI$p.1.11 * LgTRI$d11.23 * LgTRI$p23.g + 1/TIn$p111 * LgTRI$d11.23 * d23/LgTRI$sd.23) -
respvec$y1.y2.cy3 * ( -1/TIn$p110^2 * LgTRI$d.1[VC$inde] * LgTRI$p.1.10 * LgTRI$d11.23 * LgTRI$p23.g + 1/TIn$p110 * LgTRI$d11.23 * d23/LgTRI$sd.23 ) +
respvec$y1.cy2.cy3 * ( -1/TIn$p100^2 * LgTRI$d.1[VC$inde] * LgTRI$p.1.00 * LgTRI$d11.23 * LgTRI$p23.g + 1/TIn$p100 * LgTRI$d11.23 * d23/LgTRI$sd.23 ) -
respvec$y1.cy2.y3 * ( -1/TIn$p101^2 * LgTRI$d.1[VC$inde] * LgTRI$p.1.01 * LgTRI$d11.23 * LgTRI$p23.g + 1/TIn$p101 * LgTRI$d11.23 * d23/LgTRI$sd.23 )
d2l.de1.theta23 <- d2l.dF1.theta23 * LgTRI$dF1.de1[VC$inde]
d2l.dF2.theta12 <- respvec$y1.y2.y3 * ( -1/TIn$p111^2 * LgTRI$d.2 * LgTRI$p.2.11 * LgTRI$d11.12 * LgTRI$p12.g + 1/TIn$p111 * (LgTRI$d11.12 * (TIn$theta12*TIn$mar1[VC$inde] - TIn$mar2)/(1 - TIn$theta12^2) * LgTRI$p12.g - LgTRI$d11.12 * (d12/LgTRI$sd.12) * (TIn$theta23 - TIn$theta12 * TIn$theta13)/(1 - TIn$theta12^2) ) ) +
respvec$y1.y2.cy3 * ( -1/TIn$p110^2 * LgTRI$d.2 * LgTRI$p.2.10 * LgTRI$d11.12 * LgTRI$p12.g.c + 1/TIn$p110 * (LgTRI$d11.12 * (TIn$theta12 * TIn$mar1[VC$inde] - TIn$mar2)/(1 - TIn$theta12^2) * LgTRI$p12.g.c + LgTRI$d11.12 * (d12/LgTRI$sd.12) * (TIn$theta23 - TIn$theta12 * TIn$theta13)/(1 - TIn$theta12^2) ) ) -
respvec$y1.cy2.cy3 * ( 1/TIn$p100^2 * LgTRI$d.2 * LgTRI$p.2.10 * LgTRI$d11.12 * LgTRI$p12.g.c + 1/TIn$p100 * (LgTRI$d11.12 * (TIn$theta12 * TIn$mar1[VC$inde] - TIn$mar2)/(1 - TIn$theta12^2) * LgTRI$p12.g.c + LgTRI$d11.12 * (d12/LgTRI$sd.12) * (TIn$theta23 - TIn$theta12 * TIn$theta13)/(1 - TIn$theta12^2) ) ) -
respvec$y1.cy2.y3 * ( 1/TIn$p101^2 * LgTRI$d.2 * LgTRI$p.2.11 * LgTRI$d11.12 * LgTRI$p12.g + 1/TIn$p101 * (LgTRI$d11.12 * (TIn$theta12 * TIn$mar1[VC$inde] - TIn$mar2)/(1 - TIn$theta12^2) * LgTRI$p12.g - LgTRI$d11.12 * (d12/LgTRI$sd.12) * (TIn$theta23 - TIn$theta12 * TIn$theta13)/(1 - TIn$theta12^2) ) )
d2l.de2.theta12 <- d2l.dF2.theta12 * LgTRI$dF2.de2
d2l.dF2.theta13 <- respvec$y1.y2.y3 * ( -1/TIn$p111^2 * LgTRI$d.2 * LgTRI$p.2.11 * LgTRI$d11.13 * LgTRI$p13.g + 1/TIn$p111 * LgTRI$d11.13 * d13/LgTRI$sd.13 )-
respvec$y1.y2.cy3 * ( -1/TIn$p110^2 * LgTRI$d.2 * LgTRI$p.2.10 * LgTRI$d11.13 * LgTRI$p13.g + 1/TIn$p110 * LgTRI$d11.13 * d13/LgTRI$sd.13 ) +
respvec$y1.cy2.cy3 * ( -1/TIn$p100^2 * LgTRI$d.2 * LgTRI$p.2.10 * LgTRI$d11.13 * LgTRI$p13.g.c + 1/TIn$p100 * LgTRI$d11.13 * d13/LgTRI$sd.13 ) -
respvec$y1.cy2.y3 * ( -1/TIn$p101^2 * LgTRI$d.2 * LgTRI$p.2.11 * LgTRI$d11.13 * LgTRI$p13.g.c + 1/TIn$p101 * LgTRI$d11.13 * d13/LgTRI$sd.13 )
d2l.de2.theta13 <- d2l.dF2.theta13 * LgTRI$dF2.de2
d2l.dF2.theta23 <- respvec$y1.y2.y3 * ( -1/TIn$p111^2 * LgTRI$d.2 * LgTRI$p.2.11 * LgTRI$d11.23 * LgTRI$p23.g + 1/TIn$p111 * (LgTRI$d11.23 * (TIn$theta23 * TIn$mar3 - TIn$mar2)/(1 - TIn$theta23^2) * LgTRI$p23.g - LgTRI$d11.23 * d23/LgTRI$sd.23 * (TIn$theta12 - TIn$theta13 * TIn$theta23)/(1 - TIn$theta23^2) ) ) -
respvec$y1.y2.cy3 * ( -1/TIn$p110^2 * LgTRI$d.2 * LgTRI$p.2.10 * LgTRI$d11.23 * LgTRI$p23.g + 1/TIn$p110 * (LgTRI$d11.23 * (TIn$theta23 * TIn$mar3 - TIn$mar2)/(1 - TIn$theta23^2) * LgTRI$p23.g - LgTRI$d11.23 * d23/LgTRI$sd.23 * (TIn$theta12 - TIn$theta13 * TIn$theta23)/(1 - TIn$theta23^2) ) ) +
respvec$y1.cy2.cy3 * ( 1/TIn$p100^2 * LgTRI$d.2 * LgTRI$p.2.10 * LgTRI$d11.23 * LgTRI$p23.g + 1/TIn$p100 * (LgTRI$d11.23 * (TIn$theta23 * TIn$mar3 - TIn$mar2)/(1 - TIn$theta23^2) * LgTRI$p23.g - LgTRI$d11.23 * d23/LgTRI$sd.23 * (TIn$theta12 - TIn$theta13 * TIn$theta23)/(1 - TIn$theta23^2) ) ) -
respvec$y1.cy2.y3 * ( 1/TIn$p101^2 * LgTRI$d.2 * LgTRI$p.2.11 * LgTRI$d11.23 * LgTRI$p23.g + 1/TIn$p101 * (LgTRI$d11.23 * (TIn$theta23 * TIn$mar3 - TIn$mar2)/(1 - TIn$theta23^2) * LgTRI$p23.g - LgTRI$d11.23 * d23/LgTRI$sd.23 * (TIn$theta12 - TIn$theta13 * TIn$theta23)/(1 - TIn$theta23^2) ) )
d2l.de2.theta23 <- d2l.dF2.theta23 * LgTRI$dF2.de2
d2l.dF3.theta12 <- respvec$y1.y2.y3 * ( -1/TIn$p111^2 * LgTRI$d.3 * LgTRI$p.3.11 * LgTRI$d11.12 * LgTRI$p12.g + 1/TIn$p111 * LgTRI$d11.12 * d12/LgTRI$sd.12 ) -
respvec$y1.y2.cy3 * ( -1/TIn$p110^2 * LgTRI$d.3 * LgTRI$p.3.11 * LgTRI$d11.12 * LgTRI$p12.g.c + 1/TIn$p110 * LgTRI$d11.12 * d12/LgTRI$sd.12 ) +
respvec$y1.cy2.cy3 * ( -1/TIn$p100^2 * LgTRI$d.3 * LgTRI$p.3.10 * LgTRI$d11.12 * LgTRI$p12.g.c + 1/TIn$p100 * LgTRI$d11.12 * d12/LgTRI$sd.12 ) -
respvec$y1.cy2.y3 * ( -1/TIn$p101^2 * LgTRI$d.3 * LgTRI$p.3.10 * LgTRI$d11.12 * LgTRI$p12.g + 1/TIn$p101 * LgTRI$d11.12 * d12/LgTRI$sd.12 )
d2l.de3.theta12 <- d2l.dF3.theta12 * LgTRI$dF3.de3
d2l.dF3.theta13 <- respvec$y1.y2.y3 * ( -1/TIn$p111^2 * LgTRI$d.3 * LgTRI$p.3.11 * LgTRI$d11.13 * LgTRI$p13.g + 1/TIn$p111 * (LgTRI$d11.13 * ((TIn$theta13*TIn$mar1[VC$inde] - TIn$mar3)/(1 - TIn$theta13^2)) * LgTRI$p13.g - LgTRI$d11.13 * (d13/LgTRI$sd.13) * (TIn$theta23 - TIn$theta12 * TIn$theta13)/(1 - TIn$theta13^2) ) ) -
respvec$y1.y2.cy3 * ( 1/TIn$p110^2 * LgTRI$d.3 * LgTRI$p.3.11 * LgTRI$d11.13 * LgTRI$p13.g + 1/TIn$p110 * (LgTRI$d11.13 * ((TIn$theta13 * TIn$mar1[VC$inde] - TIn$mar3)/(1 - TIn$theta13^2)) * LgTRI$p13.g - LgTRI$d11.13 * (d13/LgTRI$sd.13) * (TIn$theta23 - TIn$theta12 * TIn$theta13)/(1 - TIn$theta13^2) ) ) -
respvec$y1.cy2.cy3 * ( 1/TIn$p100^2 * LgTRI$d.3 * LgTRI$p.3.10 * LgTRI$d11.13 * LgTRI$p13.g.c + 1/TIn$p100 * (LgTRI$d11.13 * ((TIn$theta13 * TIn$mar1[VC$inde] - TIn$mar3)/(1 - TIn$theta13^2)) * LgTRI$p13.g.c + LgTRI$d11.13 * (d13/LgTRI$sd.13) * (TIn$theta23 - TIn$theta12 * TIn$theta13)/(1 - TIn$theta13^2) ) ) +
respvec$y1.cy2.y3 * ( -1/TIn$p101^2 * LgTRI$d.3 * LgTRI$p.3.10 * LgTRI$d11.13 * LgTRI$p13.g.c + 1/TIn$p101 * (LgTRI$d11.13 * ((TIn$theta13 * TIn$mar1[VC$inde] - TIn$mar3)/(1 - TIn$theta13^2)) * LgTRI$p13.g.c + LgTRI$d11.13 * (d13/LgTRI$sd.13) * (TIn$theta23 - TIn$theta12 * TIn$theta13)/(1 - TIn$theta13^2) ) )
d2l.de3.theta13 <- d2l.dF3.theta13 * LgTRI$dF3.de3
d2l.dF3.theta23 <- respvec$y1.y2.y3 * ( -1/TIn$p111^2 * LgTRI$d.3 * LgTRI$p.3.11 * LgTRI$d11.23 * LgTRI$p23.g + 1/TIn$p111 * (LgTRI$d11.23 * ((TIn$theta23*TIn$mar2 - TIn$mar3)/(1 - TIn$theta23^2)) * LgTRI$p23.g - LgTRI$d11.23 * d23/LgTRI$sd.23 * (TIn$theta13 - TIn$theta12 * TIn$theta23)/(1 - TIn$theta23^2) ) ) -
respvec$y1.y2.cy3 * ( 1/TIn$p110^2 * LgTRI$d.3 * LgTRI$p.3.11 * LgTRI$d11.23 * LgTRI$p23.g + 1/TIn$p110 * (LgTRI$d11.23 * ((TIn$theta23*TIn$mar2 - TIn$mar3)/(1 - TIn$theta23^2)) * LgTRI$p23.g - LgTRI$d11.23 * d23/LgTRI$sd.23 * (TIn$theta13 - TIn$theta12*TIn$theta23)/(1 - TIn$theta23^2) ) ) +
respvec$y1.cy2.cy3 * ( 1/TIn$p100^2 * LgTRI$d.3 * LgTRI$p.3.10 * LgTRI$d11.23 * LgTRI$p23.g + 1/TIn$p100 * (LgTRI$d11.23 * ((TIn$theta23*TIn$mar2 - TIn$mar3)/(1 - TIn$theta23^2)) * LgTRI$p23.g - LgTRI$d11.23 * d23/LgTRI$sd.23 * (TIn$theta13 - TIn$theta12 * TIn$theta23)/(1 - TIn$theta23^2) ) ) -
respvec$y1.cy2.y3 * ( -1/TIn$p101^2 * LgTRI$d.3 * LgTRI$p.3.10 * LgTRI$d11.23 * LgTRI$p23.g + 1/TIn$p101 * (LgTRI$d11.23 * ((TIn$theta23*TIn$mar2 - TIn$mar3)/(1 - TIn$theta23^2)) * LgTRI$p23.g - LgTRI$d11.23 * d23/LgTRI$sd.23 * (TIn$theta13 - TIn$theta12 * TIn$theta23)/(1 - TIn$theta23^2) ) )
d2l.de3.theta23 <- d2l.dF3.theta23 * LgTRI$dF3.de3
dd11.12.dtheta12 <- (LgTRI$d11.12/(1 - TIn$theta12^2)) * ( ( (TIn$theta12 * TIn$mar2 - TIn$mar1[VC$inde]) * (TIn$theta12 * TIn$mar1[VC$inde] - TIn$mar2)/(1 - TIn$theta12^2) ) + TIn$theta12)
dd11.13.dtheta13 <- (LgTRI$d11.13/(1 - TIn$theta13^2)) * ( ( (TIn$theta13 * TIn$mar3 - TIn$mar1[VC$inde]) * (TIn$theta13 * TIn$mar1[VC$inde] - TIn$mar3)/(1 - TIn$theta13^2) ) + TIn$theta13)
dd11.23.dtheta23 <- (LgTRI$d11.23/(1 - TIn$theta23^2)) * ( ( (TIn$theta23 * TIn$mar3 - TIn$mar2) * (TIn$theta23 * TIn$mar2 - TIn$mar3)/(1 - TIn$theta23^2) ) + TIn$theta23)
dmean12.dtheta12 <- ( (-TIn$mar1[VC$inde] * TIn$theta23 - TIn$mar2 * TIn$theta13) * (1 - TIn$theta12^2) + 2 * TIn$theta12 * (TIn$mar1[VC$inde] * (TIn$theta13 - TIn$theta12 * TIn$theta23) + TIn$mar2 * (TIn$theta23 - TIn$theta12 * TIn$theta13)))/(1 - TIn$theta12^2)^2
dmean13.dtheta13 <- ( (-TIn$mar1[VC$inde] * TIn$theta23 - TIn$mar3 * TIn$theta12) * (1 - TIn$theta13^2) + 2 * TIn$theta13 * (TIn$mar1[VC$inde] * (TIn$theta12 - TIn$theta13 * TIn$theta23) + TIn$mar3 * (TIn$theta23 - TIn$theta12 * TIn$theta13)))/(1 - TIn$theta13^2)^2
dmean23.dtheta23 <- ( (-TIn$mar2 * TIn$theta13 - TIn$mar3 * TIn$theta12) * (1 - TIn$theta23^2) + 2 * TIn$theta23 * (TIn$mar2 * (TIn$theta12 - TIn$theta13 * TIn$theta23) + TIn$mar3 * (TIn$theta13 - TIn$theta12 * TIn$theta23)))/(1 - TIn$theta23^2)^2
dmean13.dtheta12 <- ( TIn$mar1[VC$inde] - TIn$mar3 * TIn$theta13 )/( 1 - TIn$theta13^2 )
dmean23.dtheta12 <- ( TIn$mar2 - TIn$mar3 * TIn$theta23 )/( 1 - TIn$theta23^2 )
dmean23.dtheta13 <- ( TIn$mar3 - TIn$mar2 * TIn$theta23 )/( 1 - TIn$theta23^2 )
dvar12.dtheta12 <- ( 2 * (TIn$theta13 * TIn$theta23 - TIn$theta12) * (1 - TIn$theta12^2) + 2 * TIn$theta12 * deno )/(1 - TIn$theta12^2)^2
dvar13.dtheta13 <- ( 2 * (TIn$theta12 * TIn$theta23 - TIn$theta13) * (1 - TIn$theta13^2) + 2 * TIn$theta13 * deno )/(1 - TIn$theta13^2)^2
dvar23.dtheta23 <- ( 2 * (TIn$theta12 * TIn$theta13 - TIn$theta23) * (1 - TIn$theta23^2) + 2 * TIn$theta23 * deno )/(1 - TIn$theta23^2)^2
dvar13.dtheta12 <- ( 2 * (TIn$theta13 * TIn$theta23 - TIn$theta12 ) )/(1 - TIn$theta13^2)
dvar23.dtheta12 <- ( 2 * (TIn$theta13 * TIn$theta23 - TIn$theta12 ) )/(1 - TIn$theta23^2)
dvar23.dtheta13 <- ( 2 * (TIn$theta12 * TIn$theta23 - TIn$theta13 ) )/(1 - TIn$theta23^2)
d2l.dtheta12.theta12 <- respvec$y1.y2.y3 * ( -1/TIn$p111^2 * (LgTRI$d11.12 * LgTRI$p12.g )^2 + 1/TIn$p111 * (( dd11.12.dtheta12 * LgTRI$p12.g ) - ((d12/LgTRI$sd.12) * (dmean12.dtheta12 + (((TIn$mar3 - LgTRI$mean.12)/(2 * LgTRI$sd.12^2)) * dvar12.dtheta12)) * LgTRI$d11.12) )) +
respvec$y1.y2.cy3 * ( -1/TIn$p110^2 * (LgTRI$d11.12 * LgTRI$p12.g.c )^2 + 1/TIn$p110 * (( dd11.12.dtheta12 * LgTRI$p12.g.c) + ((d12/LgTRI$sd.12) * (dmean12.dtheta12 + (((TIn$mar3 - LgTRI$mean.12)/(2 * LgTRI$sd.12^2)) * dvar12.dtheta12)) * LgTRI$d11.12) )) -
respvec$y1.cy2.cy3 * ( 1/TIn$p100^2 * (LgTRI$d11.12 * LgTRI$p12.g.c )^2 + 1/TIn$p100 * (( dd11.12.dtheta12 * LgTRI$p12.g.c) + ((d12/LgTRI$sd.12) * (dmean12.dtheta12 + (((TIn$mar3 - LgTRI$mean.12)/(2 * LgTRI$sd.12^2)) * dvar12.dtheta12)) * LgTRI$d11.12) )) -
respvec$y1.cy2.y3 * ( 1/TIn$p101^2 * (LgTRI$d11.12 * LgTRI$p12.g )^2 + 1/TIn$p101 * (( dd11.12.dtheta12 * LgTRI$p12.g ) - ((d12/LgTRI$sd.12) * (dmean12.dtheta12 + (((TIn$mar3 - LgTRI$mean.12)/(2 * LgTRI$sd.12^2)) * dvar12.dtheta12)) * LgTRI$d11.12) ))
d2l.dtheta13.theta13 <- respvec$y1.y2.y3 * ( -1/TIn$p111^2 * (LgTRI$d11.13 * LgTRI$p13.g)^2 + 1/TIn$p111 * (( dd11.13.dtheta13 * LgTRI$p13.g ) - ((d13/LgTRI$sd.13) * (dmean13.dtheta13 + (((TIn$mar2 - LgTRI$mean.13)/(2 * LgTRI$sd.13^2)) * dvar13.dtheta13)) * LgTRI$d11.13) )) -
respvec$y1.y2.cy3 * ( 1/TIn$p110^2 * (LgTRI$d11.13 * LgTRI$p13.g)^2 + 1/TIn$p110 * ((dd11.13.dtheta13 * LgTRI$p13.g ) - ((d13/LgTRI$sd.13) * (dmean13.dtheta13 + (((TIn$mar2 - LgTRI$mean.13)/(2 * LgTRI$sd.13^2)) * dvar13.dtheta13)) * LgTRI$d11.13) )) -
respvec$y1.cy2.cy3 * ( 1/TIn$p100^2 * (LgTRI$d11.13 * LgTRI$p13.g.c )^2 + 1/TIn$p100 * ((dd11.13.dtheta13 * LgTRI$p13.g.c) + ((d13/LgTRI$sd.13) * (dmean13.dtheta13 + (((TIn$mar2 - LgTRI$mean.13)/(2 * LgTRI$sd.13^2)) * dvar13.dtheta13)) * LgTRI$d11.13) )) +
respvec$y1.cy2.y3 * ( -1/TIn$p101^2 * (LgTRI$d11.13 * LgTRI$p13.g.c )^2 + 1/TIn$p101 * ((dd11.13.dtheta13 * LgTRI$p13.g.c) + ((d13/LgTRI$sd.13) * (dmean13.dtheta13 + (((TIn$mar2 - LgTRI$mean.13)/(2 * LgTRI$sd.13^2)) * dvar13.dtheta13)) * LgTRI$d11.13) ))
d2l.dtheta23.theta23 <- respvec$y1.y2.y3 * ( -1/TIn$p111^2 * (LgTRI$d11.23 * LgTRI$p23.g )^2 + 1/TIn$p111 * ((dd11.23.dtheta23 * LgTRI$p23.g ) - ((d23/LgTRI$sd.23) * (dmean23.dtheta23 + ((TIn$mar1[VC$inde] - LgTRI$mean.23)/(2*LgTRI$sd.23^2) * dvar23.dtheta23)) * LgTRI$d11.23) )) -
respvec$y1.y2.cy3 * ( 1/TIn$p110^2 * (LgTRI$d11.23 * LgTRI$p23.g )^2 + 1/TIn$p110 * ((dd11.23.dtheta23 * LgTRI$p23.g ) - ((d23/LgTRI$sd.23) * (dmean23.dtheta23 + ((TIn$mar1[VC$inde] - LgTRI$mean.23)/(2*LgTRI$sd.23^2) * dvar23.dtheta23)) * LgTRI$d11.23) )) +
respvec$y1.cy2.cy3 * ( -1/TIn$p100^2 * (LgTRI$d11.23 * LgTRI$p23.g )^2 + 1/TIn$p100 * ((dd11.23.dtheta23 * LgTRI$p23.g ) - ((d23/LgTRI$sd.23) * (dmean23.dtheta23 + ((TIn$mar1[VC$inde] - LgTRI$mean.23)/(2*LgTRI$sd.23^2) * dvar23.dtheta23)) * LgTRI$d11.23) )) -
respvec$y1.cy2.y3 * ( 1/TIn$p101^2 * (LgTRI$d11.23 * LgTRI$p23.g )^2 + 1/TIn$p101 * ((dd11.23.dtheta23 * LgTRI$p23.g ) - ((d23/LgTRI$sd.23) * (dmean23.dtheta23 + ((TIn$mar1[VC$inde] - LgTRI$mean.23)/(2*LgTRI$sd.23^2) * dvar23.dtheta23)) * LgTRI$d11.23) ))
d2l.dtheta12.theta13 <- respvec$y1.y2.y3 * ( -1/TIn$p111^2 * LgTRI$d11.12 * LgTRI$p12.g * LgTRI$d11.13 * LgTRI$p13.g - 1/TIn$p111 * LgTRI$d11.13 * (d13/LgTRI$sd.13) * (dmean13.dtheta12 + ((TIn$mar2-LgTRI$mean.13)/(2*LgTRI$sd.13^2) * dvar13.dtheta12) )) -
respvec$y1.y2.cy3 * ( -1/TIn$p110^2 * LgTRI$d11.12 * LgTRI$p12.g.c * LgTRI$d11.13 * LgTRI$p13.g - 1/TIn$p110 * LgTRI$d11.13 * (d13/LgTRI$sd.13) * (dmean13.dtheta12 + ((TIn$mar2-LgTRI$mean.13)/(2*LgTRI$sd.13^2) * dvar13.dtheta12) )) -
respvec$y1.cy2.cy3 * ( 1/TIn$p100^2 * LgTRI$d11.12 * LgTRI$p12.g.c * LgTRI$d11.13 * LgTRI$p13.g.c + 1/TIn$p100 * LgTRI$d11.13 * (d13/LgTRI$sd.13) * (dmean13.dtheta12 + ((TIn$mar2-LgTRI$mean.13)/(2*LgTRI$sd.13^2) * dvar13.dtheta12) )) +
respvec$y1.cy2.y3 * ( 1/TIn$p101^2 * LgTRI$d11.12 * LgTRI$p12.g * LgTRI$d11.13 * LgTRI$p13.g.c + 1/TIn$p101 * LgTRI$d11.13 * (d13/LgTRI$sd.13) * (dmean13.dtheta12 + ((TIn$mar2-LgTRI$mean.13)/(2*LgTRI$sd.13^2) * dvar13.dtheta12) ))
d2l.dtheta12.theta23 <- respvec$y1.y2.y3 * ( -1/TIn$p111^2 * LgTRI$d11.12 * LgTRI$p12.g * LgTRI$d11.23 * LgTRI$p23.g - 1/TIn$p111 * LgTRI$d11.23 * (d23/LgTRI$sd.23) * (dmean23.dtheta12 + ((TIn$mar1[VC$inde] - LgTRI$mean.23)/(2*LgTRI$sd.23^2) * dvar23.dtheta12) )) -
respvec$y1.y2.cy3 * ( -1/TIn$p110^2 * LgTRI$d11.12 * LgTRI$p12.g.c * LgTRI$d11.23 * LgTRI$p23.g - 1/TIn$p110 * LgTRI$d11.23 * (d23/LgTRI$sd.23) * (dmean23.dtheta12 + ((TIn$mar1[VC$inde] - LgTRI$mean.23)/(2*LgTRI$sd.23^2) * dvar23.dtheta12) )) +
respvec$y1.cy2.cy3 * ( 1/TIn$p100^2 * LgTRI$d11.12 * LgTRI$p12.g.c * LgTRI$d11.23 * LgTRI$p23.g - 1/TIn$p100 * LgTRI$d11.23 * (d23/LgTRI$sd.23) * (dmean23.dtheta12 + ((TIn$mar1[VC$inde] - LgTRI$mean.23)/(2*LgTRI$sd.23^2) * dvar23.dtheta12) )) -
respvec$y1.cy2.y3 * ( 1/TIn$p101^2 * LgTRI$d11.12 * LgTRI$p12.g * LgTRI$d11.23 * LgTRI$p23.g - 1/TIn$p101 * LgTRI$d11.23 * (d23/LgTRI$sd.23) * (dmean23.dtheta12 + ((TIn$mar1[VC$inde] - LgTRI$mean.23)/(2*LgTRI$sd.23^2) * dvar23.dtheta12) ))
d2l.dtheta13.theta23 <- respvec$y1.y2.y3 * ( -1/TIn$p111^2 * LgTRI$d11.13 * LgTRI$p13.g * LgTRI$d11.23 * LgTRI$p23.g - 1/TIn$p111 * LgTRI$d11.23 * (d23/LgTRI$sd.23) * (dmean23.dtheta13 + ((TIn$mar1[VC$inde]-LgTRI$mean.23)/(2 * LgTRI$sd.23^2) * dvar23.dtheta13) )) -
respvec$y1.y2.cy3 * ( 1/TIn$p110^2 * LgTRI$d11.13 * LgTRI$p13.g * LgTRI$d11.23 * LgTRI$p23.g - 1/TIn$p110 * LgTRI$d11.23 * (d23/LgTRI$sd.23) * (dmean23.dtheta13 + ((TIn$mar1[VC$inde]-LgTRI$mean.23)/(2 * LgTRI$sd.23^2) * dvar23.dtheta13) )) +
respvec$y1.cy2.cy3 * ( 1/TIn$p100^2 * LgTRI$d11.13 * LgTRI$p13.g.c * LgTRI$d11.23 * LgTRI$p23.g - 1/TIn$p100 * LgTRI$d11.23 * (d23/LgTRI$sd.23) * (dmean23.dtheta13 + ((TIn$mar1[VC$inde]-LgTRI$mean.23)/(2 * LgTRI$sd.23^2) * dvar23.dtheta13) )) -
respvec$y1.cy2.y3 * ( -1/TIn$p101^2 * LgTRI$d11.13 * LgTRI$p13.g.c * LgTRI$d11.23 * LgTRI$p23.g - 1/TIn$p101 * LgTRI$d11.23 * (d23/LgTRI$sd.23) * (dmean23.dtheta13 + ((TIn$mar1[VC$inde]-LgTRI$mean.23)/(2 * LgTRI$sd.23^2) * dvar23.dtheta13) ))
d2l.dbe1.be1 <- crossprod( VC$X1 * c(VC$weights*d2l.de1.e1), VC$X1 )
d2l.dbe2.be2 <- crossprod( VC$X2 * c(VC$weights[VC$inde]*d2l.de2.e2), VC$X2 )
d2l.dbe3.be3 <- crossprod( VC$X3 * c(VC$weights[VC$inde]*d2l.de3.e3), VC$X3 )
d2l.dbe1.be2 <- crossprod( VC$X1[VC$inde, ] * c(VC$weights[VC$inde]*d2l.de1.e2), VC$X2 )
d2l.dbe1.be3 <- crossprod( VC$X1[VC$inde, ] * c(VC$weights[VC$inde]*d2l.de1.e3), VC$X3 )
d2l.dbe2.be3 <- crossprod( VC$X2 * c(VC$weights[VC$inde]*d2l.de2.e3), VC$X3 )
if(VC$Chol == FALSE){
dtheta12.dtheta12.st <- 4 * exp( 2 * TIn$theta12.st )/( exp(2 * TIn$theta12.st) + 1 )^2
dtheta13.dtheta13.st <- 4 * exp( 2 * TIn$theta13.st )/( exp(2 * TIn$theta13.st) + 1 )^2
dtheta23.dtheta23.st <- 4 * exp( 2 * TIn$theta23.st )/( exp(2 * TIn$theta23.st) + 1 )^2
d2theta12.theta12.st <- (8 * exp( 2 * TIn$theta12.st ) - 8 * exp( 4 * TIn$theta12.st ))/( exp(2 * TIn$theta12.st) + 1 )^3
d2theta13.theta13.st <- (8 * exp( 2 * TIn$theta13.st ) - 8 * exp( 4 * TIn$theta13.st ))/( exp(2 * TIn$theta13.st) + 1 )^3
d2theta23.theta23.st <- (8 * exp( 2 * TIn$theta23.st ) - 8 * exp( 4 * TIn$theta23.st ))/( exp(2 * TIn$theta23.st) + 1 )^3
d2l.dbe1.theta12.st <- colSums(c(VC$weights[VC$inde]*d2l.de1.theta12* dtheta12.dtheta12.st) * VC$X1[VC$inde, ])
d2l.dbe1.theta13.st <- colSums(c(VC$weights[VC$inde]*d2l.de1.theta13* dtheta13.dtheta13.st) * VC$X1[VC$inde, ])
d2l.dbe1.theta23.st <- colSums(c(VC$weights[VC$inde]*d2l.de1.theta23* dtheta23.dtheta23.st) * VC$X1[VC$inde, ])
d2l.dbe2.theta12.st <- colSums(c(VC$weights[VC$inde]*d2l.de2.theta12* dtheta12.dtheta12.st) * VC$X2)
d2l.dbe2.theta13.st <- colSums(c(VC$weights[VC$inde]*d2l.de2.theta13* dtheta13.dtheta13.st) * VC$X2)
d2l.dbe2.theta23.st <- colSums(c(VC$weights[VC$inde]*d2l.de2.theta23* dtheta23.dtheta23.st) * VC$X2)
d2l.dbe3.theta12.st <- colSums(c(VC$weights[VC$inde]*d2l.de3.theta12* dtheta12.dtheta12.st) * VC$X3)
d2l.dbe3.theta13.st <- colSums(c(VC$weights[VC$inde]*d2l.de3.theta13* dtheta13.dtheta13.st) * VC$X3)
d2l.dbe3.theta23.st <- colSums(c(VC$weights[VC$inde]*d2l.de3.theta23* dtheta23.dtheta23.st) * VC$X3)
d2l.dtheta12.st <- sum( VC$weights[VC$inde]*(d2l.dtheta12.theta12 * (dtheta12.dtheta12.st)^2 + LgTRI$dl.dtheta12 * d2theta12.theta12.st) )
d2l.dtheta13.st <- sum( VC$weights[VC$inde]*(d2l.dtheta13.theta13 * (dtheta13.dtheta13.st)^2 + LgTRI$dl.dtheta13 * d2theta13.theta13.st) )
d2l.dtheta23.st <- sum( VC$weights[VC$inde]*(d2l.dtheta23.theta23 * (dtheta23.dtheta23.st)^2 + LgTRI$dl.dtheta23 * d2theta23.theta23.st) )
d2l.dtheta12.st.theta13.st <- sum(VC$weights[VC$inde]*d2l.dtheta12.theta13 * dtheta12.dtheta12.st * dtheta13.dtheta13.st )
d2l.dtheta12.st.theta23.st <- sum(VC$weights[VC$inde]*d2l.dtheta12.theta23 * dtheta12.dtheta12.st * dtheta23.dtheta23.st )
d2l.dtheta13.st.theta23.st <- sum(VC$weights[VC$inde]*d2l.dtheta13.theta23 * dtheta13.dtheta13.st * dtheta23.dtheta23.st )
}
if(VC$Chol == TRUE){
d2l.de1.theta <- matrix(c ( VC$weights[VC$inde]*d2l.de1.theta12,
VC$weights[VC$inde]*d2l.de1.theta13,
VC$weights[VC$inde]*d2l.de1.theta23), length(which(VC$inde==TRUE)), 3)
d2l.de2.theta <- matrix(c ( VC$weights[VC$inde]*d2l.de2.theta12,
VC$weights[VC$inde]*d2l.de2.theta13,
VC$weights[VC$inde]*d2l.de2.theta23), length(which(VC$inde==TRUE)), 3)
d2l.de3.theta <- matrix(c ( VC$weights[VC$inde]*d2l.de3.theta12,
VC$weights[VC$inde]*d2l.de3.theta13,
VC$weights[VC$inde]*d2l.de3.theta23), length(which(VC$inde==TRUE)), 3)
dth12.dth12.st <- 1/(1 + TIn$theta12.st^2)^(3/2)
dth12.dth13.st <- 0
dth12.dth23.st <- 0
dth13.dth12.st <- 0
dth13.dth13.st <- (1 + TIn$theta23.st^2)/(1 + TIn$theta13.st^2 + TIn$theta23.st^2)^(3/2)
dth13.dth23.st <- - (TIn$theta13.st * TIn$theta23.st)/(1 + TIn$theta13.st^2 + TIn$theta23.st^2)^(3/2)
dth23.dth12.st <- TIn$theta13.st/sqrt( (1 + TIn$theta12.st^2) * (1 + TIn$theta13.st^2 + TIn$theta23.st^2)) - (TIn$theta12.st * (TIn$theta12.st * TIn$theta13.st + TIn$theta23.st))/((1 + TIn$theta12.st^2)^(3/2) * sqrt(1 + TIn$theta13.st^2 + TIn$theta23.st^2))
dth23.dth13.st <- TIn$theta12.st/sqrt((1 + TIn$theta12.st^2) * (1 + TIn$theta13.st^2 + TIn$theta23.st^2)) - (TIn$theta13.st * (TIn$theta12.st * TIn$theta13.st + TIn$theta23.st))/(sqrt(1 + TIn$theta12.st^2) * (1 + TIn$theta13.st^2 + TIn$theta23.st^2)^(3/2))
dth23.dth23.st <- 1/sqrt((1 + TIn$theta12.st^2) * (1 + TIn$theta13.st^2 + TIn$theta23.st^2)) - (TIn$theta23.st * (TIn$theta12.st * TIn$theta13.st + TIn$theta23.st))/(sqrt(1 + TIn$theta12.st^2) * (1 + TIn$theta13.st^2 + TIn$theta23.st^2)^(3/2))
dtheta.theta.st <- matrix( c( dth12.dth12.st, dth13.dth12.st, dth23.dth12.st,
dth12.dth13.st, dth13.dth13.st, dth23.dth13.st,
dth12.dth23.st, dth13.dth23.st, dth23.dth23.st ), 3 , 3)
d2l.dbe1.theta.st <- t(VC$X1[VC$inde, ]) %*% d2l.de1.theta %*% dtheta.theta.st
d2l.dbe2.theta.st <- t(VC$X2) %*% d2l.de2.theta %*% dtheta.theta.st
d2l.dbe3.theta.st <- t(VC$X3) %*% d2l.de3.theta %*% dtheta.theta.st
d2l.dbe1.theta12.st <- d2l.dbe1.theta.st[, 1]
d2l.dbe1.theta13.st <- d2l.dbe1.theta.st[, 2]
d2l.dbe1.theta23.st <- d2l.dbe1.theta.st[, 3]
d2l.dbe2.theta12.st <- d2l.dbe2.theta.st[, 1]
d2l.dbe2.theta13.st <- d2l.dbe2.theta.st[, 2]
d2l.dbe2.theta23.st <- d2l.dbe2.theta.st[, 3]
d2l.dbe3.theta12.st <- d2l.dbe3.theta.st[, 1]
d2l.dbe3.theta13.st <- d2l.dbe3.theta.st[, 2]
d2l.dbe3.theta23.st <- d2l.dbe3.theta.st[, 3]
############################
## d2l.dtheta.st.theta.st ##
############################
d2th12.dth12.st.th12.st <- -(3 * TIn$theta12.st)/(1 + TIn$theta12.st^2)^(5/2)
d2th12.dth13.st.th12.st <- 0
d2th12.dth23.st.th12.st <- 0
d2th12.dth12.st.th13.st <- 0
d2th12.dth13.st.th13.st <- 0
d2th12.dth23.st.th13.st <- 0
d2th12.dth12.st.th23.st <- 0
d2th12.dth13.st.th23.st <- 0
d2th12.dth23.st.th23.st <- 0
d2th13.dth12.st.th12.st <- 0
d2th13.dth13.st.th12.st <- 0
d2th13.dth23.st.th12.st <- 0
d2th13.dth12.st.th13.st <- 0
d2th13.dth13.st.th13.st<- -(3 * TIn$theta13.st * (1 + TIn$theta23.st^2))/(1 + TIn$theta13.st^2 + TIn$theta23.st^2)^(5/2)
d2th13.dth23.st.th13.st <- (2 * TIn$theta23.st)/(1 + TIn$theta13.st^2 + TIn$theta23.st^2)^(3/2) - (3 * TIn$theta23.st * (1 + TIn$theta23.st^2))/(1 + TIn$theta13.st^2 + TIn$theta23.st^2)^(5/2)
d2th13.dth12.st.th23.st <- 0
d2th13.dth13.st.th23.st <- -TIn$theta23.st/(1 + TIn$theta13.st^2 + TIn$theta23.st^2)^(3/2) + (3 * TIn$theta13.st^2 * TIn$theta23.st)/(1 + TIn$theta13.st^2 + TIn$theta23.st^2)^(5/2)
d2th13.dth23.st.th23.st <- -TIn$theta13.st/(1 + TIn$theta13.st^2 + TIn$theta23.st^2)^(3/2) + (3 * TIn$theta13.st * TIn$theta23.st^2)/(1 + TIn$theta13.st^2 + TIn$theta23.st^2)^(5/2)
d2th23.dth12.st.th12.st <- -(3 * TIn$theta12.st * TIn$theta13.st + TIn$theta23.st)/((1 + TIn$theta12.st^2)^(3/2) * sqrt(1 + TIn$theta13.st^2 + TIn$theta23.st^2)) + (3 * TIn$theta12.st^2 * (TIn$theta12.st * TIn$theta13.st + TIn$theta23.st))/((1 + TIn$theta12.st^2)^(5/2) * sqrt(1 + TIn$theta13.st^2 + TIn$theta23.st^2))
d2th23.dth13.st.th12.st <- (1 + TIn$theta12.st^2 + TIn$theta23.st^2 + TIn$theta12.st^2 * TIn$theta23.st^2 + TIn$theta12.st^2 * TIn$theta13.st^2 + TIn$theta12.st * TIn$theta13.st * TIn$theta23.st)/((1 + TIn$theta12.st^2) * (1 + TIn$theta13.st^2 + TIn$theta23.st^2))^(3/2) - TIn$theta12.st^2/((1 + TIn$theta12.st^2)^(3/2) * sqrt(1 + TIn$theta13.st^2 + TIn$theta23.st^2))
d2th23.dth23.st.th12.st <- -(TIn$theta13.st * TIn$theta23.st)/(sqrt(1 + TIn$theta12.st^2) * (1 + TIn$theta13.st^2 + TIn$theta23.st^2)^(3/2)) - TIn$theta12.st/((1 + TIn$theta12.st^2)^(3/2) * sqrt(1 + TIn$theta13.st^2 + TIn$theta23.st^2)) + (TIn$theta12.st * TIn$theta23.st * (TIn$theta12.st * TIn$theta13.st + TIn$theta23.st))/((1 + TIn$theta12.st^2) * (1 + TIn$theta13.st^2 + TIn$theta23.st^2))^(3/2)
d2th23.dth12.st.th13.st <- (1 + TIn$theta13.st^2 + TIn$theta23.st^2 + TIn$theta12.st * TIn$theta13.st * (TIn$theta12.st * TIn$theta13.st + TIn$theta23.st))/((1 + TIn$theta12.st^2) * (1 + TIn$theta13.st^2 + TIn$theta23.st^2))^(3/2) - TIn$theta13.st^2/(sqrt(1 + TIn$theta12.st^2) * (1 + TIn$theta13.st^2 + TIn$theta23.st^2)^(3/2))
d2th23.dth13.st.th13.st <- -(3 * TIn$theta12.st * TIn$theta13.st + TIn$theta23.st)/(sqrt(1 + TIn$theta12.st^2) * (1 + TIn$theta13.st^2 + TIn$theta23.st^2)^(3/2)) + (3 * TIn$theta13.st^2 * (TIn$theta12.st * TIn$theta13.st + TIn$theta23.st))/(sqrt(1 + TIn$theta12.st^2) * (1 + TIn$theta13.st^2 + TIn$theta23.st^2)^(5/2))
d2th23.dth23.st.th13.st <- -(TIn$theta12.st * TIn$theta23.st + TIn$theta13.st)/(sqrt(1 + TIn$theta12.st^2) * (1 + TIn$theta13.st^2 + TIn$theta23.st^2)^(3/2)) + (3 * TIn$theta13.st * TIn$theta23.st * (TIn$theta12.st * TIn$theta13.st + TIn$theta23.st))/(sqrt(1 + TIn$theta12.st^2) * (1 + TIn$theta13.st^2 + TIn$theta23.st^2)^(5/2))
d2th23.dth12.st.th23.st <- -TIn$theta12.st/((1 + TIn$theta12.st^2)^(3/2) * sqrt(1 + TIn$theta13.st^2 + TIn$theta23.st^2)) - (TIn$theta13.st * TIn$theta23.st)/(sqrt(1 + TIn$theta12.st^2) * (1 + TIn$theta13.st^2 + TIn$theta23.st^2)^(3/2)) + (TIn$theta12.st * TIn$theta23.st * (TIn$theta12.st * TIn$theta13.st + TIn$theta23.st))/((1 + TIn$theta12.st^2) * (1 + TIn$theta13.st^2 + TIn$theta23.st^2))^(3/2)
d2th23.dth13.st.th23.st <- -(TIn$theta13.st + TIn$theta12.st * TIn$theta23.st)/(sqrt(1 + TIn$theta12.st^2) * (1 + TIn$theta13.st^2 + TIn$theta23.st^2)^(3/2)) + (3 * TIn$theta13.st * TIn$theta23.st * (TIn$theta12.st * TIn$theta13.st + TIn$theta23.st))/(sqrt(1 + TIn$theta12.st^2) * (1 + TIn$theta13.st^2 + TIn$theta23.st^2)^(5/2))
d2th23.dth23.st.th23.st <- -(3 * TIn$theta23.st + TIn$theta12.st * TIn$theta13.st)/(sqrt(1 + TIn$theta12.st^2) * (1 + TIn$theta13.st^2 + TIn$theta23.st^2)^(3/2)) + (3 * TIn$theta23.st^2 * (TIn$theta12.st * TIn$theta13.st + TIn$theta23.st))/(sqrt(1 + TIn$theta12.st^2) * (1 + TIn$theta13.st^2 + TIn$theta23.st^2)^(5/2))
d2theta12.theta12.st <- matrix( c( d2th12.dth12.st.th12.st, d2th13.dth12.st.th12.st, d2th23.dth12.st.th12.st,
d2th12.dth12.st.th13.st, d2th13.dth12.st.th13.st, d2th23.dth12.st.th13.st,
d2th12.dth12.st.th23.st, d2th13.dth12.st.th23.st, d2th23.dth12.st.th23.st ), 3 , 3)
d2theta13.theta13.st <- matrix( c( d2th12.dth13.st.th12.st, d2th13.dth13.st.th12.st, d2th23.dth13.st.th12.st,
d2th12.dth13.st.th13.st, d2th13.dth13.st.th13.st, d2th23.dth13.st.th13.st,
d2th12.dth13.st.th23.st, d2th13.dth13.st.th23.st, d2th23.dth13.st.th23.st ), 3 , 3)
d2theta23.theta23.st <- matrix( c( d2th12.dth23.st.th12.st, d2th13.dth23.st.th12.st, d2th23.dth23.st.th12.st,
d2th12.dth23.st.th13.st, d2th13.dth23.st.th13.st, d2th23.dth23.st.th13.st,
d2th12.dth23.st.th23.st, d2th13.dth23.st.th23.st, d2th23.dth23.st.th23.st ), 3 , 3)
d2l.theta <- matrix( c( sum(VC$weights[VC$inde]*d2l.dtheta12.theta12), sum(VC$weights[VC$inde]*d2l.dtheta12.theta13), sum(VC$weights[VC$inde]*d2l.dtheta12.theta23),
sum(VC$weights[VC$inde]*d2l.dtheta12.theta13), sum(VC$weights[VC$inde]*d2l.dtheta13.theta13), sum(VC$weights[VC$inde]*d2l.dtheta13.theta23),
sum(VC$weights[VC$inde]*d2l.dtheta12.theta23), sum(VC$weights[VC$inde]*d2l.dtheta13.theta23), sum(VC$weights[VC$inde]*d2l.dtheta23.theta23) ), 3 , 3)
dl.dtheta <- matrix(c ( LgTRI$dl.dtheta12,
LgTRI$dl.dtheta13,
LgTRI$dl.dtheta23), length(which(VC$inde==TRUE)), 3)
d2l.dtheta.st <- t(dtheta.theta.st) %*% d2l.theta %*% dtheta.theta.st + matrix(c(colSums(dl.dtheta) %*% d2theta12.theta12.st,
colSums(dl.dtheta) %*% d2theta13.theta13.st,
colSums(dl.dtheta) %*% d2theta23.theta23.st), 3, 3)
d2l.dtheta12.st <- d2l.dtheta.st[1, 1]
d2l.dtheta13.st <- d2l.dtheta.st[2, 2]
d2l.dtheta23.st <- d2l.dtheta.st[3, 3]
d2l.dtheta12.st.theta13.st <- d2l.dtheta.st[1, 2]
d2l.dtheta12.st.theta23.st <- d2l.dtheta.st[1, 3]
d2l.dtheta13.st.theta23.st <- d2l.dtheta.st[2, 3]
}
h1 <- cbind(d2l.dbe1.be1, d2l.dbe1.be2, d2l.dbe1.be3, d2l.dbe1.theta12.st, d2l.dbe1.theta13.st, d2l.dbe1.theta23.st)
h2 <- cbind(t(d2l.dbe1.be2), d2l.dbe2.be2, d2l.dbe2.be3, d2l.dbe2.theta12.st, d2l.dbe2.theta13.st, d2l.dbe2.theta23.st)
h3 <- cbind(t(d2l.dbe1.be3), t(d2l.dbe2.be3), d2l.dbe3.be3, d2l.dbe3.theta12.st, d2l.dbe3.theta13.st, d2l.dbe3.theta23.st)
h4 <- cbind(t(d2l.dbe1.theta12.st), t(d2l.dbe2.theta12.st), t(d2l.dbe3.theta12.st), d2l.dtheta12.st, d2l.dtheta12.st.theta13.st, d2l.dtheta12.st.theta23.st)
h5 <- cbind(t(d2l.dbe1.theta13.st), t(d2l.dbe2.theta13.st), t(d2l.dbe3.theta13.st), t(d2l.dtheta12.st.theta13.st), d2l.dtheta13.st, d2l.dtheta13.st.theta23.st)
h6 <- cbind(t(d2l.dbe1.theta23.st), t(d2l.dbe2.theta23.st), t(d2l.dbe3.theta23.st), t(d2l.dtheta12.st.theta23.st), t(d2l.dtheta13.st.theta23.st), d2l.dtheta23.st)
HTRIVec <- list(h1 = h1, h2 = h2, h3 = h3, h4 = h4, h5 = h5, h6 = h6)
HTRIVec
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.