# How can I efficiently calculate exact probability distributions for D&D dice like "4d8"?

I play a lot of tabletop role-playing games, and I have been interested recently in modeling the various game mechanics involving dice.

As I see it, a single throw of an eight-sided die can be modeled trivially:

d = DiscreteUniformDistribution[{1, 8}]


However, what I'm more often modeling is something like a 4d8, i.e. the throw of four eight-sided dice. If the dice were continuous, this could be modeled as UniformSumDistribution[4, {1, 8}], but of course they are not, and there is no built-in DiscreteUniformSumDistribution.

Instead, I can use a TransformedDistribution:

d = TransformedDistribution[
x1 + x2 + x3 + x4,
{x1, x2, x3, x4},
DiscreteUniformDistribution[{1, 8}]]]


This works, but it's rather slow, especially when the number of dice is large. Are there better existing methods to model such dice throws in Mathematica that I'm not aware of?

• You might want to consider an approach based on generating polynomials. E.g. 1d8 is $D=\sum_{i=1}^{8}x^i$ and 4d8 is then modeled by $D^4$. At the very least, the polynomials end up being a fairly workable size and relative probability information is fairly easy to extract. Jun 2 at 16:42
• From your comments on the answers below, it sounds like you're more interested in computing the probability distribution than actually generating a set of $n$ variates. I've modified the title to reflect this, but feel free to roll it back if it's not what you intended. Jun 2 at 17:49
• FYI, there are tools other than Mathematica, such as AnyDice online. anydice.com/articles/4d6-drop-lowest shows how to have it make plots for stuff like 4d6 drop lowest, or chance of rolling at least x with advantage/disadvantage on a d20. Another stack-exchange site, rpg.stackexchange.com/questions/tagged/anydice has a tag for it, and for [dice] in general. If you want to do further stuff with the distribution inside Mathematica, though, obviously you'd want to use the answers to this question. But if not, there are good tools designed specifically for dice probability. Jun 3 at 2:14

Based on eyorble's comment here's how to do it with generating functions. The main idea: with some probability distribution $$p_n: P(X=n)=p_n$$ you can associate a polynomial called the generating function, which is defined by $$G(x)=\sum_np_nx^n$$ The nice thing about this object is that when you multiply two generating functions $$G_1,G_2$$ corresponding to two distribution $$X_1,X_2$$, you will get a new generating function which corresponds to the distribution associated with $$X_1+X_2$$. This is because the distribution of the sum collects all probabilities which have the same sum $$P(X_1+X_2=n)=\sum_{i+j=m}p_ip_j$$ Multiplying out two polynomials produces the same result, you should check it out for yourself.

In Mathematica you can do this as follows:

G[x_, d_] := Sum[1/d x^n, {n, 1, d}]
dice = {8, 3, 2};
ndice = Length@dice;
Gtotal = Product[G[x, d], {d, dice}];
prob = CoefficientList[Series[Gtotal, {x, 0, Total@dice}], x];
TableForm[Table[{i, prob[[i + 1]]}, {i, ndice, Total@dice}],


$$\begin{array}{cc} \text{Total} & \text{P(Total)} \\\hline 3 & \frac{1}{48} \\ 4 & \frac{1}{16} \\ 5 & \frac{5}{48} \\ 6 & \frac{1}{8} \\ 7 & \frac{1}{8} \\ 8 & \frac{1}{8} \\ 9 & \frac{1}{8} \\ 10 & \frac{1}{8} \\ 11 & \frac{5}{48} \\ 12 & \frac{1}{16} \\ 13 & \frac{1}{48} \\ \end{array}$$ I tried this with 20 d20's and generating the final list took about 1 ms so this should be plenty fast.

• I note that this is slightly faster to do with ListConvolve[]: Drop[Fold[ListConvolve[#1, #2, {1, -1}, 0] &, Table[Prepend[ConstantArray[1/d, d], 0], {d, dice}]], ndice]. Jun 4 at 12:49

Implementing @eyorble's idea in the comments:

problist[n_, d_] := Rest[CoefficientList[(1/d Sum[x^i, {i, 1, d}])^n, x]]


This returns a list of length $$nd$$ containing the exact probabilities of rolling all possible values from 1 to $$nd$$. The entry in the $$m$$th position in this list is the probability of rolling $$m$$ with these dice.

This takes about 3.5 seconds on my computer to calculate all the probabilities of a 1000d8 roll.

Define two helper functions. One for the generator function of NdS expressed as a polynomial in terms of $$x$$:

dg[x_, n_, s_] := (x (x^s - 1)/(x - 1))^n;


And another to feed that into EmpiricalDistribution:

empiricalDistFromGenerator[g_, x_] :=
EmpiricalDistribution[
Transpose[
DeleteCases[
Transpose[{#/Total[#], Range[Length[#]] - 1} &[
CoefficientList[g, x]]], {0, _}]] /. {a_List,
b_List} :> (a -> b)];


Usage:

d6 = empiricalDistFromGenerator[dg[x,1,6],x];
{InverseCDF[d6,0], InverseCDF[d6,1], Mean[d6], StandardDeviation[d6]}


{1, 6, 7/2, Sqrt[35/3]/2}

n3d6 = empiricalDistFromGenerator[dg[x,3,6],x];
{InverseCDF[n3d6,0], InverseCDF[n3d6,1], Mean[n3d6], StandardDeviation[n3d6]}


{3, 18, 21/2, Sqrt[35]/2}

And to provide some timing info on an unreasonable case:

AbsoluteTiming[nbig = empiricalDistFromGenerator[dg[x, 200, 850], x];]
{InverseCDF[nbig, 0], InverseCDF[nbig, 1], Mean[nbig],
StandardDeviation[nbig]}


{32.7313, Null}

{200, 170000, 85100, 5 Sqrt[481666]}

• Fixed d to dg as necessary, not sure I see anywhere that needs :> though. Jun 2 at 17:58
• {a_List, b_List} -> (a -> b) Jun 2 at 18:01
• It works fine without it though, at least for me. I don't see a need for RuleDelayed's properties here, other than eliminating the one set of parentheses? Jun 2 at 18:11
• Yes, it works fine if you don't have any definition prior to the replacement, but if you have for example a=2 before replacement, you'll see the difference. With regard to dozens of people who may read your solution, I would recommend something that minimizes edge cases (just a recommendation). Jun 2 at 18:18
• Here's a slightly streamlined implementation of your idea: With[{n = 3, s = 6}, edist = EmpiricalDistribution[Normalize[Drop[Fold[ListConvolve[#1, #2, {1, -1}, 0] &, PadLeft[ConstantArray[1, {n, s}], {Automatic, s + 1}]], n], Total] -> Range[n, n s]]]; recalling that multiplying polynomials together is the same as convolving the vectors of their corresponding coefficients. Jun 4 at 15:19

I might do it like this:

RollDND[throws_, sides_] :=
Total[ Table[RandomInteger[{1, sides}], throws]]

RollDND[4, 8]

• This is a good answer for simulating the roll, but I'm asking about modeling the distribution.
– nben
Jun 2 at 17:05

A single die roll:

DieRoll[sides_Integer] := RandomInteger[{1, sides}]


Multiple rolls of same type of die:

DiceRoll[diesides_Integer, number_Integer] :=
RandomInteger[{1, diesides}, number]


Multiple rolls with different dice (I don't know if that's actually a thing in D&D):

DiceRoll[diesides : {__Integer}] := DieRoll /@ diesides


If you just want the sums:

DiceRollSum[diesides : {__Integer}] := Total[DiceRoll[diesides]];
DiceRollSum[diesides_Integer, number_Integer] :=
Total[DiceRoll[diesides, number]]