How to sample natural numbers, such that the sum is equal to a constant? Say I have $N$ items that are partitioned / clustered and I want to randomly repartition these items, such that the distribution of sizes of the clusters is 'similar' to those that I already have. I'm viewing this (perhaps unhelpfully), as trying to sample $k$ natural numbers, such that $x_i \geq 1$ and $\sum_{i=1}^{k}x_i=N$. $N$ in this context would probably be in the thousands and the distribution of the current cluster sizes would be quite skewed, with a small number of large clusters and a large number of small clusters.
I've looked at hypergeometric, binomial, geometric distributions, but none of these seem to quite fit what I'm looking for - I'm guessing that it's more complicated than a simple distribution and some type of Markov-like process would be needed. Anybody have any ideas?
 A: One easy way to achieve your goal is to permute the labels. Say you had 10 objects, with memberships defined as $\{1, 2, 3, 4 ,5, 6\}$, $\{7, 8\}$, $\{9\}$ and $\{10\}$. You take a random permutation $\sigma=(3, 7, 2, 5, 1, 8, 10, 9, 6, 4)$, and then your new clusters are $\{\sigma(1), \sigma(2), \sigma(3), \sigma(4), \sigma(5), \sigma(6)\} = \{1, 2, 3, 5, 7, 8\}$, $\{\sigma(7), \sigma(8)\}=\{9, 10\}$, $\{\sigma(9)\}=\{6\}$ and $\{\sigma(10)\} = \{4\}$.
If you want to introduce some variability in the cluster sizes, that can probably be done, too -- e.g. by drawing clusters with sizes Poi(6), Poi(2), Poi(1) and Poi(1), rejecting the samples that do not add up to 10. (Poi($\lambda$) is a Poisson random variable with rate/expected value $\lambda$.)
(I wrote this before I read @whuber's comment. Oops)
Update, based on valuable @DanielJohnson's comment: For large values of the total, the procedure becomes impractical, as it will be rejecting most samples. What you would want to do, then, is to condition on the total number of objects, $N$, and then the Poisson distribution becomes a multinomial with probabilities $\lambda_1/L, \ldots, \lambda_k/L, \ldots$, where $L=\lambda_1 + \ldots + \lambda_k + \ldots$. So your samples would simply be multinomial ones. Some clusters of size 1 or 2 may get lost though, and if you don't like that, then again you could condition on all clusters being present. This would effectively introduce variability in the size of larger clusters, while the smaller ones will remain with their original sizes. Again, if you find yourself rejecting the samples too often, you can use a combination of: (1) maintaining clusters of size 1 or 2; (2) simulating the Poisson- or multinomial-distributed clusters of larger sizes.
A: I believe what you are looking for is something like the Dirichlet Process [DP] which is a distribution on distributions. It is not an easy concept to understand, but the base measure you will use is the discrete distribution of cluster sizes you started with. The parameter $\alpha$ controls how 'close' to the original distribution your new one is. Since a sample from a DP is a probability distribution (in your case, a discrete one), you can multiply it by $N$ to get cluster sizes. The result won't be integers, but just rounding the numbers should not affect what you're trying to do in a meaningful way. 
Edit: Somehow, I am more familiar with the Dirichlet Process than the Dirichlet Distribution, which is what you are actually looking for. The DP is an infinite dimensional generalization of the DD. 
To be more algorithmically precise, consult the subsection talking about random number generation. Your parameters are going to be based on the cluster sizes that you want the random clusterings to look like $\{n_j\}_{j=0}^k$. In other words: 
$$\alpha = (\alpha_1,...,\alpha_k) = (\frac{\beta n_1}{N},...,\frac{\beta n_k}{N})$$ where $\beta$ is called a concentration parameter and it controls how close to the original distribution of cluster sizes the Dirichlet distribution will be on average. Higher values of $beta$ will mean that the resulting DD will give closer and closer values to the EXACT distribution of sample sizes you started with. 
So, the algorithm would be: (if you have access to a function that can generate Dirichlet Distribution samples, ignore steps 1 and 2.)


*

*Draw independent samples from $y_j = Gamma(\alpha_j,1)$ for each $j = 1..k$. 

*Compute $x_j = y_j/\sum^k_{j=1}y_j$ for each $j$.

*Then you have the sample $Dirichlet(\alpha_1,...,\alpha_k) = (x_1,...,x_k)$. 

*Multiply by the sample by $N$ to get the approximate new cluster sizes. $N(x_1,...,x_k)$.

*Round the result to the nearest natural number. (make sure that the new cluster sizes add to $N$ due the rounding.)

*Fill the new clusters of a random permutation of elements.


The advantage that this has over the previously discussed option of just permuting the elements is that the cluster size distribution isn't always the same. You can allow for as much or a little variation from the original cluster sizes by controlling the concentration parameter $\beta$. In simulation this will likely result in a more robust calculation. 
A: Can you think of this as $n$ balls being distributed among $k$ urns?  That seems to fit your description of clusters (where you have $k$ clusters and $n$ numbers).  If you need at least one ball in each urn, then first put 1 ball in each urn, then randomly select the urn for each of the remaining $n-k$ balls.  Here is one possible implementation in R:
> nkballs <- function(n,k) {
+ tmp <- sample(k, n-k, replace=TRUE)
+ as.numeric( table(tmp) ) + 1
+ }
> 
> nkballs(1000, 25)
 [1] 36 47 40 52 40 28 34 43 37 35 38 33 45 37 45 45 37 34 38 46 42 34 56 46 32
> sum(.Last.value)
[1] 1000

