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?

  • 4
    $\begingroup$ 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. $\endgroup$
    – eyorble
    Jun 2, 2022 at 16:42
  • 2
    $\begingroup$ 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. $\endgroup$ Jun 2, 2022 at 17:49
  • 2
    $\begingroup$ 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. $\endgroup$ Jun 3, 2022 at 2:14

5 Answers 5


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}], 
 TableHeadings -> {None, {"Total", "P(Total)"}}]

$$ \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.

  • 3
    $\begingroup$ 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]. $\endgroup$ Jun 4, 2022 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_] := 
      Transpose[{#/Total[#], Range[Length[#]] - 1} &[
        CoefficientList[g, x]]], {0, _}]] /. {a_List, 
      b_List} :> (a -> b)];


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

{32.7313, Null}

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

  • $\begingroup$ Fixed d to dg as necessary, not sure I see anywhere that needs :> though. $\endgroup$
    – eyorble
    Jun 2, 2022 at 17:58
  • $\begingroup$ {a_List, b_List} -> (a -> b) $\endgroup$
    – Ben Izd
    Jun 2, 2022 at 18:01
  • $\begingroup$ 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? $\endgroup$
    – eyorble
    Jun 2, 2022 at 18:11
  • 3
    $\begingroup$ 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). $\endgroup$
    – Ben Izd
    Jun 2, 2022 at 18:18
  • 1
    $\begingroup$ 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. $\endgroup$ Jun 4, 2022 at 15:19

I might do it like this:

RollDND[throws_, sides_] := 
 Total[ Table[RandomInteger[{1, sides}], throws]]
RollDND[4, 8]
  • 4
    $\begingroup$ This is a good answer for simulating the roll, but I'm asking about modeling the distribution. $\endgroup$
    – nben
    Jun 2, 2022 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]]
  • 2
    $\begingroup$ This is a good answer about simulating the roll, but I'm asking about modeling the distribution analytically. $\endgroup$
    – nben
    Jun 2, 2022 at 17:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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