# Random numbers that sum up to specific value

I would like to generate a list with random numbers, which add up to a specific value.

While[Total[x] == 28, x = RandomInteger[{0, 28}, 5];Print[x]]


The random number list should be returned if their sume is 28. Unforunately, this loop does not work. :( Can anybody help me?

-
How they could be random if they are related by Total[x] == 28? I mean, they are not independent. p.s. Not efficient -RandomChoice @ IntegerPartitions[28, {5}] –  Kuba Jul 10 '14 at 12:24
addressing only what is wrong with your loop, x is initially undefind, thus not 28, so you never enter the loop. You want While[Total@x=!=28,x = RandomInteger[{0, 28}, 5]];Print[x]; (note the print now outside the loop) –  george2079 Jul 10 '14 at 12:33
note there is a pretty challenging question of what exactly do you mean by random. –  george2079 Jul 10 '14 at 12:47
–  David Carraher Jul 10 '14 at 13:12
@Kuba See my last comment on the answer. –  Szabolcs Jul 10 '14 at 13:32

## 5 Answers

There are much better programming methods in Mathematica than loops.
Here is an approach based on IntegerPartitions, it chooses 5 numbers that sum up to 35:

RandomChoice[ IntegerPartitions[35, {5}]]

{12, 10, 7, 5, 1}


If we don't use RandomChoice it will write all 5-tuples, there are

IntegerPartitions[35, {5}] // Length

674


of them

One could also do it with FrobeniusSolve[{1, 1, 1, 1, 1}, 35] however the latter yields all permutations (including zeros):

FrobeniusSolve[{1, 1, 1, 1, 1}, 35] // Length

82251

-
i'm pretty sure an integer partition with repeated values should occur less frequently than one with all unique values (obviously depending on the definition of random). –  george2079 Jul 10 '14 at 12:54
+1 Nice and clean. You should perhaps wrap the function in RandomSample to avoid necessarily sorting from largest to smallest. –  David Carraher Jul 10 '14 at 13:05
@DavidCarraher Whether ordering is important should be mentioned in the question. If it is important, we can't quite just randomly permute any of the lists returned by IntegerPartitions. Instead we need to assign them weights based on how many permutations they have. For example, {2,2,2} has one permutation while {3,2,1} has 6 permutations (both sum to 6). Then choose a list based on this weight. –  Szabolcs Jul 10 '14 at 13:29
Generally, to solve these types of problems (in theory), we need to generate all possible outcomes (sets of numbers) that satisfy the OP's constraint (they sum to a certain value). Then select one randomly (i.e. assign the same probability to each). One important question will be to decide which element are considered distinguishable, e.g. does ordering matter? –  Szabolcs Jul 10 '14 at 13:31
I guess the max entropy distribution makes sense for the OP since it's the "most random" and he wanted random :P –  Rojo Jul 13 '14 at 2:27

How about...

...making 4 random 'integer cuts' or slicings of the line segment from 0 to n? I used this approach, though more crudely, in an earlier code-golf challenge: http://codegolf.stackexchange.com/questions/8574/generating-n-unique-random-numbers-with-a-specific-sum. (Chenminqi proposed a similar analysis for the present challenge but for some reason deleted it.)

riggedRandom[sum_,nPartitions_]:=
Differences@{0,Sequence@@Sort@RandomChoice[Range[0,sum],nPartitions-1],sum}


Example:

sum = 35; n = 5;
Print[n, " random numbers that sum to ", sum, ":\t" ,
numbers = riggedRandom[sum, n]]
sum == Total@numbers


5 random numbers that sum to 35: {9,1,1,14,10}
True

Analysis:

bounds = Prepend[Accumulate@numbers, 0];
sliceAt = Most@Accumulate@numbers;
intervals =
List[{{#[[1]], 1}, {#[[2]], 1}}] & /@ Partition[bounds, 2, 1];
labelPositions =
intervals /. {{a_, b_}, {c_, d_}} :> Sequence[(a + c)/2, 1.5];
labels = Text @@@ Thread[List[numbers, labelPositions]];


Below, each subsegment is indicated by a bi-directional arrow. The random numbers are the lengths of each subsegment.

Slices or cuts were made at the red points on the number line (from 0 to 35) underneath.

Graphics[{AbsolutePointSize[5], Red, Point[{#, 0}] & /@ sliceAt,
Arrowheads[{-.02, .02}], Arrow /@ intervals,
Black, labels}, BaseStyle -> 12, ImageSize -> 500,
Axes -> {True, False}, Ticks -> {bounds, {}},
PlotRange -> {{0, 35}, {0, 3}}]


-
your code does't work with me (sum_ and nPartitions_ don't appear in the function body) –  eldo Jul 10 '14 at 15:24
Oops. Parameter name change at last minute.... Fixed. –  David Carraher Jul 10 '14 at 16:34
+1: This is the approach endorsed by whuber and the statistics StackExchange. It also runs in only $O(n\log n)$ time. –  Rahul Jul 10 '14 at 16:53
This dis-favors results with zeros for some reason compared to the brute force result (It is dramatically faster than my solution though ) –  george2079 Jul 10 '14 at 19:17
Strange about disfavoring zeros. Zeros occur due to the repetition of "cuts". –  David Carraher Jul 10 '14 at 19:36

Let us arbitrarily set some details (because they were not specified): Allow all positive integers $\alpha$ and find the probability $p_\alpha$ with which each must be selected so that $\sum_{\alpha=0}^\infty \alpha\, p_\alpha = N$ for given $N$.

Applying Jaynes' maximum entropy principle (that Szabolcs mentioned), we find that the "least biased" estimate is

$$p_\alpha = \frac{1}{N+1}\left(1+\frac{1}{N}\right)^{-\alpha}$$

for which indeed $\sum_\alpha \alpha\,p_\alpha=N$ and $\sum_\alpha \,p_\alpha=1$. It looks like this:

Manipulate[
DiscretePlot[(1 + 1/n)^-α/(1 + n), {α, 0, 100},
PlotRange -> {0, .1}],
{{n, 10}, 1, 500}
]


EDIT: This solves the more general problem "here are $N$ numbers and $M$ operators, make it so we get $P$ in the end". A special case is probably this question (with only additions and integers allowed).

-
I am having trouble understanding the relationship to the question. You found a distribution whose expectation is N. I assume N is the specific value the OP wants to get as total. If one can draw each (nonnegative) number with probability palpha (so no number is allowed twice), then N would be the expectation of the total. But even a single number with this distribution could be higher than N. –  Rojo Jul 13 '14 at 5:20
It's a joke. The question shows that so little effort (undefined variables, not specifying what numbers are allowed etc) that I thought it would be funny (to me) to answer with a method that gives the least biased guess for a distribution consistent with the available (incomplete information). So, here it is. Anyway the idea with this distribution is that it's the least biased way to choose the distribution if all you know is the mean. Basically, if an ensemble of OPs were to pick lots of numbers at random for a long time, they'd get $N$ on average with a decreasing variance. –  acl Jul 13 '14 at 14:13

Brute force approach, using 9 for a somewhat smaller example:

 s93 = Select[ RandomInteger[{0, 9}, {10^6, 3}] , Total@# == 9 &];


a look at the statistics of the results:

 SortBy[Tally[Sort /@ s93], #[[1]] &]

 {{{0, 0, 9}, 2920}, {{0, 1, 8}, 6048}, {{0, 2, 7}, 5940},
{{0, 3, 6}, 6174}, {{0, 4, 5}, 6067}, {{1, 1, 7}, 2995},
{{1, 2, 6}, 6025}, {{1, 3, 5}, 5947}, {{1, 4, 4}, 3083},
{{2, 2, 5}, 2884}, {{2, 3, 4}, 5955}, {{3, 3, 3}, 1019}}


What you see is the {3,3,3} occurs relatively rarely, because obviously there is only one permutation.

Now this is a mod to @Artes IntegerPartition approach (allowing zeros) , that I think gets the statistics of the distribution correct (ie same as the brute force approach )

 n=Length@s93 (* generate same number of sets as first example yielded *)
s93prime =
RandomSample /@
RandomChoice[((Length@Permutations@#) & /@ #) -> #, n] &@
(#~Join~ConstantArray[0, {3 - Length[#]}] & /@
IntegerPartitions[9, 3]);
SortBy[Tally[Sort /@ s93prime], #[[1]] &]

 {{{0, 0, 9}, 3013}, {{0, 1, 8}, 5855}, {{0, 2, 7}, 6121},
{{0, 3, 6}, 6076}, {{0, 4, 5}, 5947}, {{1, 1, 7}, 3008},
{{1, 2, 6}, 5974}, {{1, 3, 5}, 6145}, {{1, 4, 4}, 2953},
{{2, 2, 5}, 2999}, {{2, 3, 4}, 5943}, {{3, 3, 3}, 1023}}


looks the same statistically..

( I'm sure there is a simple formula for Length@Permutations@list which would speed it up even more )

Edit: this works (Length@#)!/Product[ i! , { i , Last /@ Tally@#}] &@list , but i suppose there is a built in combinatorial function to do even better

For completeness this is the While approach (since i was a little off in my comment )

 s93w = Table[
x = {};
While[Total[x] =!= 9, x = RandomInteger[{0, 9}, 3]] ;
x, {n}];

-

The first part of this answer is probably silly, we can probably do much better using the function RandomComposition from Combinatorica instead

The package Combinatorica has the functionRandomPartition. I have checked out the code and made a compiled version out of it. It's not very nice asPartitionsP cannot be compiled and that is used all the time, but it is still faster than the function from Combinatorica.

randomPartition =
Compile[
{{n, _Integer}}
,
Module[
{mult = Table[0, {n}], j, d, r = n, z, res, counter, counter2},
res = ConstantArray[0, n];
counter = 1;
counter2 = 1;
While[
r > 0,
d = 1; j = 0; z = RandomReal[] r PartitionsP[r];
While[
z >= 0,
j++;
If[r - j d < 0, j = 1; d++;]; z -= j PartitionsP[r - j d];
];
r -= j d; mult[[j]] += d;
];

Do[
counter2 = counter + mult[[j]] - 1;
Do[
res[[i]] = j;
,
{i, counter, counter2}
];
counter = counter2 + 1
,
{j, Length[mult]}
];
res
]
,
CompilationTarget -> "C"
]


Then we can generate a random list of fixed length with the following function

randomWeakComposition =
Function[{nn},
SystemPermute[randomPartition[nn], SystemRandomPermutation[nn]]
]


Example

randomWeakComposition[10]

{0,1,1,1,0,3,1,1,1,1}


I'm not really happy about the speed, but what can you do

nn = 100;
mm = 10000;
Do[randomWeakComposition[nn], {mm}] // Timing // First

1.912838


Alternative

You can also use the multinomial distribution, which leads to a different distribution of the outcomes.

multDistr = MultinomialDistribution[nn, ConstantArray[1/nn, nn]];


Example

nn = 20;
RandomVariate[MultinomialDistribution[nn, ConstantArray[1/nn, nn]]]

{2,1,1,1,0,1,1,1,1,1,0,2,0,0,3,0,2,1,0,2}


Timing

mm = 100000;
(all = RandomVariate[multDistr, {mm}]) // Timing // First

2.286845
`

Which is a lot faster.

-