Multivariate Wallenius’ Non-Central Hypergeometric distribution in R

biased_urn.R

# https://en.wikipedia.org/wiki/Wallenius%27_noncentral_hypergeometric_distribution#Multivariate_distribution
# https://markheckmann.github.io/notes/wallenius-distribution.html
library(dplyr)
library(tibble)
library(BiasedUrn)
set.seed(1234); R = rexp(10)

ressim = t(replicate(1e5L,sort(sample.int(10, size = 3, prob = R))))
pmf = data.frame(ressim) %>% 
  group_by_all() %>% 
  tally() %>% 
  ungroup() %>% 
  mutate(p = n/sum(n))

Xs = integer(length(R))
Xs[unlist(pmf[2,1:3])] = 1L
p = dMWNCHypergeo(x = Xs, m = rep(1, length(R)), n=3, odds=R)
print(p)
#[1] 0.03290516
print(pmf[2,]$p)
#[1] 0.03297

Xs = integer(length(R))
Xs[unlist(pmf[5,1:3])] = 1L
p = dMWNCHypergeo(x = Xs, m = rep(1, length(R)), n=3, odds=R)
print(p)
#[1] 0.0128722
pmf[5,]$p
#[1] 0.01293