Generate random numbers summing to a predefined value

Generate 4 random numbers, compute their sum, divide each one by the sum and multiply by 40.

If you want Integers, then this will require a little non-randomness.

Use multinomial distribution

from numpy.random import multinomial
multinomial(40, [1/4.] * 4)

Each variable will be distributed as a binomial distribution with mean n * p equal to 40 * 1/4 = 10 in this example.

b = random.randint(2, 38)
a = random.randint(1, b - 1)
c = random.randint(b + 1, 39)
return [a, b - a, c - b, 40 - c]

(I assume you wanted integers since you said "1-40", but this could be easily generalized for floats.)

Here's how it works:

• cut the total range in two randomly, that's b. The odd range is because there are going to be at least 2 below the midpoint and at least 2 above. (This comes from your 1 minimum on each value).
• cut each of those ranges in two randomly. Again, the bounds are to account for the 1 minimum.
• return the size of each slice. They'll add up to 40.

Here's the standard solution. It's similar to Laurence Gonsalves' answer, but has two advantages over that answer.

1. It's uniform: each combination of 4 positive integers adding up to 40 is equally likely to come up with this scheme.

and

2. it's easy to adapt to other totals (7 numbers adding up to 100, etc.)

import random

def constrained_sum_sample_pos(n, total):
    """Return a randomly chosen list of n positive integers summing to total.
    Each such list is equally likely to occur."""

    dividers = sorted(random.sample(range(1, total), n - 1))
    return [a - b for a, b in zip(dividers + [total], [0] + dividers)]

Sample outputs:

>>> constrained_sum_sample_pos(4, 40)
[4, 4, 25, 7]
>>> constrained_sum_sample_pos(4, 40)
[9, 6, 5, 20]
>>> constrained_sum_sample_pos(4, 40)
[11, 2, 15, 12]
>>> constrained_sum_sample_pos(4, 40)
[24, 8, 3, 5]

Explanation: there's a one-to-one correspondence between (1) 4-tuples (a, b, c, d) of positive integers such that a + b + c + d == 40, and (2) triples of integers (e, f, g) with 0 < e < f < g < 40, and it's easy to produce the latter using random.sample. The correspondence is given by (e, f, g) = (a, a + b, a + b + c) in one direction, and (a, b, c, d) = (e, f - e, g - f, 40 - g) in the reverse direction.

If you want nonnegative integers (i.e., allowing 0) instead of positive ones, then there's an easy transformation: if (a, b, c, d) are nonnegative integers summing to 40 then (a+1, b+1, c+1, d+1) are positive integers summing to 44, and vice versa. Using this idea, we have:

def constrained_sum_sample_nonneg(n, total):
    """Return a randomly chosen list of n nonnegative integers summing to total.
    Each such list is equally likely to occur."""

    return [x - 1 for x in constrained_sum_sample_pos(n, total + n)]

Graphical illustration of constrained_sum_sample_pos(4, 10), thanks to @FM. (Edited slightly.)

0 1 2 3 4 5 6 7 8 9 10  # The universe.
|                    |  # Place fixed dividers at 0, 10.
|   |     |       |  |  # Add 4 - 1 randomly chosen dividers in [1, 9]
  a    b      c    d    # Compute the 4 differences: 2 3 4 1