R/likelihood_phi_closed.R

In BriceonWiley/IntegratedLikelihood.R:

likelihood_phi_closed <- function(phi, data) {
  # browser()
  z1 <- data[1]; z2 <- data[2]; m01 <- data[3]; m02 <- data[4]
  y01 <- data[5]; y02 <- data[6]; N0 <- data[7]; N <- data[8]

  first <- -(m01 + m02 + z1 + z2 + 0.001) * log((phi + 1) * (N0 + N) + 0.001)
  indices <- expand.grid(i = 0:z1, j = 0:z2)
  i <- indices$i; j <- indices$j
  combs <- choose(z1, i) * choose(z2, j)
  delt <- phi^(m01 + i + z2 - j)
  betas <- beta(m01 - y01 + i + 1, y01 + z2 - j + 1) *
    beta(m02 - y02 + j + 1,  y02 + z1 - i + 1)
  second <- combs * delt * betas
  exp(first) * sum(second)
}