# 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?

-
It is unclear what the sum formula has to do with the clustering or even what "clustering" is. Here's one interpretation. Suppose the cluster sizes are $n_1, \ldots, n_k$. An obvious way to generate random partitions with these sizes is to permute all $N$ items and assign the first $n_1$ to the first cluster, the next $n_2$ to the next cluster, and so on. Given its obviousness, I presume you have rejected this as a solution, but why is it not a solution? –  whuber Mar 1 '12 at 18:13
Clustering in this case means a partition of the set of items. Of course, I should have emphasised that it was the cluster sizes I am really interested in, not a particular clustering. –  Rónán Daly Mar 1 '12 at 19:39
Thanks. In that case there are many possible solutions, depending on how the clusters are generated. If you don't know that, you need to begin by characterizing the distribution of the observed cluster sizes. If you could disclose either of these--data-generation process or distribution of cluster sizes--you will likely get answers that are more appropriate for your data. –  whuber Mar 1 '12 at 19:45

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.)

1. Draw independent samples from $y_j = Gamma(\alpha_j,1)$ for each $j = 1..k$.
2. Compute $x_j = y_j/\sum^k_{j=1}y_j$ for each $j$.
3. Then you have the sample $Dirichlet(\alpha_1,...,\alpha_k) = (x_1,...,x_k)$.
4. Multiply by the sample by $N$ to get the approximate new cluster sizes. $N(x_1,...,x_k)$.
5. Round the result to the nearest natural number. (make sure that the new cluster sizes add to $N$ due the rounding.)
6. 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.

-
This seems to be just what I was looking for. I had a look at the Dirichlet Process material as well, which leads me to believe I could generate an algorithm for $k$ not fixed. Thanks a lot for your detailed answer. –  Rónán Daly Mar 1 '12 at 22:53
Yep, a great thing about the Dirichlet Process is that it is a nonparametric model. –  Daniel Johnson Mar 2 '12 at 3:01

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.

-
You need to be careful with your Poisson rejection scheme. You want $\sum^k_{j=1}Pois(n_j) = N$ or the scheme will reject the proposal. Remember that the sum of Poissons is Poisson with the parameter being the sum of the individual problems. In this case, $P(\sum^k_{j=1}Pois(n_j) = N) = N^Ne^-N/N! \sim \frac{1}{\sqrt{2\pi N}}$ by Sterling. For N = 10,000 which is on the order of what @Ronan had in mind this probability is about .004 meaning that an average of about 10000/.004 = 250,000 Poisson random variables need to be computed for each trial. –  Daniel Johnson Mar 1 '12 at 23:07
@DanielJohnson, right. Let me think about it. –  StasK Mar 2 '12 at 1:43

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

-
I like this idea, but as it is now it doesn't fit the requirement that the new cluster sizes be close to the original. One way to achieve this would be to, instead of uniformly choosing an urn for a ball, assign each urn a probability proportional to the size of the original cluster and then use that categorical distribution when choosing which urn to put a ball into. This is intimately connected to the Dirichlet type process i described below, with the only disadvantage here is the lack of a parameter describing how 'close' the sample is to the original distribution on average. –  Daniel Johnson Mar 1 '12 at 19:41