Dice Mechanics: Equal Probability Outcomes from D6

Multinomial Distribution for dice

We have two ways to tackle a distribution of dice.

1. Exhaustive Enumeration of the dice combinations.
2. Direct computation of the results.

Exhaustive enumeration is easy to implement. We can use the itertools.product function to repeat the values from \(\{1, 2, 3, 4, 5, 6\}\) the right number of times.

For 2D, this is \(\{1, 2, 3, 4, 5, 6\} \times \{1, 2, 3, 4, 5, 6\}\). It has \(6^2 = 36\) distinct pairs of values. It has 11 distinct sums between 2 and 12.

For 10D, there are \(6^{10} = 60,466,176\) distinct 10-tuples of values, with sums from 10 to 60.

Why would anyone roll 10D to make a fair choice among alternatives? While it doesn’t seem terribly practical, it seems prudent to at least explore the possibility.

The problem is that enumerating over 60 million dice alternatives, computing over 60 million sums, and then accumulating the values into a Counter object takes a lot of time. It only gets worse for more dice.

My box of 6-sided die has a dozen dice. \(6^{12} = 2,176,782,336\). Ouch, that’s a lot of combinations.

from collections import Counter
from collections.abc import Iterator
import itertools

def all_rolls(n: int, k: int = 6) -> Iterator[int]:
    """
    Create *all* alternatives that will define a distribution with n samples of k categories.
    Sum the n-tuples to create totals that range in value between n and n*k.
    """
    yield from (
        sum(x) for x in itertools.product(range(1, k + 1), repeat=n)
    )

def distribution_hard(n: int, k: int=6) -> Counter[int]:
    """Create a counter to summarize """
    return Counter(all_rolls(n, k))

sorted(distribution_hard(2).items())

[(2, 1),
 (3, 2),
 (4, 3),
 (5, 4),
 (6, 5),
 (7, 6),
 (8, 5),
 (9, 4),
 (10, 3),
 (11, 2),
 (12, 1)]

For direct computation of the distribution of dice, we need a lesson in multinomials.

See https://towardsdatascience.com/modelling-the-probability-distributions-of-dice-b6ecf87b24ea/. From this we learn the following (amongst other things). For \(n\) dice of \(s\) sides (6 in our case), the probability of getting the target value \(T\) is this:

equation GIF

(Above is the GIF from the original article.)

This is the formula, rewritten here in LaTeX:

\(P(n, s, T) = \Bigl(\sum\limits_{k=0}^{\lfloor \frac{T-n}{s} \rfloor}\bigl(-1\bigr)^k \frac{n!}{(n-k)!k!} \frac{(T-sk-1)!}{(T-sk-n)!(n-1)!}\Bigr)\Bigl(\frac{1}{s}\Bigr)^n\)

Yes, they’re the same.

We’ll use the following, with two optimizations. First, Python has the math.comb() function to compute binomials directly. We can leverage \(\frac{n!}{(n-k)!k!} = \binom{n}{k}\). Second, we’ll drop the \(\bigl(\frac{1}{s}\bigr)^n\) term because we don’t want a probability, we want the actual number of combinations of dice with the given total, \(T\).

\(P(n, s, T) = \sum\limits_{k=0}^{\lfloor \frac{T-n}{s} \rfloor}\bigl(-1\bigr)^k \binom{n}{k} \binom{T-sk-1}{n-1}\)

We can compute the distribution using \(P(n, s, T)\) instead of enumerating the product of \(n\) dice.

Also, as a Python implementation choice, we prefer to have the values of n and s last, since these don’t change as frequently as T. For our purposes, s has a default value of 6, which is not going to change for this analysis.

from math import comb, factorial
from functools import cache

@cache
def dice_prob(T, n, s=6):
    upper = (T - n) // s
    series = (
        (-1 if k % 2 == 1 else 1) * comb(n, k) * comb(T - s * k - 1, n - 1)
        for k in range(0, upper + 1)
    )
    return int(sum(series))

def distribution_2(n, s=6):
    return dict((t, dice_prob(t, n)) for t in range(n, n*6+1))

sorted(distribution_2(2).items())

[(2, 1),
 (3, 2),
 (4, 3),
 (5, 4),
 (6, 5),
 (7, 6),
 (8, 5),
 (9, 4),
 (10, 3),
 (11, 2),
 (12, 1)]

def test_distribution():
    """Takes over 3 minutes 40 seconds to demonstrate the two algorithms produce the same result."""
    for dice in range(2, 13):
        closed_form = distribution_2(dice, 6)
        hard_way = distribution_hard(dice, 6)
        assert closed_form == hard_way, f"mismatch for {dice=} {closed_form=} {hard_way=}"
    print("Identical!")

# test_distribution()

Mechanic Details

This data structure records the details of a dice-lookup table mechanic.

The Mechanic object will contain one or more DieRange instances.

A DieRange object lists the dice values from the source distribution used to create (approximately) equal-sized outcome bins. It also has the total number of dice combinations for the collection of values. When summarizing, the low and high attribute are important to summarize the range of values that define the outcome bin.

The Mechanic object is created with a list of DieRange instances. These define all the individual outcomes. The match with ideal bin-sizes isn’t going to be absolutely correct. The score and max_error are two ways to evaluate how will the mechanic fits the ideal.

from dataclasses import dataclass, field
from statistics import mean, variance

@dataclass
class DieRange:
    map_to: int
    values: list[int] = field(default_factory=list)
    total: int = field(default=0)

    @property
    def low(self):
        return min(self.values)

    @property
    def high(self):
        return max(self.values)

@dataclass
class Mechanic:
    result_size: int
    n_dice: int
    total_size: int
    mapping: list[DieRange] = field(default_factory=list)

    _actuals: list[float] = field(default_factory=list, init=False)

    @property
    def actuals(self):
        if not self._actuals:
            self._actuals = [dr.total / self.total_size for dr in self.mapping]
        return self._actuals

    @property
    def total(self):
        return sum(self.mapping)

    @property
    def score(self):
        expected = 1 / self.result_size
        return sum(
            (actual - expected)**2
            for actual in self.actuals
        )

    @property
    def max_error(self):
        expected = 1 / self.result_size
        errors = (
            actual - expected
            for actual in self.actuals
        )
        return max(errors, key=lambda v: abs(v))

The goal is to produce Mechanic instances for a given number of equal-probability outcomes, and a selected number of dice.

We can vary the number of dice for a given number of outcomes, looking for one that has a minimal error.

from collections import defaultdict

def make_mechanic(result_size: int, number_dice: int) -> Mechanic:
    """
    Computes the multinomial distribution for nD6.
    Formally, n samples of k=6 categories.
    Then, breaks the distribution into `result_size` distinct buckets.
    Return a summary of the Mechanic to make it easier to compare alternatives.

    :param result_size: number of outcomes
    :param number_dice: number of dice to try and map to outcomes
    :returns: a Mechanic with details of a table to provide (approximately) equal probability outcomes
    """
    distribution = distribution_2(number_dice)

    n_tile_width = sum(distribution.values()) / result_size
    # print(f"{n_tile_width:.1f} out of {sum(distribution.values())}")
    rolls = defaultdict(list)
    tile_num = 0
    accum_fq = 0
    for k in sorted(distribution):
        accum_fq = accum_fq + distribution[k]
        if accum_fq >= n_tile_width:
            tile_num += 0 if tile_num == result_size-1 else 1
            accum_fq = accum_fq - n_tile_width
        rolls[tile_num].append(k)
    m = Mechanic(
        result_size, number_dice, sum(distribution.values()),
        [
            DieRange(r, rolls[r], sum(distribution[n] for n in rolls[r]))
            for r in rolls
        ]
    )
    return m

The report() function explores the alternative space. It creates a bunch of alternative mechanics, and then ranks them by the score computation.

In some cases, the best fit also involves a big handful of dice. A slightly worse fit may involve substantially fewer dice. It’s not a simple matter of picking the “best” fit; it also has to be reasonably useful by actual players.

def report(outcomes: int = 5, max_dice: int | None = None):
    """
    Locate optimal Mechanics for a given number of outcomes.
    This prints a summary of each mechanic considered.

    :param outcomes: Number of distinct, equal-probability outcomes
    :param max_dice: Maximum number of dice to consider, the default is the number of outcomes.
    """
    mechanics = [
        make_mechanic(outcomes, n_dice)
        for n_dice in range(2, max_dice or outcomes+1)
    ]
    for m in sorted(mechanics, key=lambda m: m.score):
        print(f"{m.result_size} choices with {m.n_dice}D")
        print(f"| roll | result | frequency | error |")
        print(f"| ---- | ------ | --------- | ----- |")
        expected = 1 / outcomes
        for dr in m.mapping:
            actual = dr.total / m.total_size
            print(f"| {dr.low}-{dr.high} | {dr.map_to} | {actual:.1%} | {actual - expected:.1%} |")
        print(f"{m.score=:.3f} {m.max_error=:.1%}")
        print()

Let’s use this report() function to compute a number of different kinds of equal-probability mechanics.

report(5)

5 choices with 4D
| roll | result | frequency | error |
| ---- | ------ | --------- | ----- |
| 4-10 | 0 | 15.9% | -4.1% |
| 11-12 | 1 | 17.7% | -2.3% |
| 13-14 | 2 | 22.1% | 2.1% |
| 15-16 | 3 | 20.4% | 0.4% |
| 17-24 | 4 | 23.9% | 3.9% |
m.score=0.004 m.max_error=-4.1%

5 choices with 3D
| roll | result | frequency | error |
| ---- | ------ | --------- | ----- |
| 3-7 | 0 | 16.2% | -3.8% |
| 8-9 | 1 | 21.3% | 1.3% |
| 10-10 | 2 | 12.5% | -7.5% |
| 11-12 | 3 | 24.1% | 4.1% |
| 13-18 | 4 | 25.9% | 5.9% |
m.score=0.012 m.max_error=-7.5%

5 choices with 5D
| roll | result | frequency | error |
| ---- | ------ | --------- | ----- |
| 5-13 | 0 | 15.2% | -4.8% |
| 14-16 | 1 | 24.8% | 4.8% |
| 17-17 | 2 | 10.0% | -10.0% |
| 18-20 | 3 | 27.9% | 7.9% |
| 21-30 | 4 | 22.1% | 2.1% |
m.score=0.021 m.max_error=-10.0%

5 choices with 2D
| roll | result | frequency | error |
| ---- | ------ | --------- | ----- |
| 2-4 | 0 | 16.7% | -3.3% |
| 5-5 | 1 | 11.1% | -8.9% |
| 6-7 | 2 | 30.6% | 10.6% |
| 8-8 | 3 | 13.9% | -6.1% |
| 9-12 | 4 | 27.8% | 7.8% |
m.score=0.030 m.max_error=10.6%

We'll omit all the other experiments we ran to fill in the TL;DR table up front.