How to easily determine the results distribution for multiple dice?

I want to calculate the probability distribution for the total of a combination of dice.

I remember that the probability of is the number of combinations that total that number over the total number of combinations (assuming the dice have a uniform distribution).

What are the formulas for

• The number of combinations total
• The number of combinations that total a certain number
-

Exact solutions

The number of combinations in $n$ throws is of course $6^n$.

These calculations are most readily done using the probability generating function for one die,

$$p(x) = x + x^2 + x^3 + x^4 + x^5 + x^6 = x \frac{1-x^6}{1-x}.$$

(Actually this is $6$ times the pgf--I'll take care of the factor of $6$ at the end.)

The pgf for $n$ rolls is $p(x)^n$. We can calculate this fairly directly--it's not a closed form but it's a useful one--using the Binomial Theorem:

$$p(x)^n = x^n (1 - x^6)^n (1 - x)^{-n}$$

$$= x^n (\sum_{k=0}^{n} {n \choose k} (-1)^k x^{6k} )( \sum_{j=0}^{\infty} {-n \choose j} (-1)^j x^j).$$

The number of ways to obtain a sum equal to $m$ on the dice is the coefficient of $x^m$ in this product, which we can isolate as

$$\sum_{6k + j = m - n} {n \choose k}{-n \choose j}(-1)^{k+j}.$$

The sum is over all nonnegative $k$ and $j$ for which $6k + j = m - n$; it therefore is finite and has only about $(m-n)/6$ terms. For example, the number of ways to total $m = 14$ in $n = 3$ throws is a sum of just two terms, because $11 = 14-3$ can be written only as $6 \cdot 0 + 11$ and $6 \cdot 1 + 5$:

$$-{3 \choose 0} {-3 \choose 11} + {3 \choose 1}{-3 \choose 5}$$

$$= 1 \frac{(-3)(-4)\cdots(-13)}{11!} + 3 \frac{(-3)(-4)\cdots(-7)}{5!}$$

$$= \frac{1}{2} 12 \cdot 13 - \frac{3}{2} 6 \cdot 7 = 15.$$

(You can also be clever and note that the answer will be the same for $m = 7$ by the symmetry 1 <--> 6, 2 <--> 5, and 3 <--> 4 and there's only one way to expand $7 - 3$ as $6 k + j$; namely, with $k = 0$ and $j = 4$, giving

$${3 \choose 0}{-3 \choose 4} = 15 \text{.)}$$

The probability therefore equals $15/6^3$ = $5/36$, about 14%.

By the time this gets painful, the Central Limit Theorem provides good approximations (at least to the central terms where $m$ is between $\frac{7 n}{2} - 3 \sqrt{n}$ and $\frac{7 n}{2} + 3 \sqrt{n}$: on a relative basis, the approximations it affords for the tail values get worse and worse as $n$ grows large).

I see that this formula is given in the Wikipedia article Srikant references but no justification is supplied nor are examples given. If perchance this approach looks too abstract, fire up your favorite computer algebra system and ask it to expand the $n^{\text{th}}$ power of $x + x^2 + \cdots + x^6$: you can read the whole set of values right off. E.g., a Mathematica one-liner is

With[{n=3}, CoefficientList[Expand[(x + x^2 + x^3 + x^4 + x^5 + x^6)^n], x]]

-
Will that mathematica code work with wolfram alpha? –  user28 Oct 14 '10 at 23:56
@Srikant I didn't even try it because usually alpha does not process Mathematica code correctly except for simple expressions. However, see what happens when you type "Expand[(x + x^2 + x^3 + x^4 + x^5 + x^6)^3]" ! ( wolframalpha.com ) –  whuber Oct 15 '10 at 14:15
That works. I tried your earlier version but could not make any sense of the output. –  user28 Oct 15 '10 at 14:24
@Srikant: Expand[Sum[x^i,{i,1,6}]^3] also works in WolframAlpha –  A. N. Other Oct 18 '10 at 7:04

Yet another way to quickly compute the probability distribution of a dice roll would be to use a specialized calculator designed just for that purpose.

Torben Mogensen, a CS professor at DIKU has an excellent dice roller called Troll.

The Troll dice roller and probability calculator prints out the probability distribution (pmf, histogram, and optionally cdf or ccdf), mean, spread, and mean deviation for a variety of complicated dice roll mechanisms. Here are a few examples that show off Troll's dice roll language:

Roll 3 6-sided dice and sum them: sum 3d6.

Roll 4 6-sided dice, keep the highest 3 and sum them: sum largest 3 4d6.

Roll an "exploding" 6-sided die (i.e., any time a "6" comes up, add 6 to your total and roll again): sum (accumulate y:=d6 while y=6).

Troll's SML source code is available, if you want to see how its implemented.

Professor Morgensen also has a 29-page paper, "Dice Rolling Mechanisms in RPGs," in which he discusses many of the dice rolling mechanisms implemented by Troll and some of the mathematics behind them.

A similar piece of free, open-source software is Dicelab, which works on both Linux and Windows.

-

Approximate Solution

I explained the exact solution earlier (see below). I will now offer an approximate solution which may suit your needs better.

Let:

$X_i$ be the outcome of a roll of a $s$ faced dice where $i=1, ... n$.

$S$ be the total of all $n$ dice.

$\bar{X}$ be the sample average.

By definition, we have:

$\bar{X} = \frac{\sum_iX_i}{n}$

In other words,

$\bar{X} = \frac{S}{n}$

The idea now is to visualize the process of observing ${X_i}$ as the outcome of throwing the same dice $n$ times instead of as outcome of throwing $n$ dice. Thus, we can invoke the central limit theorem (ignoring technicalities associated with going from discrete distribution to continuous), we have as $n \rightarrow \infty$:

$\bar{X} \sim N(\mu, \sigma^2/n)$

where,

$\mu = (s+1)/2$ is the mean of the roll of a single dice and

$\sigma^2 = (s^2-1)/12$ is the associated variance.

The above is obviously an approximation as the underlying distribution $X_i$ has discrete support.

But,

$S = n \bar{X}$.

Thus, we have:

$S \sim N(n \mu, n \sigma^2)$.

Exact Solution

The wiki has a brief explanation as how to calculate the required probabilities. I will elaborate a bit more as to why the wiki's explanation makes sense. To the extent possible I have used similar notation as the wiki.

Suppose that you have $n$ dice each with $s$ faces and you want to compute the probability that a single roll of all $n$ dice the total adds up to $k$. The approach is as follows:

Define:

$F_{s,n}(k)$: Probability that you get a total of $k$ on a single roll of $n$ dices with $s$ faces.

By definition, we have:

$F_{s,1}(k) = \frac{1}{s}$

The above states that if you just have one dice with $s$ faces the probability of obtaining a total $k$ between 1 and s is the familiar $\frac{1}{s}$.

Consider the situation when you roll two dice: You can obtain a sum of $k$ as follows: The first roll is between 1 to $k-1$ and the corresponding roll for the second one is between $k-1$ to $1$. Thus, we have:

$F_{s,2}(k) = \sum_{i=1}^{i=k-1}{F_{s,1}(i) F_{s,1}(k-i)}$

Now consider a roll of three dice: You can get a sum of $k$ if you roll a 1 to $k-2$ on the first dice and the sum on the remaining two dice is between $k-1$ to $2$. Thus,

$F_{s,3}(k) = \sum_{i=1}^{i=k-2}{F_{s,1}(i) F_{s,2}(k-i)}$

Continuing the above logic, we get the recursion equation:

$F_{s,n}(k) = \sum_{i=1}^{i=k-n+1}{F_{s,1}(i) F_{s,n-1}(k-i)}$

See the wiki link for more details.

-
@Srikant Excellent answer, but does that function resolve to something arithmetic (ie: not recursive)? –  C. Ross Oct 14 '10 at 20:06
@C. Ross Unfortunately I do not think so. But, I suspect that the recursion should not be that hard as long as are dealing with reasonably small n and small s. You could just build-up a lookup table and use that repeatedly as needed. –  user28 Oct 14 '10 at 20:10
The wikipedia page you linked has a simple nonrecursive formula which is a single sum. One derivation is in whuber's answer. –  Douglas Zare Jul 17 '12 at 23:44
The wiki link anchor is dead, do you know of a replacement? –  Midnighter Apr 13 at 20:21

Characteristic functions can make computations involving the sums and differences of random variables really easy. Mathematica has lots of functions to work with statistical distributions, including a builtin to transform a distribution into its characteristic function.

I'd like to illustrate this with two concrete examples: (1) Suppose you wanted to determine the results of rolling a collection of dice with differing numbers of sides, e.g., roll two six-sided dice plus one eight-sided die (i.e., 2d6+d8)? Or (2) suppose you wanted to find the difference of two dice rolls (e.g., d6-d6)?

An easy way to do this would be to use the characteristic functions of the underlying discrete uniform distributions. If a random variable $X$ has a probability mass function $f$, then its characteristic function $\varphi_X(t)$ is just the discrete Fourier Transform of $f$, i.e., $\varphi_X(t) = \mathcal{F}\{f\}(t) = E[e^{i t X}]$. A theorem tells us:

If the independent random variables $X$ and $Y$ have corresponding probability mass functions $f$ and $g$, then the pmf $h$ of the sum $X + Y$ of these RVs is the convolution of their pmfs $h(n) = (f \ast g)(n) = \sum_{m=-\infty}^\infty f(m) g(n-m)$.

We can use the convolution property of Fourier Transforms to restate this more simply in terms of characteristic functions:

The characteristic function $\varphi_{X+Y}(t)$ of the sum of independent random variables $X$ and $Y$ equals the product of their characteristic functions $\varphi_{X}(t) \varphi_{Y}(t)$.

This Mathematica function will make the characteristic function for an s-sided die:

MakeCf[s_] :=
Module[{Cf},
Cf := CharacteristicFunction[DiscreteUniformDistribution[{1, s}],
t];
Cf]


The pmf of a distribution can be recovered from its characteristic function, because Fourier Transforms are invertible. Here is the Mathematica code to do it:

RecoverPmf[Cf_] :=
Module[{F},
F[y_] := SeriesCoefficient[Cf /. t -> -I*Log[x], {x, 0, y}];
F]


Continuing our example, let F be the pmf that results from 2d6+d8.

 F := RecoverPmf[MakeCf[6]^2 MakeCf[8]] 

There are $6^2 \cdot 8 = 288$ outcomes. The domain of support of F is $S=\{3,\ldots,20\}$. Three is the min because you're rolling three dice. And twenty is the max because $20 = 2 \cdot 6 + 8$. If you want to see the image of F, compute

In:= F /@ Range[3, 20]

Out= {1/288, 1/96, 1/48, 5/144, 5/96, 7/96, 13/144, 5/48, 1/9, 1/9, \
5/48, 13/144, 7/96, 5/96, 5/144, 1/48, 1/96, 1/288}


If you want to know the number of outcomes that sum to 10, compute

In:= 6^2 8 F[10]

Out= 30


If the independent random variables $X$ and $Y$ have corresponding probability mass functions $f$ and $g$, then the pmf $h$ of the difference $X - Y$ of these RVs is the cross-correlation of their pmfs $h(n) = (f \star g)(n) = \sum_{m=-\infty}^\infty f(m) g(n+m)$.

We can use the cross-correlation property of Fourier Transforms to restate this more simply in terms of characteristic functions:

The characteristic function $\varphi_{X-Y}(t)$ of the difference of two independent random variables ${X,Y}$ equals the product of the characteristic function $\varphi_{X}(t)$ and $\varphi_{Y}(-t)$ (N.B. the negative sign in front of the variable t in the second characteristic function).

So, using Mathematica to find the pmf G of d6-d6:

 G := RecoverPmf[MakeCf[6] (MakeCf[6] /. t -> -t)] 

There are $6^2 = 36$ outcomes. The domain of support of G is $S=\{-5,\ldots,5\}$. -5 is the min because $-5=1-6$. And 5 is the max because $6-1=5$. If you want to see the image of G, compute

In:= G /@ Range[-5, 5]

Out= {1/36, 1/18, 1/12, 1/9, 5/36, 1/6, 5/36, 1/9, 1/12, 1/18, 1/36}

-
Of course, for discrete distributions, including distributions of finite support (like those in question here), the cf is just the probability generating function evaluated at x = exp(i t), making it a more complicated way of encoding the same information. –  whuber Oct 17 '10 at 21:24
@whuber: As you say, the cf, mgf, and pgf are more-or-less the same and easily transformable into one another, however Mathematica has a cf builtin that works with all the probability distributions it knows about, whereas it doesn't have a pgf builtin. This makes the Mathematica code for working with sums (and differences) of dice using cfs particularly elegant to construct, regardless of the complexity of dice expression as I hope I demonstrated above. Plus, it doesn't hurt to know how cfs, FTs, convolutions, and cross-correlations can help solve problems like this. –  A. N. Other Oct 18 '10 at 2:38
@Elisha: Good points, all of them. I guess what I wonder about the most is whether your ten or so lines of Mathematica code are really more "elegant" or efficient than the single line I proposed earlier (or the even shorter line Srikant fed to Wolfram Alpha). I suspect the internal manipulations with characteristic functions are more arduous than the simple convolutions needed to multiply polynomials. Certainly the latter are easier to implement in most other software environments, as Glen_b's answer indicates. The advantage of your approach is its greater generality. –  whuber Oct 18 '10 at 3:35

Here's another way to calculate the probability distribution of the sum of two dice by hand using convolutions.

To keep the example really simple, we're going to calculate the probability distribution of the sum of a three-sided die (d3) whose random variable we will call X and a two-sided die (d2) whose random variable we'll call Y.

You're going to make a table. Across the top row, write the probability distribution of X (outcomes of rolling a fair d3). Down the left column, write the probability distribution of Y (outcomes of rolling a fair d2).

You're going to construct the outer product of the top row of probabilities with the left column of probabilities. For example, the lower-right cell will be the product of Pr[X=3]=1/3 times Pr[Y=2]=1/2 as shown in the accompanying figure. In our simplistic example, all the cells equal 1/6.

Next, you're going to sum along the oblique lines of the outer-product matrix as shown in the accompanying diagram. Each oblique line passes through one-or-more cells which I've colored the same: The top line passes through one blue cell, the next line passes through two red cells, and so on.

Each of the sums along the obliques represents a probability in the resulting distribution. For example, the sum of the red cells equals the probability of the two dice summing to 3. These probabilities are shown down the right side of the accompanying diagram.

This technique can be used with any two discrete distributions with finite support. And you can apply it iteratively. For example, if you want to know the distribution of three six-sided dice (3d6), you can first calculate 2d6=d6+d6; then 3d6=d6+2d6.

There is a free (but closed license) programming language called J. Its an array-based language with its roots in APL. It has builtin operators to perform outer products and sums along the obliques in matrices, making the technique I illustrated quite simple to implement.

In the following J code, I define two verbs. First the verb d constructs an array representing the pmf of an s-sided die. For example, d 6 is the pmf of a 6-sided die. Second, the verb conv finds the outer product of two arrays and sums along the oblique lines. So conv~ d 6 prints out the pmf of 2d6:

d=:$% conv=:+//.@(*/) |:(2+i.11),:conv~d 6 2 0.0277778 3 0.0555556 4 0.0833333 5 0.111111 6 0.138889 7 0.166667 8 0.138889 9 0.111111 10 0.0833333 11 0.0555556 12 0.0277778  As you can see, J is cryptic, but terse. - +1 for the J example. – whuber Oct 26 '10 at 17:26 There's a very neat way of computing the combinations or probabilities in a spreadsheet (such as excel) that computes the convolutions directly. I'll do it in terms of probabilities and illustrate it for six sided dice but you can do it for dice with any number of sides (including adding different ones). (btw it's also easy in something like R or matlab that will do convolutions) Start with a clean sheet, in a few columns, and move down a bunch of rows from the top (more than 6). 1. put the value 1 in a cell. That's the probabilities associated with 0 dice. put a 0 to its left; that's the value column - continue down from there with 1,2,3 down as far as you need. 2. move one column to the right and down a row from the '1'. enter the formula "=sum(" then left-arrow up-arrow (to highlight the cell with 1 in it), hit ":" (to start entering a range) and then up-arrow 5 times, followed by ")/6" and press Enter - so you end up with a formula like =sum(c4:c9)/6 (where here C9 is the cell with the 1 in it). Then copy the formula and paste it to the 5 cells below it. They should each contain 0.16667 (ish). Don't type anything into the empty cells these formulas refer to! 3. move down 1 and to the right 1 from the top of that column of values and paste ... ... a total of another 11 values. These will be the probabilities for two dice. It doesn't matter if you paste a few too many, you'll just get zeroes. 4. repeat step 3 for the next column for three dice, and again for four, five, etc dice. We see here that the probability of rolling$12$on 4d6 is 0.096451 (if you multiply by$4^6$you'll be able to write it as an exact fraction). If you're adept with Excel - things like copying a formula from a cell and pasting into many cells in a column, you can generate all tables up to say 10d6 in about a minute or so (possibly faster if you've done it a few times). If you want combination counts instead of probabilities, don't divide by 6. If you want dice with different numbers of faces, you can sum$k$(rather than 6) cells and then divide by$k\$. You can mix dice across columns (e.g. do a column for d6 and one for d8 to get the probability function for d6+d8):

-