Try multivariate noncentral hypergeometric distribution in BiasedUrn package (without replacement)

MVNCNHypergeo.R

library(BiasedUrn)
library(dplyr)
K=10
set.seed(1234); OR = rexp(K)
M=3
ressim = replicate(100000, sort(sample.int(K, size = M, prob = OR, replace = FALSE))) %>% 
  t() %>% 
  as.data.frame() %>% 
  group_by_all() %>% 
  tally()

ressim = ungroup(ressim) %>% 
  mutate(prob = n/sum(n))

head(ressim)
theo_F = numeric(nrow(ressim))
theo_W = numeric(nrow(ressim))
i=1
for(i in 1:nrow(ressim)){
  x = numeric(K)
  x[unname(unlist(ressim[i,1:3]))] = 1
  theo_F[i] = dMFNCHypergeo(x = x, m = rep(1, K), n=M, odds = OR)
  theo_W[i] = dMWNCHypergeo(x = x, m = rep(1, K), n=M, odds = OR)  
}

plot(theo_F, theo_W)
abline(0,1,lty=2)

plot(ressim$prob, theo_F)
abline(0, 1, lty=2)

plot(ressim$prob, theo_W)
abline(0, 1, lty=2)