Algorithm to generate N random numbers between A and B which sum up to X

Question

This problem seemed like something which should be solvable with but a few lines of code.

Unfortunately, once I actually started to write the thing, I've realized it's not as simple as it sounds.

What I need is a set of X random numbers, each of which is between A and B and they all add up to X. The exact variables for the problem I'm facing seem to be even simpler: I need 5 numbers, between -1 and 1 (note: these are rational (floating point) numbers), which add up to 1.

My initial "few lines of code, should be easy" approach was to randomize 4 numbers between -1 and 1 (which is simple enough), and then make the last one 1-(sum of previous numbers). This quickly proved wrong, as the last number could just as well be larger than 1 or smaller than -1.

What would be the best way to approach this problem?

PS. Just for reference: I'm using C#, but I don't think it matters. I'm actually having trouble creating a good enough solution for the problem in my head.

---

I also wanted to provide my current solution to the problem, but do remember that it's quite imperfect and it was created as a quick fix to the initial problem!

• Generate 4 random numbers between -1 and 1
• Create a "final" number X=SUM(previousNumbers)
• If the final number is >1 or <-1, then:
  • Get the "amount" over 1 / under -1 and make the last number a 1 / -1
  • Find another number which can accept this amount and still inside the brackets
  • If no number can take the amount (it's too large / too small) then divide the amount in half and try again for each half
• Randomize the order of the generated numbers and return them

This works in the sense that this algorithm generates 5 numbers which are between -1 and 1 and their sum is 1. However the downside is that quite often one of the generated numbers is a 1 (or -1), which doesn't feel very random.

Answer 1

As said before, this question actually doesn't have an answer: The restrictions imposed on the numbers make the randomness questionable at best.

However, you could come up with a procedure that returns a list of numbers like that:

Let's say we have picked the first two numbers randomly as -0.8 and -0.7. Now the requirement is to come up with 3 'random' numbers that sum up to 2.5 and are all in the range [-1,1]. This problem is very similar to the starting problem, only the dimensions have changed. Now, however, if we take a random number in the range [-1,1] we might end up with no solution. We can restrict our range to make sure that solutions still exist: The sum of the last 2 numbers will be within the range [-2,2]. This means we need to pick a number in the range [0.5,1] to make sure we can reach a total of 2.5.

The section above describes one step in the process.

In general: Determine the range for the next number by applying the range of the rest of the numbers to the required sum. Pseudo-code:

function randomNumbers (number, range, sum) {
  restRange = range * (number - 1)
  myRange = intersection ([sum - restRange.upper, sum - restRange.lower], range)

  myNumber = random(myRange)

  rest = randomNumbers (number-1, range, sum - myNumber)

  return [myNumber, rest]
}

So for the case described above

randomNumbers (3, [-1,1], 2.5)
  restRange = [-1,1] * (3-1) = [-2,2]
  myRange = intersection ([2.5-2,2.5-(-2)], [-1,1]) = intersection ([0.5,4.5],[-1,1]) = [0.5,1]

A quick-and-dirty implementation in Java:

public class TestRandomNumberList
{

    @Test
    public void test()
    {
        double[] numbers = new double[5];
        randomNumbers(numbers, 0, -1, 1, 1);
        assertEquals(sum(numbers), 1.0, 0.00001);
        for (double d : numbers)
        {
            assertTrue(d >= -1 );
            assertTrue(d <= 1);
        }
    }

    private void randomNumbers(double[] numbers, int index, double lowerBound, double upperBound, double sum)
    {
        int next = index + 1;
        if (next == numbers.length)
        {
            numbers[index] = sum;
        }
        else
        {
            int rest = numbers.length - next;  

            double restLowerBound = lowerBound * rest;
            double restUpperBound = upperBound * rest;

            double myLowerBound = Math.max(lowerBound, sum - restUpperBound);
            double myUpperBound = Math.min(upperBound, sum - restLowerBound);

            numbers[index] = random(myLowerBound, myUpperBound);
            randomNumbers(numbers, next, myLowerBound, myUpperBound, sum - numbers[index]);
        }
    }

    private double random(double myLowerBound, double myUpperBound)
    {
        double random = Math.random();
        return myLowerBound + random * (myUpperBound - myLowerBound);
    }

    private double sum(double[] numbers)
    {
        double result = 0;
        for (double num : numbers)
        {
            result += num;
        }
        return result;
    }

}

Answer 2

Simple, as long as you know how many.

1. You need N numbers called V1 to Vn. The required sum is S.
2. Generate N random numbers (in any convenient range). They are R1 to Rn.
3. Calculate their sum as SR.
4. Scale each number so Vn = Rn * S / SR.

You may produce a tiny rounding error, but I doubt this will be a problem.

If the number N is supposed to be random, then choose that first.

---

My apologies, I missed the requirement for the numbers to be between A and B. The algorithm is exactly the same. Pick N random numbers, then scale them based on A, B, actual sum and required sum. I leave the rest as an implementation detail.

Answer 3

It is possible for each variable to be uniformly distributed on an interval while the joint distribution is satisfies a codimension 1 constraint. For example, if you pick (x,y,z) so that it is uniformly distributed on a unit sphere, then each coordinate is uniformly distributed on the interval [-1,1]. The coordinates are just not independent.

Even if you give up independence, it is not possible for 5 numbers uniformly distributed on [-1,1] to have a constant sum of 1. This is because expectation is linear for any random variables, not just for independent ones. If you have 5 random variables which are each uniformly distributed on [-1,1], their sum has average value 0, so it can't be the constant 1.

Some other answers suggest picking the first numbers to be uniform on [-1,1], and then fix the last few numbers. This usually gives up the symmetry between the numbers. You might be able to tell which numbers were generated without constraints, and which were used to make the sum fit, since the 5th number might be in [-1,1] but it might not have a uniform distribution.

Instead, you might take a conditional distribution. Because the problem is symmetric, the conditional distribution is symmetric. Imagine that you choose some epsilon>0, sample from [-1,1] 5 times independently, and then reject the 5-tuple if the sum is more than epsilon away from 1. This gives you some joint distribution within the unit 5-cube that is concentrated near the constraint hyperplane x0+x1+x2+x3+x4=1. Then let epsilon go to 0, and you get a limiting distribution. Since the probability that you get a point on the plane is too low (mathematically 0, but positive for, say, random 32 bit floats) to use rejection sampling directly, you need another method for producing this distribution.

It is equivalent to produce {xi/2 + 1/2} uniform on [0,1] with sum 3, or {xi/6+1/6} uniform on [0,1/3] with sum 1. So, produce 5 positive numbers summing to 1, and then reject the sample (repeat if any of these are greater than 1/3. (Afterwards, multiply these by 6 and subtract 1 to get numbers on [-1,1] summing to 1.) To produce 5 positive numbers summing to 1, one way is to generate 4 random numbers on [0,1], sort them (y0,y1,y2,y3), add 0 and 1 to the front and back, (0,y0,y1,y2,y3,1), and then take the differences (y0,y1-y0,y2-y1,y3-y2,1-y3). (Another method is to create 5 exponentially distributed random variables by taking -log(uniform), and renormalize by their sum.)

You have to be wary of rejection sampling in high dimensions since the probability of acceptance might be low, and then you have to iterate too many times. The probability that among the 5 numbers summing to 1, none of them is at least 1/3 can be found with inclusion-exclusion: 1-5(2/3)^4+10(1/3)^4 = 11/81 > 1/8. So, the average number of 5-tuples you need to generate by this method before you find one with no number greater than 1/3 is under 8, so on average, it takes fewer than 40 uniformly random numbers to generate a 5-tuple with sum 1 drawn from the conditional distribution. If the parameters of the problem change, you might prefer a more complicated method that avoids rejection sampling.

Answer 4

This question is maybe old, but here is my solution (written in C#), which is based on Maarten Winkels Java code:

public static double[] GenerateRandomNumbers(uint values, double minimum, double maximum, double sum, Random generator = null)
    {
        if (values == 0)
            throw new InvalidOperationException($"Cannot create list of zero numbers.");
        if (minimum * values > sum)
            throw new InvalidOperationException($"The minimum value ({minimum}) is too high.");
        if (maximum * values < sum)
            throw new InvalidOperationException($"The maximum value ({maximum}) is too low.");
        if (minimum > maximum)
            throw new InvalidOperationException($"The maximum value ({maximum}) is lower than the minimum value ({minimum}).");
        if (generator == null)
            generator = new Random();

        var numberList = new double[values];

        for (var index = 0; index < values - 1; index++)
        {
            var rest = numberList.Length - (index + 1);

            var restMinimum = minimum * rest;
            var restMaximum = maximum * rest;

            minimum = Math.Max(minimum, sum - restMaximum);
            maximum = Math.Min(maximum, sum - restMinimum);

            var newRandomValue = generator.NextDouble(minimum, maximum);
            numberList[index] = newRandomValue;
            sum -= newRandomValue;
        }

        numberList[values - 1] = sum;

        return numberList;
    }

Code for generating random double values between a minimum and maximum value:

public static double NextDouble(this Random generator, double minimum, double maximum)
    {
        if (minimum > maximum)
            throw new InvalidOperationException($"The maximum value ({maximum}) is lower than the minimum value ({minimum}).");

        return generator.NextDouble() * (maximum - minimum) + minimum;
    }

And here is how it is used:

var min = -1.0;
var max = 1.0;
var sum = 0.0;
uint values = 100000;
var seed = 123;
var generator = new Random(seed);

var randomNumbers = Extensions.GenerateRandomNumbers(values, min, max, sum, generator);

Debug.WriteLine($"Distinct Values: {randomNumbers.Distinct().Count()}");
Debug.WriteLine($"Min: {randomNumbers.Min()}");
Debug.WriteLine($"Max: {randomNumbers.Max()}");
Debug.WriteLine($"Average: {randomNumbers.Average()}");
Debug.WriteLine($"Median: {randomNumbers.Median()}");
Debug.WriteLine($"Sum: {randomNumbers.Sum()}");

Debug.WriteLine("\nFirst 10 values:");
randomNumbers.Take(10).ToList().ForEach(v => Debug.WriteLine(v));

The output:

Distinct Values: 99800
Min: -0,999962684698385
Max: 1
Average: 0
Median: 0,00128587102577371
Sum: 0

First 10 values:
0,969113830462617
0,815630646336652
0,487091036274606
0,623283306892628
0,477558290342595
-0,903369966849391
-0,965998261219821
-0,701281160908416
-0,610592191857562
0,26017893536956

One problem I faced with Winkels' code was that it was a recursive method and for large lists (> 30 000 numbers) the program would throw StackOverflow exceptions. This is why I wrote it as a iterative function. I also tried to clean up and removed unnecessary operations. I also added support for seeding the random number generator, and some error handling.

There is still room for improvement. I haven't looked too much into how truly random this is, so investigate further if you need to do anything scientific.

Answer 5

You can start with any 5 random numbers in the range, but then you just need to pay attention to which are negative and which are positive.

Case (A) : All 5 numbers are positive.

Here you can just divide by their sum. The sum is guaranteed to be greater than each individual number, since each number N is positive. Then you have

0 < N < 1 

for each individual number N, and you have the normalized sum equal to 1.

Case (B) : Positives and negatives.

First normalize the negatives so that their sum is

-1 < (Sum of normalized negatives) < 0

At this point we are sure that each negative number is still within the range

-1 < N_negative < 0.

for each individual negative number N_negative.

Next normalize the positives so that

1 - (Sum of normalized positives) = (Sum of normalized negatives)

Now we might have got into trouble taking one or more of the positive numbers outside of the range. Check if that happens. If not, then we are done. If so, then renormalize the positive numbers so that the maximum individual value is 1. We know that the sum of renormalized positives will be less than the sum of the first positive normalization, but also the renormalized sum will be greater than 1. Therefore we can renormalize the negatives so that

1 - (Sum of *re*normalized positives) = (Sum of *re*normalized negatives)

And now we are definitely done with case B.

Case (C) : All 5 numbers are negative.

This has no solution.

Answer 6

One option to solve this problem is to build a search tree. The goal state will be the sum of all nodes along your current path equal X. Each node is a different value between the range, so the number of nodes at each level will be abs(A) + abs (B). Explore all the paths and get all the possible solutions. Then choose one at random from the list of solutions.

The search space will become too large to calculate when N and abs(A) + abs (B) are large. You can limit the search space by creating rules. A simple rule here would be that if the path contains the same values, but in a different order, then they are considered the same. That would eliminate a lot of paths on the search tree. You can think up a lot more to limit the search space, but I will leave that exercise to you.

Note: This solution is overkill, but it's a fun exercise.

Answer 7

I'd start out by simplifying the problem. If you have N numbers from [A,B], they will add up to a value between N*A and N*B. Just subtract A from each number and N*A from the target X.

So, without loss of generality we can look for N numbers between [0, B'] which add up to X'. Again, this can be simplified by dividing X' by B'.

The resulting problem is looking for N numbers in [0,1] which add up to X" = (X-A)/(B-A).

What's the distribution of these random numbers? The requirement that they add up to X" means that they can't be independent distributions. I'll assume that we want the same distribution for all N numbers. Otherwise, it's trivial. If we know all distributions are the same, we do know that their mean has to be X"/N. This is not necessarily 1/2, so we have to come up with a possibly asymmetric distribution over [0,1] - a symmetric distribution does have mean 0,5 since p(x)==p(1-x). Yet all numbers should be possible if X"/N is not 0 or 1.

At this point, you can freely choose any parametrized distribution which allows for means in (0,1). An exponential distribution should work fine.

Finally, once you've drawn one number A(0) from this distribution, you should recurse and calculate a new distribution for the next number as you now need to find N-1 numbers which add up to X"-A(0). And of course for the last number you have no choice left.