Generating N uniform random numbers that sum to M

Question

This question has been asked before, but I've never really seen a good answer.

1. I want to generate 8 random numbers that sum to 0.5.

2. I want each number to be randomly chosen from a uniform distribution (ie, the simple function below will not work because the numbers will not be uniformly distributed).

   def rand_constrained(n,tot):
       r = [random.random() for i in range(n)]  
       s = sum(r)
       r = [(i/s*tot) for i in r] 
       return r
   

The code should be generalizable, so that you can generate N uniform random numbers that sum to M (where M is a positive float). If possible, can you please also explain (or show with a plot) why your solution generates random numbers uniformly in the appropriate range?

Related questions that miss the mark:

Generate multiple random numbers to equal a value in python (current accepted answer isn't uniform--another answer which is uniform only works with integers)

Getting N random numbers that the sum is M (same question in Java, currently accepted answer is just plain wrong, also no answers with uniform distribution)

Generate N random integers that sum to M in R (same question, but in R with a normal--not uniform--distribution)

Any help is greatly appreciated.

Answer 1

What you are asking for seems to be impossible.

However, I will re-interpret your question so that it makes more sense and is possible to solve. What you need is a probability distribution on the seven-dimensional hyperplane x_1 + x_2 + ... + x_8 = 0.5. Since the hyperplane is infinite in extent, a uniform distribution on the whole hyperplane will not work. What you probably(?) want is the chunk of the hyperplane where all the x_i>0. That region is a simplex, a generalization of a triangle, and the uniform distribution on the simplex is a special case of the Dirichlet Distribution.

You may find this section of the Dirichlet Distribution Wikipedia article, string cutting, particularly illuminating.

Implementation

The Wikipedia article gives the following implementation in Python in the Random Number Generation section:

params = [a1, a2, ..., ak]
sample = [random.gammavariate(a,1) for a in params]
sample = [v/sum(sample) for v in sample]

What you probably(?) want is the case where all ai=1 which results in a uniform distribution on the simplex. Here k corresponds to the number N in your question. To get the samples to sum to M instead of 1, just multiply sample by M.

Update

Thanks to Severin Pappadeux for pointing out that gammavariate can return infinity under rare circumstances. That is mathematically "impossible" but can occur as an artifact of the implementation in terms of floating point numbers. My suggestion for handling that case is to check for it after sample is first calculated; if any of the components of sample are infinity, set all the non-infinity components to 0 and set all the infinity components to 1. Then when the xi are calculated, outcomes like xi=1, all other x's=0, or xi=1/2, xj=1/2, all other x's=0 will result, collectively "corner samples" and "edge samples".

Another very-low-probability possibility is for the sum of the gammavariates to overflow. I would guess that we could run through the entire underlying pseudo-random number sequence and not see that happen, but theoretically it is possible (depending on the underlying pseudorandom number generator). The situation could be handled by rescaling sample, e.g., dividing all elements of sample by N, after the gammavariates have been calculated but before the x's are calculated. Personally, I wouldn't bother because the odds are so low; program crashes due to other reasons would have higher probability.