R/likelihood_p_closed.R

In BriceonWiley/IntegratedLikelihood.R:

likelihood_p_closed <- function(p, data) {
  n00 <- data[1]; n01 <- data[2]; n10 <- data[3]
  n11 <- data[4]; X <- data[5]; Y <- data[6]

  indices <- expand.grid(j = 0:Y, k = 0:X)
  j <- indices$j; k <- indices$k

  probs <-
    (n00 + n01 + Y + X - j - k) * log(1 - p) +
    (n10 + n11 + j + k) * log(p)
  combinations <-
    lchoose(Y, j) +
    lchoose(X, k)
  betas <-
    lbeta(n00 + Y - j + 1, n01 + X - k + 1) +
    lbeta(n10 + j + 1, n11 + k + 1)
  sum(exp(probs + combinations + betas))
}