Calculate density of the multivariate Fisher's noncentral hypergeometric distribution

Density_FMVNCHypergeo.R

library(BiasedUrn)

softmax <- function(x){
  mx = max(x)
  ex = exp(x-mx)
  ex/sum(ex)
}


MVNCH_Fisher <-  function (x, m, weight) {
  x <- seq_len(x)
  n <- length(x)
  x0 <- x
  x <- as.integer(x)
  m <- as.integer(m)
  e <- 0
  h <- m
  a <- seq_len(m)
  count <- round(choose(n, m))
  out <- matrix(x[a], nrow = m, ncol = count)
  prob <- numeric(count)
  r <- x[a]
  prob[1] <- sum(weight[r])
  i <- 2L
  nmmp1 <- n - m + 1L
  while (a[1L] != nmmp1) {
  if (e < n - h) {
      h <- 1L
      e <- a[m]
      j <- 1L
    }else {
      e <- a[m - h]
      h <- h + 1L
      j <- 1L:h
    }
    a[m - h + j] <- e + j
    r <- x[a]
    out[, i] <- r
    prob[i] <- sum(weight[r])
    i <- i + 1L
  }
  list(X=out, prob=softmax(prob))
}


K=10
M=3
set.seed(123); OR = rexp(K)
dist_F = MVNCH_Fisher(K,M,log(OR))
theo_F = numeric(nrow(dist_F$X))
for(i in 1:ncol(dist_F$X)){
  x = numeric(K)
  x[dist_F$X[,i]] = 1
  theo_F[i] = dMFNCHypergeo(x = x, m = rep(1, K), n = M, odds = OR)
}

plot(dist_F$prob, theo_F, col=hcl.colors(length(dist_F$prob)))
abline(0, 1, lty=2)