R/likelihood_gam_closed.R

In BriceonWiley/IntegratedLikelihood.R:

likelihood_gam_closed <- function(gam, data1, data2) {
  # browser()
  n000 <- data1[1]; n001 <- data1[2]; n010 <- data1[3]; n011 <- data1[4]
  X0 <- data1[5]; Y0 <- data1[6]; n100 <- data2[1]; n101 <- data2[2]
  n110 <- data2[3]; n111 <- data2[4]; X1 <- data2[5]; Y1 <- data2[6]
  N0 <- sum(data1); N1 <- sum(data2)

  indices <- expand.grid(i = 0:X0, j = 0:Y0, k = 0:X1, l = 0:Y1)
  i <- indices$i; j <- indices$j; k <- indices$k; l <- indices$l
  ijkl <- i + j + k + l

  first <- exp(gam)^(-n110 - n111 - k - l)
  combs <-
    choose(X0, i) * choose(Y0, j) * choose(X1, k) * choose(Y1, l)
  betas <-
    beta(Y0 + n000 - j + 1, X0 + n001 - i + 1) *
    beta(Y1 + n100 - l + 1, X1 + n101 - k + 1) *
    beta(n010 + j + 1, n011 + i + 1) *
    beta(n110 + l + 1, n111 + k + 1)
  gammas <-
    gamma(n000 + n001 + n100 + n101 + X0 + Y0 + X1 + Y1 - ijkl + 1) *
    gamma(n010 + n011 + n110 + n111 + ijkl + 1)

  hyp <-
    Vectorize(hypergeometric2F1, "a")(
      n010 + n011 + n110 + n111 + ijkl + 1, N1, 2 + N0 + N1,
      1 - cosh(gam) + sinh(gam), method = "Laplace", log = FALSE
      )
  sum(first * combs * betas * gammas * hyp)
}