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 sum is 28. Unfortunately, this loop does not work.
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
RandomSample to avoid necessarily sorting from largest to smallest.
$\endgroup$
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.
$\endgroup$
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: https://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}}]
sum_ and nPartitions_ don't appear in the function body)
$\endgroup$
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).
As mentioned in comments and answer, you have not specified exactly how you want the random numbers to be distributed. One distribution that fits is the multinomial distribution
nn=28;
kk=5;
multDistr = MultinomialDistribution[nn, ConstantArray[1/kk, kk]];
Example
RandomVariate[multDistr]
{6, 7, 6, 4, 5}
Resurrecting an old question because I think we should start with:
Needs["Combinatorica`"]
RandomPartition[28]
RandomComposition[28, 5]
The RandomPartition command seems to do what Artes was doing (except not limited to 5 summands, more on that below). The RandomComposition seems to do what others were trying to do.
Although I'm not finding the right documentation to verify, I assume these choose a partition (composition) uniformly at random from all possible partitions (compositions) of 28 (into 5 parts, for the composition). I'm not sure if this was a feature that we didn't have at the time this question was asked (Artes didn't mention, so perhaps not, I'm not seeing dates in documentation). Maybe RandomPartition was just passed over because it doesn't allow us to specify the length of the partition?
So I will grant that the question is vague, but clearly, we should be able to improve Artes answer using a built-in command like RandomPartition if that command is written well at all. In principle, we should be able to do this faster than something that requires generating all possible partitions. And if the goal was a composition, RandomComposition seems to take care of that out-of-the-box.
Now, if you just need it to be exactly five, and you don't care about increasing the value 28 or 5, there is a crude way of doing it, similar to the original code:
x = {}; While[Length[x] != 5, x = RandomPartition[28]]; x
That will also be uniform (if you want a random integer 1 to 5, just roll a 6-sided die and throw away the 6s -- it's still uniform).
Before we go any further, I'll just define "partition" and "part." A partition of $n$ is a set of positive integers that add up to $n$ (taken "in no particular order", so $5=1+1+3$ is the same as $5=3+1+1$). Each of the summands is a part, so $5=3+1+1$ is a partition of 5 into 3 parts.
If you want, for example, a partition of 100 into 5 parts (or into 90 parts), that While loop takes much longer (on expectation) because there are relatively few partitions of 100 into exactly five (or 90) parts.
So let's talk a bit about efficiency. Say you want $n$ (e.g. $n=28$) to be partitioned randomly into $k$ (e.g. $k=5$) parts. It's been a while and I don't remember precisely, but for relatively feasible values of $n$, if $k$ is something like $\sqrt n$, you should be fine using the While loop.
Note: For small values of $n$, or relatively few computations, this doesn't take too long. So if you just want to do this for $n=28$ or something, you can stop here and you're fine. The While loop is fine -- use that. Or use Artes's solution. (Or if you're looking for slightly different type of output, use RandomComposition or perhaps Jacob Akkerboom's answer.) The speed is negligible at such low values, comparing thousandths of seconds to ten-thousandths.
For our own use, let's piece together a function to count the number of partitions of $n$ into $k$ parts:
myP[n_, k_] := myP[n, k] =
If[n <= 0 || k <= 0||k>n, 0,
If[k == 1 || k == n, 1,
myP[n - 1, k - 1] + myP[n - k, k]
]
];
We'll use that shortly to help us balance out our probabilities, assuming uniformity is still the goal. This is built up on a recurrence $p(n,k)=p(n-1,k-1)+p(n-k,k)$. Any partition of $n$ into $k$ parts ends in 1, in which case you can remove it and obtain a partition of $n-1$ into $k-1$ parts, or else you can remove 1 from each part, obtaining a partition of $p-k$ into $k$ parts. It's not hard to prove this is bijective, so the recurrence is exact.
One thing we're going to capitalize on is this idea of removing 1 from each part.
If $k$ is significantly smaller or larger than this nice middle ground (which is probably not exactly $\sqrt n$, but off the top of my head, that's the best estimate you'll get out of me), you'll want to be more careful.
First, say $k$ is very large. Many of the parts should be small. For example, if $n=101$ and $k=58$, many of the parts have to be 1s (in fact, 15 of them do, because $15\times 1 + 43\times 2 = 101$). So a random partition of 100 into 58 parts is likely to have 1s we want to ignore.
If we could ignore the 1s, we could subdivide the problem by randomly generating the (hopefully much smaller) remaining parts randomly. The issue is that there could be more than 15 of them that are 1. If it were always exactly 15, this would be easier, but in fact, k-1 of them could be 1.
But if we could remove all the 1s, we could take all those parts of 2 or more, remove 1 from them, and repeat this process. Remember, we want to do it randomly! And we want to be careful to weight things so that we don't have some issue with a non-uniform sampling.
Say we have $n=101$ and $k=58$. There are at least 15 parts that are 1, as we mentioned, but there could be more. If $k_1$ of them are 1s, the remaining $k-k_1$ make up a partition of $n-k_1$. We break that down into a new, smaller partition.
For example, say we have a smaller partition like: $$3+2+2+2+1+1+1+1+1+1+1+1+1+1+1+1+1+1.$$
We split off the parts that are 1: $$3+2+2+2((+1+1+1+1+1+1+1+1+1+1+1+1+1+1)).$$
We discard those to obtain: $$3+2+2+2,$$ a partition in which no part is 1.
We can make this correspond to the same partition with every part decreased by 1 (since we got rid of the 1s): $$2+1+1+1.$$
If we could randomly generate such a partition, we can use this to randomly generate our original partition, except that we need to bias the sample. We can't choose this parameter $k_1$ uniformly (indeed, $k_1$ can't be less than 15 for sure, but even so, it's not uniform among feasible choices).
We need to know the probability of each possible value of $k_1$ (which could be anything, in general, from none, all the way to $k_1$, although as we noticed, there can be some restrictions making $k_1$ large -- which is good for us, to have $k_1$ bigger).
But we know how to count this: There are $p(101-k,58-k_1)$ such partitions to build upon, and each of those will correspond with one of our possible random partitions. So we can tally up this list by making a table:
myCount=Table[myP[101-k,58-j],{j,0,57}]
We need to choose according to these, and we need them to be probabilities, so we normalize to 1:
myProb=N[Normalize[myCount, Total]]
We can use that to randomly sample values (using RandomChoice) of $k_1$. From there we can recursively build our partition that has high $k$ value.
To pack that up as a Module, we can do:
normalK[n_, k_] := Module[{part = {}},
While[Length[part] != k, part = RandomPartition[n]];
Return[part]
]
highK[n_, k_] := Module[{k1, prob, subpart = {}},
If[k < 6 Sqrt[n], Return[normalK[n, k]]];
prob = Normalize[Table[myP[n - k, k - j], {j, 0, k - 1}], Total];
k1 = RandomChoice[prob -> Range[0, k - 1]];
If[k - k1 < 6 Sqrt[n - k1] && n - k1 > 20,
subpart = 1 + normalK[n - k, k - k1],
subpart = 1 + highK[n - k, k - k1]
];
Return[Join[subpart, Table[1, {k1}]]];
];
That's the best I can come up with for $k$ relatively large. We can test it for speed, but I think we are optimistic already:
AbsoluteTiming[Do[normalK[100, 60];, {10}]]
(* Output: {44.2695, Null} *)
AbsoluteTiming[Do[highK[100, 60];, {10}]]
(* Output: {0.185296, Null} *)
It's way faster.
If $k$ is very small instead and the While loop runs too slowly, you'll have to sort of do the "opposite" -- write a similar module you could call lowK that looks at the largest part of the partition. If you know more about partitions, I'm suggesting essentially that lowK operates on the "conjugate" of the ideas in highK. Having one very large part is similar to having many 1s based on this concept, see: https://en.wikipedia.org/wiki/Partition_(number_theory)#Conjugate_and_self-conjugate_partitions
Ideally, you'd combine these two modules so that it's a single recursive module that can handle low, high, or normal ranges of k based on each of these ideas, but I'll leave those details to anyone who's come this far in my answer and wants to implement it.
So now does this stack up against other answers or other perspectives? Well, clearly, if order matters (or doesn't matter, depending on what "matters" means to you), this will not work -- and randomly permuting it will hurt our carefully balanced uniformity. It may be possible to sneak that into the RandomChoice probability distribution in highK, but at that point we'd also need to work up our own normalK too, since we wouldn't be able to rely on RandomPartition at all. But we've already seen RandomComposition, and Jacob Akkerboom gives one way to do it, coming from a different perspective.
If we are happy with the ordering but want 0s, it's quite easy to do one of:
normalK[n+k,k]-1
highK[n+k,k]-1
And just to prove that this wasn't a completely vacuous exercise, let's compare the timing:
AbsoluteTiming[Do[highK[100, 60];, {100}]]
(* Output: {0.74234, Null} *)
AbsoluteTiming[Do[RandomChoice[IntegerPartitions[100, {60}]];, {100}]]
(* Output: {3.28749, Null} *)
multDistr = MultinomialDistribution[100, ConstantArray[1/60, 60]];
AbsoluteTiming[Do[Reverse[Sort[RandomVariate[multDistr]]];, {100}]]
(* Output: {0.653411, Null} *)
AbsoluteTiming[Do[RandomComposition[100, 60], {100}];]
(* Output: {0.0142456, Null} *)
AbsoluteTiming[Do[RandomPartition[100], {100}];]
(* Output: {0.113123, Null} *)
So of course we're quite faster than normalK, as shown above, but also faster than Artes' method (saving some time not constructing all the possible partitions) and surprisingly competitive with Jacob Akkerboom's method (which surprised me). But RandomComposition trumps Jacob's algorithm, as we would expect, and it performs better than our highK (I certainly woudln't expect highK to be better, but it's at least comparable).
The Reverse and Sort are to make the output "the same" in some sense -- but they add negligible time anyway. And I'm letting Artes's code produce the table of values every single time because I want a good head-to-head comparison if, for example, you wanted to use it repeatedly with $n$ and $k$ of a consistent size. If you are literally using the exact same $n,k$ values over and over and Artes's method is feasible for those values, it may be cheaper to generate the table once and then make multiple random choices from it. But that's not my assumption.
If we increase $k$ a bit (e.g. 80), Jacob's method ends up being slower, and Artes's is faster. That makes sense because the list is very short. But if we increase $n$ a bit too (e.g. 400, with $k=320$), highK is still about as fast as Jacob's method, and Artes's seems infeasible.
If we increase $n$ substantially, we continue to perform well:
AbsoluteTiming[Do[highK[1000, 900];, {10}]]
(* Out= {0.348483, Null} *)
multDistr2 = MultinomialDistribution[1000, ConstantArray[1/900, 900]];
AbsoluteTiming[Do[Reverse[Sort[RandomVariate[multDistr2]]];, {10}]]
(* Out= {1.01725, Null} *)
AbsoluteTiming[Do[RandomComposition[1000, 900], {100}];]
(* Output: {0.0224496, Null} *)
AbsoluteTiming[Do[RandomPartition[1000], {100}];]
(* Output: {3.29139, Null} *)
I recommend no one attempt Artes's solution for $n=1000$. It does not scale well, but of course, I think everyone understood that was not the intention. Artes's solution would have provided a way to do highK (and lowK) but in a way that won't necessarily scale very well. But RandomPartition does, and so does our highK (and our theoretical lowK and whatever we build to tie all three together).
And there's a big surprise -- we're now beating RandomPartition even though we're doing more work (abstractly) than it.
If we fix $n-k=100$ and increase $n$ it seems like Jacob's timing increases (linearly?), as does RandomComposition at its much slower values, but ours doesn't seem to change at all. Of course, those two do something very different than our highK.
Anyway, the TLDR is that I've suggested using RandomComposition (if you want a composition, which seems to be the general consensus or at least plurality of interpretation), and that RandomPartition might be useful if you want a partition like in Artes's popular, but non-scalable, answer. These two Combinatorica functions are nice aesthetically, for low $n,k$, because they use built-in functionality, and they are very practical for higher $n,k$. And because we still don't get everything out of the box, I've whipped up a way of making RandomPartition work when what we really want is a nonexistent version that accepts a second argument, like RandomPartition[28,{5}].
"Combinatorica`" is a nice package, but to refer to it as built-in is perhaps a bit unsual and to expect the same performance from it as functions built into the kernel may lead to disappointment. (...)
$\endgroup$
Commented
Jun 29, 2017 at 8:34
"Combinatorica`" functions have been implemented in an older version of Mathematica, when there were fewer options to increase performance. Nowadays loading the "Combinatorica`" package gives error messages, so that is an additional argument against using "Combinatorica`" functions.
$\endgroup$
Commented
Jun 29, 2017 at 8:35
Combinatorica (besides the one about some of the functions being deprecated)? I don't see any, but perhaps it is different for you. As for my use of the term "built-in," I mostly mean "no need to write any code of our own." I know it's developed by a third-party, but it's certainly fast, and as the most popular Mathematica package ever (I certainly would estimate, at least), it is well-traveled and its use seems quite routine. You are right, of course, and I do understand, but it is no cause for alarm.
$\endgroup$
Commented
Jun 29, 2017 at 15:19
I didn't see anything in the original question that restricted the random numbers to integers (other than the original poster using RandomInteger). If you allow Reals, then a simple answer is
sum = 28;
num = 5;
res = {0, RandomReal[sum, num - 1], sum} // Flatten // Sort // Differences
(* {9.59902, 0.623308, 3.77447, 10.0473, 3.9559} *)
Total[res]
(* 28. *)
If the number of random numbers is itself a random number, then you could randomize that too.
You can replace RandomReal with RandomInteger too.
res = {0, RandomInteger[sum, num - 1], 28} // Flatten // Sort // Differences
(* {1, 2, 23, 1, 1} *)
res//Total
(* 28 *)
This problem has echos of this topic.
Uniformly distributed n-dimensional probability vectors over a simplex
Edit
Here's a histogram of the RandomInteger version, having generated 100,000 examples each with 5 numbers that sum to 28. Generating the samples is almost instantaneous.
Here it is with sum = 28, 5 numbers, and using RandomReal, 1,000,000 runs.
Histogram[
Table[{0, RandomReal[sum, 5 - 1], 28} // Flatten // Sort //
Differences, {1000000}] // Flatten, 100]
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}];
Adapting the code from RandomComposition from Combinatorica` which we can see see by using PrintDefinitions from GeneralUtilities` (by doing PrintDefinitions@RandomComposition), we can define
myRandomComposition[n_Integer, k_Integer] :=
Differences@
Join[{0},
Sort@RandomSample[Range[(n + k) - 1], k - 1], {n + k}] - 1;
I see that this is (pretty much) the riggedRandom function in the answer by DavidC and I almost deleted my answer. The discussion in the comments of that answer is relevant, as the function RandomKSubset used in the "Combinatorica`" function is equivalent to RandomSample, as opposed to RandomChoice. I did not delete my answer because my function is faster because it handles packed arrays better. For a proper comparison, we also test the version with RandomChoice
myRiggedRandom[n_Integer, k_Integer] :=
Differences@
Join[{0},
Sort@RandomChoice[Range[(n + k) - 1], k - 1], {n + k}] - 1;
We then have
nn = 100000;
kk = 90000;
displayer =
Function[Null, # //
RepeatedTiming // {First@#, Length@Last@#, Total@Last@#} &,
HoldAll];
myRandomComposition[nn, kk] // displayer
RandomComposition[nn, kk] // displayer
riggedRandom[nn, kk] // displayer
myRiggedRandom[nn, kk] // displayer
{0.012,90000,100000} {0.027,90000,100000} {0.018,90000,100000} {0.011,90000,100000}
If you want $n$ random values that sum to a total $t$:
f[n_, t_] :=
t Normalize[RandomVariate[UniformDistribution[], n], Total]
f[5,28]
(*
{5.91145, 5.82431, 8.2661, 3.8773, 4.12084}
*)
Check: Total@%
(*
28
*)
As pointed out by @george2079 and @JacobAkkerboom you need a definition as to what you mean by "random" rather than just giving some code.
It appears from your code that you want a set of integers (randomly chosen from 0 to 28: meaning each integer has a 1/29 chance of being selected on any one draw) that sum to 28. There might be just one integer in a "random set" (if 28 is chosen as the first selection) or there could be hundreds if lots of zeros are chosen.
Note that something similar to the above description is what you need to give so that we don't have to guess and you won't get as many answers that don't quite fit. For example, I'm assuming that the length of the random set could be anywhere from 1 to infinity. Maybe you had a stopping rule in mind or that there are a fixed number of draws. The point is that such details are needed.
Here is some code to generate a random sample with the above definition:
n = 28;
x = {};
While[Total[x] < n,
x = Append[x, RandomInteger[{0, n - Total[x]}]]];
Print[x]
(* {24,1,2,0,1} *)
This could be put into a loop to generate multiple samples:
n = 28;
nsamples = 5;
Do[
x = {};
While[Total[x] < n,
x = Append[x, RandomInteger[{0, n - Total[x]}]]];
Print[x],
{i, nsamples}]
{27,0,0,1}
{13,6,1,0,3,2,1,2}
{1,6,14,0,1,0,1,5}
{21,4,1,0,2}
{14,9,5}
This is separate from my other answer and maybe the one that the OP actually requested (because I didn't read carefully enough the first time).
If independent random samples are taken 5 times with for each time the integers 0 through 28 are equally likely (i.e., a 1/29 chance of being selected for each of the 5 selections). The 29^5 = 20,511,149 equally likely possibilities are then screened to just those samples which sum to 28.
The complete sample space of equally likely outcomes where the sum is 28 is generated by
n = 28;
sampleSpace = Select[Tuples[Range[0, n], 5], Total[#] == n &];
This takes about 20 to 30 seconds to generate. Because each outcome in the complete set of sample values is equally likely to occur, any subset (such as the subset that sums to 28) is also equally likely to occur. So we can take a "random sample" in the following manner:
RandomChoice[sampleSpace]
(* {13, 10, 1, 1, 3} *)
This is a bit "brute-force" but it is simple and is constructed based on a known set of probabilities so we know what is meant by a "random sample".
Total[x] == 28? I mean, they are not independent. p.s. Not efficient -RandomChoice @ IntegerPartitions[28, {5}]$\endgroup$While[Total@x=!=28,x = RandomInteger[{0, 28}, 5]];Print[x];(note the print now outside the loop) $\endgroup$