It’s easy to sample uniformly from the set of partitions of a set: you pick a number of bins using an appropriate exponential distribution, then randomly i.i.d. toss each element of the set into one of those bins. At the end of this procedure, the nonempty bins will constitute your uniformly sampled random partition. [literature ref: “Generation of a random partition of a finite set by an urn model” by Stam, 1983; pretty pictures ref: see this blog post].
Is there a similarly efficient and simple algorithm for sampling uniformly from the set of partitions into a fixed number of nonempty sets? The only algorithm I’m aware of takes advantage of the fact that the number of such partitions is given by \(\left\{n \atop p \right\}\), a Stirling number of the second kind, if the set has \(n\) elements and we desire \(p\) nonempty subsets in the partition. In particular, we have the identity
\[
\left\{ {n \atop p} \right\} = \left\{ {n-1 \atop p-1} \right\} + p \left\{ {n-1 \atop p} \right\}
\]
that comes from observing that the only two ways to construct a partition of \(n\) elements into \(p\) nonempty sets are to: either partition the first \(n-1\) elements into \(p-1\) nonempty sets and take the remaining singleton as our final set in the partition, or we partition the first \(n-1\) elements into \(p\) nonempty sets and place the remaining element into any of these \(p\) sets.
This observation leads to a straightforward recursive sampling procedure: with probability \(\left\{ {n-1 \atop p-1} \right\}/\left\{ {n \atop p} \right\}\), use the first procedure with a randomly sampled partition of the first \(n-1\) elements into \(p-1\) nonempty sets, otherwise use the second procedure with a randomly sampled partition of the first \(n-1\) elements into \(p\) nonempty sets.
Unfortunately, this is not an efficient procedure for several reasons. In particular, it requires computing Stirling numbers, and taking their ratio. When \(n,k\) are large, this will require both computational time and, more of a practical impediment, arbitrary-precision arithmetic. A straightforward implementation also relies on recursion, which is infeasible when \(n\) is large. Clearly one can implement this algorithm without using recursion; one can also use an asymptotic expansion of the Stirling numbers to approximate the ratio \(\left\{ {n-1 \atop p-1} \right\}/\left\{ {n \atop p} \right\}\) when \(n,k\) are large … but this comes at the cost of some unspecified inaccuracy and just doesn’t feel right.
So the question remains: is there an efficient and simple way to sample uniformly from the set of partitions into a fixed number of nonempty sets?