R/trinom.R

In Exact: Unconditional Exact Test

trinom <-
function(x12, x21, n, p12, p21, delta) {
  # See Hsueh et al. (2001) for trinomial with delta.  Look at H0u
  outLog <- (lgamma(n+1)-lgamma(x12+1)-lgamma(x21+1)-lgamma(n-x12-x21+1)) + 
    log((p12+delta)^x12) + log(p21^x21) + log((1-p12-p21-delta)^(n-x12-x21))
  return(exp(outLog))
  
  # One version that works, but has issues with large datasets:
  #exp(lgamma(n+1)-lgamma(x12+1)-lgamma(x21+1)-lgamma(n-x12-x21+1))*(p12+delta)^x12*p21^x21*(1-p12-p21-delta)^(n-x12-x21)
  
  # Another version that works, but is slower:
  #matTbls <- cbind(n-x12-x21, x12, x21, rep(0, length(x12)))
  #unlist(lapply(split(matTbls, seq(nrow(matTbls))), dmultinom, size=n, prob=c(1-p12-p21,p12,p21,0)), use.names = FALSE)
}