12
$\begingroup$

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?

$\endgroup$
  • 4
    $\begingroup$ 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}] $\endgroup$ – Kuba Jul 10 '14 at 12:24
  • $\begingroup$ 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) $\endgroup$ – george2079 Jul 10 '14 at 12:33
  • 3
    $\begingroup$ note there is a pretty challenging question of what exactly do you mean by random. $\endgroup$ – george2079 Jul 10 '14 at 12:47
  • 1
    $\begingroup$ Related: codegolf.stackexchange.com/questions/8574/… $\endgroup$ – DavidC Jul 10 '14 at 13:12
  • $\begingroup$ @Kuba See my last comment on the answer. $\endgroup$ – Szabolcs Jul 10 '14 at 13:32

11 Answers 11

24
$\begingroup$

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
$\endgroup$
  • $\begingroup$ 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). $\endgroup$ – george2079 Jul 10 '14 at 12:54
  • 2
    $\begingroup$ +1 Nice and clean. You should perhaps wrap the function in RandomSample to avoid necessarily sorting from largest to smallest. $\endgroup$ – DavidC Jul 10 '14 at 13:05
  • 1
    $\begingroup$ @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. $\endgroup$ – Szabolcs Jul 10 '14 at 13:29
  • 1
    $\begingroup$ 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? $\endgroup$ – Szabolcs Jul 10 '14 at 13:31
  • 1
    $\begingroup$ I guess the max entropy distribution makes sense for the OP since it's the "most random" and he wanted random :P $\endgroup$ – Rojo Jul 13 '14 at 2:27
14
$\begingroup$

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}}]

picture

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

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}
 ]

enter image description here

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).

$\endgroup$
  • 2
    $\begingroup$ 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. $\endgroup$ – Rojo Jul 13 '14 at 5:20
  • $\begingroup$ 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. $\endgroup$ – acl Jul 13 '14 at 14:13
5
$\begingroup$

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}
$\endgroup$
5
$\begingroup$

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}].

$\endgroup$
  • 1
    $\begingroup$ Your answer looks nice, though I haven't had the time to read it in full detail. Skimming your answer, it looks like you probably have a stronger grasp on combinatorics than me, but maybe I can still provide some useful comments. First of all "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$ – Jacob Akkerboom Jun 29 '17 at 8:34
  • 1
    $\begingroup$ Most, if not all, the "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$ – Jacob Akkerboom Jun 29 '17 at 8:35
  • 1
    $\begingroup$ Which error messages are you seeing when you load 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$ – Kellen Myers Jun 29 '17 at 15:19
  • $\begingroup$ I thought I saw a "real" error message today, rather than just the warning message and an additional message about a symbol that was shadowed (which does not usually always happen when I load it and is totally reasonable). But I could not reproduce it, maybe it was my own mistake (I thought it said something along the lines of "Set::write: Tag a in a[___List] is Protected."). $\endgroup$ – Jacob Akkerboom Jun 29 '17 at 17:02
4
$\begingroup$

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}];
$\endgroup$
4
$\begingroup$

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.

enter image description here

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]

enter image description here

$\endgroup$
  • 1
    $\begingroup$ Very clever. Any idea what the resulting distribution is? $\endgroup$ – John Doty Jun 29 '17 at 14:29
  • 1
    $\begingroup$ Not sure. The method when applied to populating a simplex results in a uniform distribution over the simplex face. If you then project it down to a single axis, you get this shape. I'll show a histogram of it. $\endgroup$ – MikeY Jun 29 '17 at 14:36
  • $\begingroup$ This approach seems not to echo but to literally ask the same question as the linked topic (Uniformly distributed n-dimensional probability vectors over a simplex) in one dimension, yes? Or is there a subtlety that I've missed in your answer? $\endgroup$ – Kellen Myers Jun 29 '17 at 15:21
  • $\begingroup$ The other topic had the added subtlety of looking for a uniform distribution of results on the face. No such requirement here, and the histogram above shows it's not close. Now requiring them to be uniform would be a tricky complication... $\endgroup$ – MikeY Jun 29 '17 at 16:25
2
$\begingroup$

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}
$\endgroup$
1
$\begingroup$

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}

$\endgroup$
1
$\begingroup$

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".

$\endgroup$
1
$\begingroup$

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

*)

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.