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

Question

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,
  Thread@Distributed[
     {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?

Answer 1

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.

Answer 2

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.

Answer 3

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

Answer 4

I might do it like this:

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

RollDND[4, 8]

Answer 5

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