Let $n,p \in \mathbb{N}_+$ with $p \leq n.$ Let $\mathcal{P}$ denote the set of partitions of $\{1, \ldots, n\}$ into $p$ nonempty sets. How can I efficiently sample uniformly from $\mathcal{P}$?

$\begingroup$ I would imagine that a Markov chain with rapid mixing might be useful here, though it would give only approximately uniform distribution. $\endgroup$– DmitryZSep 12, 2013 at 21:24

$\begingroup$ How big n are you interested in? Do you ask because you want one method, or because you actually want to do computations on random ones? It might be worth to change computer language, for speed, even. $\endgroup$– Per AlexanderssonMar 27 at 17:32

1$\begingroup$ There is a substantial literature on this problem. See, for example, djalil.chafai.net/blog/2012/05/03/…. $\endgroup$– Ira GesselMar 27 at 19:16
3 Answers
The sets you're interested in are counted by Stirling Numbers of the Second Kind, which satisfy the recursion $$\left\{{n \atop k}\right \}=\left\{{n1 \atop k1}\right \}+k \left\{{n1 \atop k}\right \}$$ Here the first term represents those partitions where $n$ is its own set, and the remaining term represents inserting $n$ into a partition of $\{1, \dots, n1\}$. This recursion can also be used to generate a set partition recursively:
With probability $\left\{{n1 \atop k1}\right \}/\left\{{n \atop k}\right \}$ put $n$ in its own set, and make the rest a partition of $n1$ elements into $k1$ sets chosen uniformly at random. Otherwise generate a uniform random partition of $n1$ elements into $k$ sets, and insert $n$ into a set uniformly chosen from those sets.
The algorithm would run in time on the order of $nk$, with the main overhead being computing (or looking up) all of the Stirling Numbers up to $\left\{{n \atop k}\right \}$ at the start of the algorithm.

$\begingroup$ Hmmm. Efficiency issues aside, this algorithm isn't practical for large n and k unless you're using a CAS that can compute the requisite Stirling numbers and then compute the quotient to reasonable accuracy. $\endgroup$– AatGSep 25, 2013 at 20:32

1$\begingroup$ @AatG The desired output is a partition of $\{1,\ldots,n\}$, which is going to take on the order of $n$ units of memory even to write down. That puts an upper limit on how big $n$ can be. If you have a computer that has $n$ units of memory then it will likely be powerful enough to carry out this computation. A CAS is not really needed; you just need support for highprecision arithmetic. $\endgroup$ Mar 27 at 15:21

1$\begingroup$ It might be enough to only compute this approximately $\endgroup$ Mar 27 at 16:23
I have provided some Python code for a quick and dirty implementation of the algorithm described by Kevin in the hope that it will save someone some time. The randomPartitionsFixedP() function is the top level function.
Code:
import numpy as np
from scipy.special import stirling2 as s2
# 
def randomPartitionsFixedP(n,p, num_partitions):
# This function returns num_partitions random partitions of n elements into p nonempty subsets.
# These partitions are sampled uniformly randomly from the set of all such partitions
partitions = []
for i in range(num_partitions):
partitions.append(randomPartition(n,p))
return partitions
# 
def randomPartition(n,p):
# Returns a single random partition according to the uniformly random algorithm
# Procedure 0:
#  Place the first n1 objects into p1 nonempty sets, and the remaining singleton in its own set.
#  This procedure terminates the descent into the tree when p=2
#  When working up from the bottom, add next object to its own set
# Procedure 1:
#  Place the first n1 objects into p nonempty sets, and the remaining singleton into one of these at random.
#  This procedure terminates the descent into the tree when n=p+1
#  When working up from the bottom add this object to an existing set uniformly randomly
# Recursively select procedures according to the stirlingbased probabilities
n_tracker = n
p_tracker = p
procedure_list = []
while (p_tracker>=2 and n_tracker>=p_tracker+1):
procedure = gumbelSampling(n_tracker, p_tracker)
procedure_list.append(procedure)
if procedure == 0:
p_tracker = 1
n_tracker = 1
else:
n_tracker = 1
# Reverse this list as we work from the bottom of the decision tree up.
procedure_list.reverse()
# Initialise the partition assignment list
object_partition_assignments = []
if procedure_list[0] == 0:
object_partition_assignments.extend([0]*n_tracker+[1])
max_partition = 1
else:
object_partition_assignments.extend(list(range(n_tracker)))
object_partition_assignments.extend(np.random.choice(object_partition_assignments, size=1))
max_partition = np.max(object_partition_assignments)
# Now continue back up the decision tree
for procedure in procedure_list[1:]:
if procedure==0:
object_partition_assignments.append(max_partition+1)
max_partition += 1
else:
object_partition_assignments.append(np.random.choice(np.unique(object_partition_assignments), size=1)[0])
return object_partition_assignments
# 
def gumbelSampling(n,p):
# returns 0 or 1 (selecting procedure 0 or 1 in the recursive algorithm according to stirling number based probabilities)
procedures = [0,1]
unnormalised_log_probs = computeLogProbabilities(n,p)
new_probs = np.random.gumbel(loc=0, scale=1, size=len(procedures)) + unnormalised_log_probs
return procedures[np.argmax(new_probs)]
# 
def computeLogProbabilities(n,p):
# computes the log probabilities for p and 1p, where p=S(n1,p1)/S(n,p)
# returns a list [p, 1p] which will be provided to the Gumbel max sampling function
S_a = s2(n1,p1)
S_b = s2(n,p)
return [np.log(S_a)np.log(S_b), np.log(S_bS_a)np.log(S_b)]
Example (Code):
n_parts = 100000
n = 4
p = 3
parts = randomPartitionsFixedP(n = n,p = p, num_partitions=n_parts)
youneek = np.unique(parts, axis=0, return_counts=True)
print("Number of Unique Partitions: ", s2(4,3))
print("The Unique Partitions: \n", youneek[0])
print("Frequency of Each Partition: ", youneek[1]/n_parts)
Example (Output):
Number of Unique Partitions: 6.0
The Unique Partitions:
[[0 0 1 2]
[0 1 0 2]
[0 1 1 2]
[0 1 2 0]
[0 1 2 1]
[0 1 2 2]]
Frequency of Each Partition: [0.16595 0.16816 0.16616 0.16611 0.16696 0.16666]
K.C. Locey, Random integer partitions with restricted numbers of parts
An algorithm is presented to generate uniform random samples of integer partitions for a total $Q$ with $N$ parts from the set of all $P(Q, N)$ partitions.
This algorithm was developed by Ken Locey in response to a 2010 question on StackOverflow.

$\begingroup$ Thanks, but I want to sample set partitions, not integer partitions. $\endgroup$– AatGSep 13, 2013 at 17:57

$\begingroup$ If $N=a_1+\dots + a_k$ is a random integer partition, and $x_1,x_2,\dots,x_N$ is a random permutation, then $\{ \{ x_1,\dots,x_{a_1}\}, \{x_{a_1+1},\dots,x_{a_1+a_2} \}, \dots \}$ is a random set partition. It isn't clear to me that it won't be uniformly random. $\endgroup$ Sep 13, 2013 at 21:41

3$\begingroup$ @KevinO'Bryant it can't be uniformly random, because otherwise we would have $P(n)B_n$ ($P(n)$ being the number of partitions and $B_n$ the Bell numbers, i.e. the number of set partitions), and it's easy to see this doesn't hold  in fact, it already breaks down for $n=3$ when there are $3$ partitions ($3$, $2+1$, $1+1+1$) and $5$ set partitions ($123$, $123$, $213$, $312$, $123$). $\endgroup$ Sep 14, 2013 at 0:31