Random Total
June 12, 2015
Our task today is to generate a set of k random positive integers that add up to a given totaln. For instance, if we want 4 random numbers that add up to 9, there are six possible results (not counting permutations of them): {6,1,1,1}, {5,2,1,1}, {4,3,1,1}, {4,2,2,1}, {3,3,2,1} and {3,2,2,2}.
An easy way to do that is to choose k−1 random numbers r with 0 ≤ r < n−k, sort them, calculate the differences between them, calculate the difference between 0 and the smallest, calculate the difference between n−k and the largest, shuffle the differences, and add 1 to each; subtracting k and adding 1 ensures that all the numbers are positive. For our example above, choose three random non-negative integers less than n−k = 5, say 1, 3, and 3, the differences are 1, 2, 0 and 2, and the four resulting numbers are 2, 3, 1 and 3, which form the fifth of the six sets shown above.
Your task is to write the program that generates a random set of integers that adds to a given total, as described above. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.
Programming Praxis 
I may be incorrect (I’m feeling slow due to a cold running through my house),
but I think it can be done with a simple recursion if you don’t care about
order. Details in Scheme follow.
First, here’s a function that returns a list of
ncopies ofx:(define (repeat n x) (let loop ((n n) (acc (list))) (if (< n 1) acc (loop (-1+ n) (cons x acc)))))We’ll use that in a base case of our recursion. Next, the main event:
(define (random-total k n) (cond ((or (< k 1) (> k n)) (error "1 <= k <= n only")) ((< n 1) (error "1 <= n only")) (else (let loop ((k k) (n n) (acc (list))) (cond ((= 1 k) (cons n acc)) ((= k n) (append (repeat k 1) acc)) (else (let ((r (1+ (random (- n k 1))))) (loop (-1+ k) (- n r) (cons r acc)))))))))I’m on a computer without a decent editor, so please pardon any typos or errors in
my counting parentheses. As always, let me know if my “solution” really isn’t one.
@Graham: It worked fine, with no typos. But it’s not very good. Here’s a typical run:
> (random-total 17 1000)
(3 1 1 1 1 1 4 3 10 3 6 6 48 92 29 714 77)
The problem is that once you have a “big” number, like 714 in the example above, all the rest of the numbers have to be “little”, and you end up with 5 ones, 3 threes, 1 four, two sixes, and a ten. The function did indeed produce 17 random numbers that add up to 1000, but it’s probably not what you want.
Don’t feel bad. I did the same thing the first time I solved this problem.
@ProgrammingPraxis: Thanks for responding!
Why would something like this not be what you want, though? Shouldn’t all (permutations of) k-partitions of n be equally likely?
In that case, random totals where all the summands are “close” would have less probability than ones with “big” numbers.
Scratch that, my “solution” looks much more biased in experimentation than I expected.
@Graham: I did 10000 trials of the suggested solution on (rand-total 4 9) and came up with these counts:
> (uniq-c equal? (sort lt? (map (lambda (x) (sort < x)) (map (lambda (x) (rand-total 4 9)) (range 10000))))) (((1 1 1 6) . 88) ((1 1 2 5) . 1180) ((1 1 3 4) . 1125) ((1 2 2 4) . 2868) ((1 2 3 3) . 2883) ((2 2 2 3) . 1856))And here’s the same thing with your version of rand-total:
So there’s no difference. But when n is very much larger than k, the difference between the two methods is obvious. I’m not sure which is “right”, but I think the suggested solution looks “better” based on my expectations of what the program should do.
I guess the right answer is to choose a random number between 1 and n! / (n-k)! — I don’t think I got that right — and map that to a particular permutation. I’ll have to think about that. But I can tell you that the algorithm of the suggested solution is the “computing folklore” algorithm.
@ProgrammingPraxis: cool stuff, and thanks for running the experiments. I need to pick up more of the folklore.
As for the other route you’re thinking of, check out the combinadic number system.
@Praxis: Here’s a paper describing your idea: http://files.figshare.com/484635/Partitioning_paper.pdf
In summary:
1. Calculate p(n,k) = the number of integer partitions of n, having k parts.
2. Generate a random number r, 1 <= r <= p(n,k), then the r-th partition in lex-order is the randomly selected partition.
3. Map r onto a partition:
a. let s = ceil(n / k), the smallest big number that can be in the partition. For p(9, 4), s = 3, corresponding to the partition (3, 2, 2, 2).
b. let t = n – (k – 1), the largest number that can be in the partition. For p(9, 4), t = 6, corresponding to the partition (6, 1, 1, 1).
c. for j in range [s, t]
d. calculate, x = the number of partitions in which the largest number is j.
e. If x = r, you know j is the first digit in the partition.
e. now let n -= j and k -=1 and goto step a to find the next number in the partition
Here’s a solution based on the paper I linked above.
import random from math import ceil pk_cache = {} def pk(n, k): """return number of k-partitions of n. """ if n == 0 == k: return 1 elif n <= 0 or k <= 0: return 0 else: key = (n, k) if key not in pk_cache: pk_cache[key] = pk(n - k, k) + pk(n - 1, k - 1) return pk_cache[key] pkm_cache = {} def pkm(n, k, m): """return number of k-partitions of n, in which the largest element is less than or equal to m. """ if 0 < k <= n <= m*k: if n == m*k: return 1 key = (n, k, m) if key not in pkm_cache: s = 0 for i in range(n//m + 1): nn = n - i*m nk = k - i nm = min(m - 1, nn - nk + 1) s += pkm(nn, nk, nm) pkm_cache[key] = s return pkm_cache[key] return 0 def rand_partition(n, k): """ select a uniformly random partition from the k-partitions of n. """ partition = [] min_m = ceil(n/k) max_m = n - k + 1 which = random.randrange(pk(n, k)) # count the number of k-partitions and select a random index while n: for m in range(min_m, max_m + 1): # determine the value of m for which pkm(n,k,m) contains the index count = pkm(n, k, m) # to determine the largest element of the partition if which < count: count = pkm(n, k, m - 1) break partition.append(m) n -= m # repeat for the next rest of the partition if n == 0: break k -= 1 which -= count min_m = (n//k + 1) if n % k else (n//k) max_m = m return partitionTest:
from collections import Counter part = Counter(tuple(rand_partition(9,4)) for _ in range(12000)) for k,v in sorted(part.items()): print(k, ':', v) (3, 2, 2, 2) : 1994 (3, 3, 2, 1) : 1952 (4, 2, 2, 1) : 2004 (4, 3, 1, 1) : 1976 (5, 2, 1, 1) : 2098 (6, 1, 1, 1) : 1976My python-fu is still very lacking, but hopefully I got in the ballpark of the right answer. Long time follower, first time contributor, so hopefully I get the formatting and whatnot correct.
def random_total(desired_total, list_size): if desired_total < list_size: raise Exception("list_size must be less than desired_total") if( list_size < 1): raise Exception("list size must be greater than 0") random_numbers = sorted([randint(0, desired_total-list_size) for x in range(list_size-1)], reverse=True) # Choose the random numbers, sort in descending order. differences = [random_numbers[i] - random_numbers[i+1] for i in range(len(random_numbers)-1)] # Find the differences between neighbor values. differences.append(random_numbers[-1]) # Find the difference of the smallest and 0 differences.append((desired_total-list_size) - max(random_numbers)) # Find the difference of the n-k and the biggest value. shuffled_differences = [x + 1 for x in differences] #Shuffle the values and add one return shuffled_differencesTo test this, I just ran the algorithm a bunch of times on (9,4), assuming that if run enough times I could get all possible combinations, and printed those to screen. If they matched, at least as a litmus test it should be approximately correct.
for x in range(100): simulation_result = sorted(random_total(9,4), reverse=True) if simulation_result not in outcomes: outcomes.append(simulation_result) for values in sorted(outcomes, reverse=True): print(",".join(map(str, values)))Output:
6,1,1,1
5,2,1,1
4,3,1,1
4,2,2,1
3,3,2,1
3,2,2,2
Here’s a run with a larger n:
rand_partition(1000, 17) => [281, 133, 85, 83, 71, 66, 66, 45, 40, 34, 31, 26, 15, 11, 8, 4, 1]
Of interest, there are pk(1000, 17) ~ 3.3 x 1020 partitions. The python docs say the random number generator only has about 53 bits of precision, whereas approx 69 bits of precision are needed to index all 3.3 x 1020 partitions . As written, the function above would not be able to choose some of the partitions.
I don’t think the folklore solution generates numbers with the proper distribution. For example, looking at the case for rand_partition(9,4), here is a comparison odds for the largest number in a partition:
largest #true probfolklore prob
60.1670.006
50.1670.124
40.3330.405
30.3330.467
table tags didn’t work
largest true folklore
number prob prob
6 0.167 0.006
5 0.167 0.124
4 0.333 0.405
3 0.333 0.467
@Mike: No, the folklore solution doesn’t generate numbers with the proper distribution; you saw that in my experiment with Graham. But it does generate solutions that “look good” in a way similar to your solution based on a random partition.
If you need random numbers with lots of bits, use a random number generator to generate as many coin flips as necessary, one at a time.
Here’s another way: a k-partition of n either has smallest member 1, or all elements are > 1. If the smallest member is 1, we can remove this, and get a k-1 partition of n-1; if all elements are > 1, we can subtract 1 from each element and get a k-partition of n-k. If p(n,k) is the number of partitions for a given n and k, this gives a recurrence relation:
p(n,k) = p(n-1,k-1) + p(n-k,k)
The base cases are p(0,0) = 1, otherwise p(n,k) = 0 if n <= 0 or k <= 0.
This enables us to both calculate the number of partitions and also map an integer i, 0 <= i < p(n,k) to a partition. Here's some Haskell:
p n k
| k == 0 && n == 0 = 1
| k <= 0 || n <= 0 = 0
| otherwise = p (n-k) k + p (n-1) (k-1)
gen n k i
| k <= 0 || n <= 0 = []
| otherwise =
if i < t then 1:gen(n-1)(k-1) i
else map (+1) (gen(n-k) k (i-t))
where t = p(n-1)(k-1)
parts n k = map (gen n k) [0..(p n k)-1]
main = print (parts 10 5)
— [[1,1,1,1,6],[1,1,1,2,5],[1,1,1,3,4],[1,1,2,2,4],[1,1,2,3,3],[1,2,2,2,3],[2,2,2,2,2]]
p would benefit from memoization, and I've left out actually selecting a random number in the appropriate range.
p n k | k == 0 && n == 0 = 1 | k <= 0 || n <= 0 = 0 | otherwise = p (n-k) k + p (n-1) (k-1) gen n k i | k <= 0 || n <= 0 = [] | otherwise = if i < t then 1:gen(n-1)(k-1) i else map (+1) (gen(n-k) k (i-t)) where t = p(n-1)(k-1) parts n k = map (gen n k) [0..(p n k)-1] main = print (parts 10 5) -- [[1,1,1,1,6],[1,1,1,2,5],[1,1,1,3,4],[1,1,2,2,4],[1,1,2,3,3],[1,2,2,2,3],[2,2,2,2,2]]And we can turn this in to a nice iterative solution – rather than building up the partition directly, we generate the differences between the elements of the partition. C++ for this one:
#include <vector> #include <stdlib.h> #include <stdio.h> int p(int n, int k) { if (n == 0 && k == 0) return 1; if (k <= 0 || n <= 0) return 0; return p(n-k,k) + p(n-1,k-1); } void gen(int *a, int n, int k, int i) { int j = 0; while (j < k && n > 0) { int t = p(n-1,k-j-1); if (i < t) { n--; j++; } else { a[j]++; n-=k-j; i-=t; } } } int main(int argc, char *argv[]) { int n = atoi(argv[1]); int k = atoi(argv[2]); for (int i = 0; i < p(n,k); i++) { std::vector<int>a(k); gen(&a[0],n,k,i); int t = 1; for (int j = 0; j < k; j++) { t += a[j]; printf("%d%s", t, (j<k-1)?" ":"\n"); } } }@matthew — nice!
public class RandomNumbers { public static void main(String[] args){ System.out.println("Enter the number of random positive integers: "); Scanner scanner = new Scanner(System.in); int k = scanner.nextInt(); System.out.println("Enter total to be scored"); int n = scanner.nextInt(); List<Integer> numList = new ArrayList<Integer>(); int[] input = new int[k-1]; Random t = new Random(); int count = k -1 ; for(int i = 0; i<input.length; i++){ input[i] = t.nextInt(n-k); } Arrays.sort(input); numList.add(input[0]-0); for(int i = 1; i<=input.length-1;i++){ numList.add(input[i] - input[i-1]); } numList.add((n-k)-input[input.length-1]); int[] finallist = new int[numList.size()]; for(int i = 0;i<finallist.length; i++){ finallist[i] = numList.get(i).intValue(); finallist[i] = finallist[i]+1; } for(int i =0;i<finallist.length;i++){ System.out.print(" "+finallist[i]); } } }