# Algorithm for uniform sampling with bounded replacement

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

• This sounds like a special case of the recently-asked question at stats.stackexchange.com/questions/129644. – whuber Dec 19 '14 at 19:35
• @whuber I don't think it is. In that question, the OP seems to be asking for something different: "the number of balls in the bucket has no influence on the outcome as long as there is a ball in the bucket". So for example, given $n=3,m=2,r=2$, it sounds like they want an initial 0 to be followed by a 1 with probability $1/2$. I'm asking for all legal sequences to have the same probability, which implies that an initial 0 is followed by a 1 with probability $2/3$. – Chris Culter Dec 19 '14 at 19:41
• Thanks. I'm still trying to get my head around what you want. It would help to fix the inconsistency in the last two bullets, which claim to apply to the same triple $(n,m,r)=(3,2,2)$. – whuber Dec 19 '14 at 19:44
• @whuber Yep, done. – Chris Culter Dec 19 '14 at 19:50

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.

• Hmm, Algorithm 1 doesn't work for my test case of $n=3,m=3,r=2$. It generates a permutation of 012 with probability $8\times3!\times3!/6!=2/5$, not $1/4$. – Chris Culter Dec 20 '14 at 0:37
• Algorithm 2 doesn't reduce to sampling with replacement for the simpler test case $n=2,m=2,k=2$. After an initial 0, it draws a 1 with probability $2/3$, not $1/2$. – Chris Culter Dec 20 '14 at 0:39