Random generation from multivariate hypergeometric distribution

Question

James E. Gentle (2003) in Random Number Generation and Monte Carlo Methods describes the following algorithm for random generation from multivariate hypergeometric distribution:

To generate a multivariate hypergeometric random deviate, a simple way is to work with the marginals. The generation is done sequentially. Each succeeding conditional marginal is hypergeometric. To generate the deviate, combine all classes except the first in order to form just two classes. Next, generate a univariate hypergeometric deviate $x_1$. Then remove the first class and form two classes consisting of the second one and the third through the last combined, and generate a univariate hypergeometric deviate based on $N − x_1$ draws. This gives $x_2$, the number of the second class. Continue in this manner until the number remaining to be drawn is $0$ or until the classes are exhausted. For efficiency, the first marginal used would be the one with the largest probability.

This procedure however seems to give biased results for the class that is at the end of the queue, since it is possible to sample only from what is left from the previous draws, or am I mistaken? Is there any better way of sampling from this distribution?

Below I paste the naive implementation of the algorithm:

rng_mvhyper <- function(n, k) {
  N <- sum(n)
  m <- length(n)
  n_otr <- N - n[1]

  x <- rhyper(1, n[1], n_otr, k)
  for (i in 2:m) {
    n_otr <- n_otr - n[i]
    k <- k - x[i-1]
    x[i] <- rhyper(1, n[i], n_otr, k)
  }
  x
}

And the results of simulation:

> table(replicate(1e6, which.max(rng_mvhyper(c(2,2,2), 4))))

     1      2      3 
400246 333303 266451 

While classes have the same size in the population, being in the beginning of the queue seems to make values of such class more common in the sample drawn.

Answer 1

Answer #1 (March 04, 15:15GMT)

Here is a corrected R code

rng_mvhyper <- function(n, k) {
  N <- sum(n)
  m <- length(n)
  n_otr <- N - n[1]

  x <- rep(0,m)
  x[1] <- rhyper(1, n[1], n_otr, k)
  for (i in 2:(m-1)) {
    n_otr <- n_otr - n[i]
    k <- k - x[i-1]
    x[i] <- rhyper(1, n[i], n_otr, k)
  }
  x[m]=k-x[m-1]
  x
}

with a simulation experiment:

> apply(replicate(1e4, rng_mvhyper(c(27,27,27), 6)),1,mean)/6
[1] 0.3360667 0.3302500 0.3336833

Answer #2 (March 04, 21:14GMT)

The original code of the question has been modified into the current version and does not show any bias: when running the code rng_mvhyper given in the question above, there is no discrepancy with the expected frequencies. For instance

> apply(replicate(1e6, rng_mvhyper(c(27,27,27), 6)),1,mean)/6
[1] 0.3332098 0.3334538 0.3333363

The issue with the observed bias is the use of which.max in the verification, which induces the bias rather than the formal solution proposed by James Gentle, which has nothing to do with bias. It is a valid representation of the multivariate Hypergeometric distribution, which is only defined when $k\le\sum_i n_i$. In which case it is impossible that "the remainder of marbles to be drawn is greater than number of marbles in the last category".