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:


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.


1 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([
for z in y:
    td = 0.
    for i in range(c):
        td += calc_num(m[i],z[i],w[i])
    denom += np.exp(td)


numerator / denom

yields 0.1329621176812538


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.