Author: michaelnhowes

The hypergeometric distribution describes the number of white balls in a sample of \(n\) balls draw without replacement from an urn with \(m_W\) white balls and \(m_B\) black balls. The probability of having \(x\) white balls in the sample is:

\(\displaystyle{f(x; m_W, m_B, n) = \frac{\binom{m_W}{x} \binom{m_B}{n-x}}{\binom{m_W+m_B}{n}}}\)

This is because there are \(\binom{m_W+m_B}{n}\) possible ways to sample \(n\) balls from the urn which has \(m_W+m_B\) black balls in total. Of these possibilities, there are \(\binom{m_W}{x} \binom{m_B}{n-x}\) ways of selecting exactly \(x\) white balls.

The multivariate hypergeometric distribution is a generalization of the hypergeometric distribution. In the multivariate hypergeometric distribution, there are \(k\) colours instead of \(2\). Suppose that the urn contains \(m_1\) balls of the first colour, \(m_2\) balls of the second colour and so on. The total number of balls in the urn is \(M = m_1+m_2+\cdots + m_k\). Next we draw \(n\) balls from the urn. Let \(X_i\) be the number of balls of colour \(i\) in our sample. The vector \(X=(X_1,X_2,\ldots, X_k)\) follows the multivariate hypergeometric distribution. The vector \(X\) satisfies the constraints \(X_i \in \{0,1,\ldots,M_i\}\) and \(\sum_{i=1}^k X_i = n\). The probability of observing \(X=x\) is

\(\displaystyle{f(x; m, n) = \frac{\binom{m_1}{x_1} \binom{m_2}{x_2}\cdots \binom{m_k}{x_k}}{\binom{M}{n} }, \quad x_1+x_2+\cdots+x_k = n}\)

This is because there are \(\binom{M}{n}\) possible ways of sampling \(n\) balls from the urn that has \(M\) balls in total. Likewise, The number of ways of choosing \(x_i\) balls of colour \(m_i\) is \(\binom{m_i}{x_i}\).

Sampling

There are two common methods for sampling from the multivariate hypergeometric distribution. The first option is to sample \(n\) numbers from \(1,\ldots,M\) without replacement. The numbers \(1,\ldots,M\) represent the balls in the urn. We can let \(1,\ldots,m_1\) represent the numbers with colour 1, and let \(m_1+1,\ldots,m_1+m_2\) represent the numbers with colour 2 and so on. Suppose that \(Y_1,\ldots,Y_n\) is the sample of size \(n\) without replacement from \(1,\ldots,M\). Next define:

\(\displaystyle{X_i = \#\{j : X_1+X_2+\cdots+X_{i-1} < Y_j \le X_1+X_2+\cdots+X_{i-1}+X_i\}}\)

Then the vector \(X=(X_1,\ldots,X_k)\) will be distributed according to the multivariate hypergeometric distribution. This method of sampling has complexity on the order of \(n\log(M)\). This is because sampling a random number of size \(M\) has complexity \(\log(M)\) and we need to sample \(n\) of these numbers.

There is a second method of sampling that is more efficient when \(n\) is large and \(k\) is small. This is a sequential method that samples each \(X_i\) from a (univariate) hypergeometric distribution. The idea is that when sampling \(X_1\) we can treat the balls of colour \(1\) as “white” and all the other colours as “black”. This means that marginally, \(X_1\) is hypergeometric with parameters \(m_W = m_1\) and \(m_B = m_2+\cdots + m_k\) and \(n = n\). Furthermore, conditional on \(X_1\), the second entry, \(X_2\), is also hypergeometric. The parameters for \(X_2\) are \(m_W = m_2\), \(m_B = m_3+\cdots+m_k\) and \(n=n-X_1\). In fact, all of the conditionals of \(X\) are hypergeometric random variables and so can sampled sequentially. The python code at the end of this post implements the sequential sampler.

The complexity of the sequential method is bounded above by \(k\) times the complexity of sampling from the hypergeometric distribution. By using the ratio of uniform distribution this can be done in complexity on the order of \(\log(m_w+m_b)\). The overall complexity of the sequential method is therefore on the order of \(k\log(M)\) which can be much smaller than \(n\log(M)\).

import numpy as np

def multivariate_hypergeometric_sample(m, n):
    x = np.zeros_like(m)
    k = len(m)
    for i in range(k-1):
        x[i] = np.random.hypergeometric(m[i], np.sum(m[(i+1):]), n)
        n -= x[i]
    x[k-1] = n
    return x

I was working on a problem involving the maximum of a collection of geometric random variables. To state the problem, let \((X_i)_{i=1}^n\) be i.i.d. geometric random variables with success probability \(p=p(n)\). Next, define \(M_n = \max\{X_i : 1 \le i \le n\}\). I wanted to know the limiting distribution of \(M_n\) as \(p \to 0\) and \(n \to \infty\). I was particularly interested in the case when \(p\) was on the order of \(1/n\).

It turns out that this problem has a very nice solution: If \(p\to 0\) as \(n \to \infty\), then for all real numbers \(c\)

\(\displaystyle{\lim_{n \to \infty}\mathbb{P}\left(M_n \le \frac{\log n}{p}+\frac{c}{p}\right) \to \exp(-e^{-c})}\)

That is, \(M_n\) is on the order of \(\frac{\log n}{p}\) and has fluctuations of order \(1/p\). Furthermore, these fluctuations are distributed according to the Gumbel distribution. Importantly, this result places no restrictions on how quickly \(p\) goes to \(0\).

In the rest of this post, I’ll derive the above identity. The first step relates \(M_n\) to the maximum of a collection of exponential random variables. The second step uses a coupling to prove that the exponential and geometric random variables are close to each other.

An easier version

When \(p\) goes to \(0\), the scaled geometric random variables \(p X_i\) converge to exponential random variables. Exponential random variables are nicer to work because they have a continuous distribution.

This suggests that we could try replacing all the geometric random variables with exponential random variables. More precisely, we will replace \(X_i\) with \(Y_i/p\) where \(Y_i\) are exponentially distributed. Next, we will derive the main result for the exponential random variables. Then, in the next section, we will prove that \(X_i\) is sufficiently close to \(Y_i/p\).

Let \(\widetilde{M}_n = \max\{ Y_i/p : 1 \le i \le n\}\) where \(Y_i\) are i.i.d and exponentially distributed. This means that for \(y \ge 0\), \(\mathbb{P}(Y_i \le y)=1-e^{-y}\). The condition \(\widetilde{M}_n \le y\) is equivalent to \(Y_i \le py\) for all \(i =1,\ldots,n\). Since the \(Y_i\)’s are independent, this implies that for \(y \ge 0\)

\(\displaystyle{\mathbb{P}(\widetilde{M}_n \le y) = \prod_{i=1}^n \mathbb{P}(Y_i \le py) = \left(1-e^{-py}\right)^n}\)

Now fix a real number \(c\) and let \(y = \frac{\log n}{p}+\frac{c}{p}\). Since \(n \to \infty\), we know that \(y \ge 0\) for large enough \(n\). Thus, for sufficiently large \(n\)

\(\displaystyle{\mathbb{P}\left(\widetilde{M}_n \le \frac{\log n}{p}+\frac{c}{p}\right) = \left(1-e^{-\log n – c})\right)^n = \left(1-\frac{e^{-c}}{n}\right)^n\to \exp(-e^{-c})}\)

And so the main result holds for \(\widetilde{M}_n = \max\{Y_i/p : 1 \le i \le n\}\). We want to show that it holds for \(M_n = \max\{X_i : 1 \le i \le n\}\). The next section proves a close relationship between \(X_i\) and \(Y_i\).

Exponential and geometric coupling

It is true that as \(p \to 0\), \(pX_i\) converges in distribution to an exponential random variable. This is a result that holds for a single random variable. For an approximation like \(M_n \approx \widetilde{M}_n\), we need a result that works for every \(X_i\) at once. This seems tricky since it involves the convergence of an infinite collection of random variables. Fortunately, there is a simple approach based on coupling.

For each geometric random variable \(X_i\) it is possible to produce an exponential random variables \(Y_i\) with

\(\displaystyle{\frac{1}{\lambda} Y_i \le X_i \le \frac{1}{\lambda} Y_i + 1},\)

where \(\lambda = -\log(1-p)\). This is because \(\mathbb{P}(X_i \le x) = 1-(1-p)^{[x]}\) and \(\mathbb{P}(Y_i/\lambda \le x) = 1-e^{-x/\lambda} = 1-(1-p)^{x}\). This means that

\( \displaystyle{\mathbb{P}(Y_i/\lambda \le x-1) \le \mathbb{P}(X_i \le x)\le \mathbb{P}(Y_i/\lambda \le x-1) }.\)

If we sample both \(X_i\) and \(Y_i\) with inverse transform sampling, and use the same uniform random variables, then we will have

\(\displaystyle{ \frac{1}{\lambda} Y_i \le X_i \le \frac{1}{\lambda} Y_i + 1},\)

as claimed. If we use this inequality for every single pair \((X_i,Y_i)\) we will get

\(\displaystyle{\max\{Y_i/ \lambda\} \le M_n \le \max\{Y_i/\lambda\} +1}\)

As \(p \to 0\), we have that \(\lambda \to 0\) and \(\frac{p}{\lambda}=\frac{p}{-\log(1-p)} \to 1\). This means that the random variable \(\max\{Y_i/\lambda\} \) has the same limit as \(\widetilde{M}_n = \max\{ Y_i/p : 1 \le i \le n\}\). Furthermore, \(\max\{Y_i/\lambda\} +1\) also has the same limit as \(\widetilde{M}_n\). Since \( M_n\) is “sandwiched” between \(\max\{Y_i/\lambda\}\) and \(\max\{Y_i/\lambda\} +1\), the result must also hold for \(M_n\). And so

\(\displaystyle{\lim_{n \to \infty}\mathbb{P}\left(M_n \le \frac{\log n}{p}+\frac{c}{p}\right) \to \exp(-e^{-c})}\)