# Improve the method to calculate the frequency distribution of sums for n 6-sided fair dice

I was tring to excute the code below for $n\geq 10$, but my computer somehow aborts automatically. Is there a way to improve the code?

(#/Total@Values@#) &[Counts@Flatten@
Array[
Evaluate[
Plus @@ Table[ToExpression["#" <> ToString[i]], {i, n}]] &,
Table[6, n]]]

• Monitor your memory usage as you run this using an external program, like Activity Monitor on mac os. Mine spiked above 2.5 GB, so it is very likely that is what is causing the abort. I have not taken the code apart to figure out what it is doing, so I cannot suggest alternatives at the moment. – rcollyer Nov 16 '16 at 4:21
• I have version 7, so I can't try out your code. Can you explain in words what it does? – Kiro Nov 16 '16 at 5:44
• @rcollyer I tried when $n=11$ and the maximum memory usage was about 5 GB… Also the timing command returned 0.000107 s but it was at least 30 s . Maybe it's memory leak or something? Before $n\geq 9$ all the results returned almost instantaneously. – Taptic Nov 16 '16 at 6:56
• Just to make sure I understood you correctly, does it more or less do this? ListPlot[ Tally[Sort[Total /@ Tuples[Range[6], n]]], Filling -> Axis, PlotMarkers -> {Automatic, Medium}] – Kiro Nov 16 '16 at 7:16
• @Kiro Yep. And now I know why it took so long… – Taptic Nov 16 '16 at 7:23

From your comment you seem to have re-discovered probability generating functions.

For $n=1$, there is an equal probability (1/6) to get a 1, 2, 3, 4, 5, or a 6. The probability generating function for $n=1$ is

pgf = Sum[t^i/6, {i, 6}]


For $n>1$, the probability generating function of the sum is pgf^n. The probability of the sum being some value s is given by the coefficient of t^s in pgf^n. So the probability of the sum being 12 with n=11 is given by

Coefficient[pgf^11, t^12]
(* 11/362797056 *)


Using CoefficientList[pgf^n,t] gets you a list of the probabilities of the sum being 0 to 6n (from which you can easily remove the first n items in the list as those probabilities are zero). The frequencies can be obtained by multiplying the probabilities by 6^n.

I will go ahead and extend my comment into an answer. There are $6^n$ different combinations of results from throwing $n$ 6-sided fair dice, i.e. the number of cases grows explosively. To get all combinations, one can simply use the built-in function Tuples which does just that, given correct input.

Since there are considerably less different total sums of these, one option would be to try to find an algorithm that can find the sums without considering all $6^n$ combinations. If you have or find one and need help implementing that in Mathematica, you can edit your question or perhaps ask a new one.

EDIT: I gave this a little bit more thought and came up with what should be a good starting point to arrive at a full answer. Consider the code below.

Table[
IntegerPartitions[i, {n}, Range[6]],
{i, n, 6 n}]


This generates all possible unordered combinations (there will be a great deal less of them than $6^n$). To get the frequencies, the output would have to be processed a bit further, say, by replacing each element with a list of two elements: one for the sum and the other for the number of different combinations of the elements, with repeated elements treated as identical. The latter task should be possible to do without actually calling Permutations, by just looking at how many distinct elements there are.

The brute-force solution using Tuples could be used to check the answer for small n.

• Apparently the coefficients are related to Sextinomial link. Manipulate[ ListPlot@Transpose@{Range[n, 6 n], CoefficientList[(1 + x + x^2 + x^3 + x^4 + x^5)^n, x]}, {n, 1, 100, 1}] works like a charm – Taptic Nov 16 '16 at 8:28
• Nicely done! You might want to answer your own question if you want to make it 'official', I for one would like to upvote that. – Kiro Nov 16 '16 at 9:01
• I'll modify the question so that the answer is related to it. – Taptic Nov 17 '16 at 20:17

Here are two ways calculate the frequency distribution of the sums for $n$ six-sided fair dice:

1. To calculate the frequency of each possible sum (or a combination of $n$ sides), each dice can be thought of as a polynomial $$p(x)=\left(x^1+x^2+x^3+x^4+x^5+x^6\right)\\\text{(the powers represent the numbers on different sides)}$$ Also, $n$ dice can be represented by $$P(x)=\left(x^1+x^2+x^3+x^4+x^5+x^6\right)^n$$ The frequency of a particular sum say $m$ is equal to the sum of the coefficients of the terms that are divisible by $\prod_{i=1}^6(x^i)^{k_i}$ after $P(x)$ is expanded, with the condition that $\sum_{i=1}^6 k_i=n$ and $\sum_{i=1}^6 i k_i=m$. $k_i$ represents the number of times that side $i$ appears in $n$ dice. The coefficient can be written as $$\sum_{(\sum_{i=1}^6 k_i=n)\wedge(\sum_{i=1}^6 i k_i=m)} \frac{n!}{\prod_{i=1}^6 k_i !}=\sum_{(\sum_{i=1}^6 k_i=n)\wedge(\sum_{i=1}^6 i k_i=m)}\binom{n}{k_1,\ldots,k_6}$$ The expression on the RHS is the multinomial coefficient.

Even though the above expression looks complicated, notice that $\prod_{i=1}^6(x^i)^{k_i}=x^m$, therefore we only need to find the coefficient of $m$th term when the expansion of $P(x)$ is simplified.

The code for Mathematica is
Transpose@ {Range[n, 6 n], CoefficientList[(1 + x + x^2 + x^3 + x^4 + x^5)^n, x]}.
(the coefficients are the same, and starting with $x^0=1$ has the advantage that there's no need to drop first $n$ elements)
2. Thanks to Kiro's answer that relates partition of integers with combinations, there's another way to calculate frequency that is closely related to the first method.

Partitions of an integer $i$ with exactly $n$ numbers are that ranging from 1 to 6 is given by IntegerPartitions[i, {n}, Range[6]], this can be thought of as the ways to create a number $i$ with $n$ 6-sided fair dice. However the function only gives unordered combinations. For example, the number of ways to create number 5 with 4 dice is 4 ((2,1,1,1), (1,2,1,1), (1,1,2,1), (1,1,1,2)), but IntegerPartitions only gives 1 result ((2,1,1,1)).

To fix the problem we use
Table[{i, Total[ (Multinomial @@ Last@Transpose@Tally@#) & /@ IntegerPartitions[i, {n}, Range[6]] ]}, {i, n, 6 n}]
It first gives the tally of a particular partition, then it isolates the part that gives the information about how the numbers are distributed. (This essentially give a particular combination of $k_i$ in method 1) Finally it uses Multinomial to give the frequency that one particular partition will contribute, and it uses Total to find the final number.

This method is significantly slower than the first one, but it provides some insights to how should one think about this problem.