Biased urn experiment and Fisher's noncentral hypergeometric distribution

Question

Let's have an urn with $m_i= 2$ balls for $c=3$ different colours with different weights $w_i$.

$n=2$ balls are taken randomly but the probability of sampling a particular coloured ball is proportional to its weight. This is a biased urn problem and the probability of taking a given sample follows the Multivariate Fisher's Noncentral Hypergeometric Distribution with the following probability mass function ($S$ is the support of the PMF):

$$ \frac{1}{P_0}\prod_{i=1}^{c} \binom{m_i}{x_i}\omega_i^{x_i}\\ P_0 = \sum_{(y_0,\ldots,y_c)\in \mathrm{S}}\prod_{i=1}^{c} \binom{m_i}{y_i}\omega_i^{y_i}\\ \mathrm{S} = \left\{ \mathbf{x} \in \mathbb{Z}_{0+}^c \, : \, \sum_{i=1}^{c} x_i = n \right\} $$

I would like to calculate the probability of the realization $x = \{1, 0, 1\}$.

The R's package BiasedUrn exposes the Multivariate Fisher's Noncentral Hypergeometric Distribution. Here a reproducible example:

library(BiasedUrn)

set.seed(123)
n = 2 # balls taken 
c = 3 # different colors
m = rep(2, c) # nuber of balls in the urn for each colour

w = runif(n = c)
w = w/sum(w) #weights of the balls (biasd urn)

# a realization
x = c(1, 0, 1)

prob = dMWNCHypergeo(x, m, n, odds = w)

The calculated probability is: $prob = 0.1209264$.

when I try to calculate it with the formula above I get different results:

numerator = 1.0
for (i in 1:c) {
    numerator = numerator * choose(m[i], x[i]) * (w[i] ^ x[i])
}

S = combn(x=rep(c(1:c),n), m = n, tabulate, nbins = c)
S = S[, !duplicated(t(S))]

denominator = 0.
for (y in S) {
    prod = 1
    for (i in 1:c) {
        prod = prod * choose(m[i], y) * (w[i]^y)
    }
    denominator = denominator + prod
}

prob_2 = numerator / denominator

and $prob\_2 = 0.1329621$. So, it's wrong!

While I am rather confident that the numerator is calculated in the right way, my instinct tells me that's something wrong in the denominator and probably I cannot identify in the right way the support of the PMF. I defined $S$ as the 6 possible (distinct) outcomes: $\{1,1,0\}$, $\{1,0,1\}$, $\{2,0,0\}$, $\{0,1,1\}$, $\{0,2,0\}$, $\{0,0,2\}$.

So, please, could someone help me find what I am doing wrong here?

I would like to generalize the problem and implement the Multivariate Fisher's Noncentral Hypergeometric Distribution in stan to make inference on $\omega_i$ (an old question of mine) and so I need help to understand what's I am doing wrong here.

Answer 1

The issue is in the numerical accuracy you are using to compute the values.

The numerator should be computed using logarithms and then exponentiate at the end.

So the numerator becomes (in python)

import numpy as np
from scipy.special import gammaln
def calc_num(mi,xi,wi):
    return gammaln(mi+1)-gammaln(xi+1)-gammaln(mi-xi+1) + xi*np.log(wi)
log_numerator = 0.0
for i in range(c):
    log_numerator += calc_num(m[i], x[i], w[i])
numerator = np.exp(log_numerator)

and the denominator

y = np.asarray([
    [2,0,0],
    [1,1,0],
    [1,0,1],
    [0,2,0],
    [0,1,1],
    [0,0,2]
])
denom=0.
for z in y:
    td = 0.
    for i in range(c):
        td += calc_num(m[i],z[i],w[i])
    denom += np.exp(td)

then

numerator / denom

yields 0.1329621176812538