R/hessian_matrix.R

In BriceonWiley/IntegratedLikelihood.R:

hessian_matrix <- function(params, data) {
  phi <- params[1]; lam2 <- params[2]; theta1 <- params[3]; theta2 <- params[4]
  z1 <- data[1]; z2 <- data[2]; m01 <- data[3]; m02 <- data[4]
  y01 <- data[5]; y02 <- data[6]; N0 <- data[7]; N <- data[8]
  matrix(
    c(
      -(m01 / phi^2) -
        (z1 * (1 - theta1)^2) / (theta2 + (1 - theta1) * phi)^2 -
        (z2 * theta1^2) / (1 - theta2 + theta1 * phi)^2,
      -N - N0,
      (z1 * (1 - theta1) * phi) / (theta2 + (1 - theta1) * phi)^2 -
        z1 / (theta2 + (1 - theta1) * phi) -
        (z2 * theta1 * phi) / (1 - theta2 + theta1 * phi)^2 +
        z2 / (1 - theta2 + theta1 * phi),
      -((z1 * (1 - theta1)) / (theta2 + (1 - theta1) * phi)^2) +
        (z2 * theta1) / (1 - theta2 + theta1 * phi)^2,
      -N - N0,
      -((m01 + m02 + z1 + z2) / lam2^2),
      0,
      0,
      (z1 * (1 - theta1) * phi) / (theta2 + (1 - theta1) * phi)^2 -
        z1 / (theta2 + (1 - theta1) * phi) -
        (z2 * theta1 * phi) / (1 - theta2 + theta1 * phi)^2 +
        z2 / (1 - theta2 + theta1 * phi),
      0,
      -((m01 - y01) / (1 - theta1)^2) -
        y01 / theta1^2 -
        (z1 * phi^2) / (theta2 + (1 - theta1) * phi)^2 -
        (z2 * phi^2) / (1 - theta2 + theta1 * phi)^2,
      (z1 * phi) / (theta2 + (1 - theta1) * phi)^2 +
        (z2 * phi) / (1 - theta2 + theta1 * phi)^2,
      -((z1 * (1 - theta1)) / (theta2 + (1 - theta1) * phi)^2) +
        (z2 * theta1) / (1 - theta2 + theta1 * phi)^2,
      0,
      (z1 * phi) / (theta2 + (1 - theta1) * phi)^2 +
        (z2 * phi) / (1 - theta2 + theta1 * phi)^2,
      -((m02 - y02) / (1 - theta2)^2) -
        y02 / theta2^2 -
        z1 / (theta2 + (1 - theta1) * phi)^2 -
        z2 / (1 - theta2 + theta1 * phi)^2
    ), nrow = 4, ncol = 4
  )
}