Algorithm for uniform sampling with bounded replacement

Question

Is there a simple algorithm to sample from the uniform distribution on sequences of $n$ numbers, each taking one of $m$ integer values from $0$ to $m-1$, where each value can be repeated at most $r$ times? (For this to make sense, we need $rm\geq n$.)

Some examples:

• If $r \geq n$ then the task is sampling with replacement. The elements of the sequence are independent, and each one is uniformly distributed, so we can just call something like rand() % m, $n$ times.
  Note: This test case rules out any algorithms that unconditionally decrease the probability of drawing a value twice in a row. When $r\geq n$, the elements must be drawn independently.
• If $r=1$ then the task is sampling without replacement. We can select the $i$th element of the sequence in order, choosing uniformly over the $n-i$ values that have not been chosen already.
• If $rm=n$, then it suffices to choose a random permutation on $n$ symbols
• If $n=3, m=2, r=2$, then we should generate each of the following six sequences with equal probability: 001, 010, 011, 100, 101, 110. In particular, an initial 0 should be followed by a 1 with probability 2/3.
• If $n=3, m=3, r=2$, then we should generate each of the following 24 sequences with equal probability: 001, 002, 010, 011, 012, 020, 021, 022, 100, 101, 102, 110, 112, 120, 121, 122, 201, 202, 210, 211, 212, 220, 221. The sequence should include all three values with probability $6/24=1/4$.
  Note: This test case is pretty strong. It can rule out algorithms that generate sequences in a symmetric fashion, roughly speaking, and pass the preceding tests.

It would be nice to be able to choose each value in order from the conditional distribution given the previous values (and marginalized over future values), but I don't know how to express that.

We can perform rejection sampling: sample with full replacement and reject sequences that don't fit the $m$ criterion. But for $rm\approx n$, this method will be very slow.

My question is motivated by this stackoverflow question: How to generate evenly distributed n-permutation from range 0 to m

Answer 1

The answer is yes.

Tabulating such a sample produces a sequence of $m$ integers with values between $0$ and $r$ inclusive. "Uniform" suggests making the counts exchangeable in the sense that any sample in which the integer $i$ appears $k_i$ times, which can be summarized as the vector $\mathbf{k} = (k_0,k_1,k_2, \ldots, k_{m-1}),$ is just as likely as a sample summarized by any permutation of the components of $\mathbf{k}.$

We may picture such a vector as a tableau in which $n$ cells of an $m\times k$ checkerboard are filled, with $k_i$ filled cells in row $i.$ Given such a tableau, create a list consisting of $k_0$ $0$'s, $k_1$ $1$'s, and so on, then randomly permute that list: there is a random sample.

It remains to create a random tableau. Simply sample $n$ of its cells without replacement.

In R this can be accomplished with this function:

rtableau <- function(n, m, r)  (sample.int(m * r, n) - 1) %% m

It samples without replacement from the integers $0,1,\ldots, mr-1$ and reduces them modulo $m.$ The properties of this procedure are so obvious they scarcely need stating or proving.

Here are some tableaux for random samples with $n=11, m=5, r=3.$ Red cells are the selected ones:

Answer 2

Would this work (Algorithm 1):

1. Take the set $\{ 0, 1, \ldots, m-1\}$ and replicate it $r$ times, so that you have $\{ 0, 0, \ldots \mbox{$r$ times} \ldots 0, 1, 1, \ldots \mbox{$r$ times} \ldots, 1, \ldots, m-1\}$.
2. Permute this set of size $r \times m$ randomly.
3. Take the first $n$ items.

This is a functional procedure that produces the samples, but

1. It is probably woefully inefficient, and
2. It does not provide the full listing of all possible samples.

Another version which I also think is going to work is (Algorithm 2):

1. Initialize the step counter $k=0$. Initialize the probabilities of selection as $(p^{(0)}_0, p^{(0)}_1, \ldots, p^{(0)}_{m-1}) = ( r/(rm), r/(rm), \ldots, r/(rm)$ so that all entries have the same probability.
2. Draw the $k$-th element of your sequence with probabilities $(p^{(k)}_0, p^{(k)}_1, \ldots, p^{(k)}_{m-1})$. Let's say the $k$-th element thus drawn is $v_k$.
3. Increment $k \leftarrow k+1$
4. Update the selection probabilities for the next step: $p^{(k)}_{v_k}=p^{(k-1)}_{v_k}-1/(rm)$, and rescale all probabilities so that they sum up to 1, ${\bf p} \times={\bf p} \, (rm-k+1)/(rm-k) $. That is, if you drew $v_k=0$, go to the first entry (well... zeroth entry in C indexing notation), and reduce the probability that it is taken for the next entry. After it is taken $r$ times, this category is exhausted, and you cannot take from that category anymore.
5. If $k=n$, stop, otherwise go back to step 2.

I think (although I can't prove it) that this algorithm will produce the right probabilities. It is possible that with a smart choice of the ways to draw the next element, you can enumerate all possible sequences: instead of random numbers, draw entries in a deterministic cycle.

Answer 3

This is an unnatural selection method in that the probabilities of the early draws depends on the total to be drawn. So your $n=3, m=3, r=2$ example gives a probability that the first two to be drawn are $00$ of $\frac1{12}$ while reducing $n$ to $2$ would give a probability that the first two to be drawn being $00$ of $\frac19$; your future intentions are affecting the past probabilities.

It is as if you want to generate all the possible permitted sequences ($24$ in your final example) and then just select one of them. That will actually work as an algorithm (let's call it A1), though it will be inefficient with large numbers, probably inevitable given your condition that the future affects the past.

Two alternative approaches which do not have future intentions affecting the past, so giving different probabilities, would involve

• A2: starting with $m$ bins with $r$ balls each, and successively choosing a non-empty bin and taking a ball from that: the probability of choosing a ball from bin $i$ is $0$ if that bin is empty and $1$ divided by the number of non-empty bins if not.

• A3: starting with $r$ balls of each of the $m$ types, and successively choosing a ball without replacement: the probability of choosing a ball of type $i$ is the the number of remaining balls of type $i$ divided by the total number of remaining balls.

Both of these allow you to go on and increase $n$ later without having to start again; it means that the distributions for the early selections do not depend on $n$.

A3 has the feature that if $n=mr$ (so a full selection without replacement) then is gives the same distribution as A1 without having to generate all $\frac{(mr)!}{(r!)^m}$ possibilities in advance.