Dice probabilities

[dm319]

So I came across this question while playing a new board game at the weekend.

This board game has multiple dice which are unusual - they have 2 blank faces, 2 'one dot' faces and 2 'two dot' faces. I.e. a 1/3rd chance to roll a zero, one or two.

In this game we roll multiple dice and then sum the number of dots to score that roll.

What is the probability of rolling a score of more than or equal to 7 when I have 5 of these dice in my hand?

Does anyone know of a good way to calculate this?

I've thought of two approaches, both of which are a bit inelegant. I'm curious to see if there is a neater way especially for the general case, and whether a program could be written to calculate it.

[Thomas Klemm]

\(
\begin{align}
(x^2+x+1)^5 &= x^{10} + 5 x^9 + 15 x^8 + 30 x^7 + 45 x^6 + 51 x^5 + 45 x^4 + 30 x^3 + 15 x^2 + 5 x + 1 \\
\\
p &= \frac{1 + 5 + 15 + 30}{3^5} \\
&= \frac{51}{243} \\
&= \frac{17}{81} \\
&\approx 0.20988
\end{align}
\)

[dm319]

That looks good to me, but can you explain for someone who is not so bright?

[Thomas Klemm]

We might rather write it: \(x^2+x^1+x^0\)
Thus for an ordinary die we would use: \(x^6+x^5+x^4+x^3+x^2+x^1\)
If these terms get multiplied, their exponents are added.
We're using the distributive law.
Every possible case gets combined with all the others.
And the exponents are added up for us.

I assume that the reason we add up the coefficients of \(x^{10} \cdots x^7\) is clear.

[dm319]

Oh that is fascinating.

So you are using the powers of x and the distributive law to 'add up' the powers of x and therefore essentially generate the crossing of 0, 1, 2.

Are you just using the convenience of this in order to calculate that - because I don't think there is anything inherently 'squared' about a die with two dots, or I may be missing something.

Also I imagine it took a while to expand out?

[dm319]

This was my simple method simply generating every possible roll using the function `expand.grid` in R, then summing, and finding the number with scores >= 7 (and other scores as well as other number of dice).

Code:


library(tidyverse)

faces <- c(0, 1, 2)
max_dice_needed <- 8
max_score_needed <- 12

results <- 
data.frame(dice = numeric(),
           score_needed = numeric(),
           prob = numeric()) 

for(k in 1:max_dice_needed){

  ndice <- k
  dice <- list()

  # add number of dice to list
  for (i in 1:ndice){
    dice <- append(dice, list(faces))
  }

  # find every permutation of dice faces
  rolls <- expand.grid(dice)

  # calculate total score for the throw
  scores <-
  rolls %>%
    rowwise %>%
    mutate(score = sum(c_across(starts_with("Var")))) %>%
    pull(score)

  message("dice:", ndice, " scores:", length(scores))

  # calculate probability of getting >= required score
  for(j in 0:max_score_needed){

    score_needed <- j


    results <- 
    rbind(results, 
    data.frame(dice = ndice,
               score_needed = j,
               prob = mean(scores >= score_needed)
               )
    )
  }
}

# plot results
results %>%
  mutate(dice = factor(dice)) %>%
  ggplot(aes(x = score_needed, y = prob, colour = dice))+
  geom_point()+
  geom_line()+
  geom_label(aes(label = dice))+
  scale_x_continuous(breaks = scales::pretty_breaks(n = 12))+
  scale_y_continuous(breaks = scales::pretty_breaks(n = 10),
                     labels = scales::percent)+
  labs(y = "probability", x = "score needed")+
  theme(legend.position = "none")

ggsave("probabilities.png")

Another method I had was to consider the number of combinations of dice that made up 7 or higher, i.e.:

7: 2, 2, 2, 1, 0 - 5!/3!
7: 2, 2, 1, 1, 1 - 5!/(2! x 3!)
8: 2, 2, 2, 2, 0 - 5!/4!
8: 2, 2, 2, 1, 1 - 5!(2! x 3!)
9: 2, 2, 2, 2, 1 - 5!/4!
10: 2, 2, 2, 2, 2 - 1

then (5!/3! + 5!/(2! x 3!) + 5!/4! + 5!(2! x 3!) + 5!/4! + 5!/4! + 1) / 3^5 = 20.99%

[Albert Chan]

Dice = (x^2 + x^1 + x^0) / 3 = 100%      // see Probability-generating function
Dice^5 = 100% probability too.

From this distribution, just pick what we needed.

>10101^5                     // (x^2+x+1)^5, for x=100
 1.05153045515E20
>(1+5+15+30+45+51+45+30+15+5+1)/3^5
 1
>(1+5+15+30)/3^5      // P(sum >= 7)
 .20987654321

[Gil]

For a program with HP50G

Just simulate, one time, the case of having such 5 dice:

Code:

\<< 0. 0. 1. 5.
  START RAND \-> n
    \<<
      CASE n 3. INV \<=
        THEN 0.
        END n 3. INV 2. * \<=
        THEN 1.
        END 2.
      END
    \>> +
  NEXT 7. \>=
  IF
  THEN 1. +
  END
\>>

and call that subroutine P.

Just simulate now the above subroutine P, not once, but 10000 times as follows:
\<< 0. 1. 10000.
START P +
NEXT 10000. /

and you should find about 0.21.
\>>

To speed up the subroutine P, you can change the code into:

Code:

\<< 0. 0. 1. 5.
  START RAND DUP
    CASE 3. INV \<=
      THEN DROP 0.
      END 3. INV 2. * \<=
      THEN 1.
      END 2.
    END +
  NEXT 7. \>=
  IF
  THEN 1. +
  END
\>>

[Thomas Klemm]

(04-29-2024 07:14 PM)dm319 Wrote:  Also I imagine it took a while to expand out?

Not really: (x^2+x+1)^5

---

(04-29-2024 04:23 PM)dm319 Wrote:  (…) whether a program could be written to calculate it.

Here's a program for the HP-42S:

Code:

00 { 43-Byte Prgm }
01 1
02 3
03 NEWMAT
04 1
05 +
06 ×
07 STO "A"
08 DIM?
09 +
10 LASTX
11 DIM "A"
12 RCL "A"
13 TRANS
14 STO "A"
15 DIM?
16 1
17 -
18 DIM "A"
19 RCL "A"
20 TRANS
21 RSUM
22 END

Initialisation

3 ENTER 1
NEWMAT
1 +

x: [ 3×1 Matrix ]

Loop

R/S

x: [ 5×1 Matrix ]

R/S

x: [ 7×1 Matrix ]

R/S

x: [ 9×1 Matrix ]

R/S

x: [ 11×1 Matrix ]

The result is:

\(
\begin{bmatrix}
1 \\
5 \\
15 \\
30 \\
45 \\
51 \\
45 \\
30 \\
15 \\
5 \\
1 \\
\end{bmatrix}
\)

---

We can now also ask:

• What about dice with 1 blank faces, 2 'one dot' faces and 3 'two dot' faces?
• What about dice with 3 blank faces, 2 'one dot' faces and 1 'two dot' faces?
  

Can you come up with similar formulas for their corresponding probabilities?

[Thomas Klemm]

Also the following links might be of interest:

• A027907
• Trinomial Triangle
• Trinomial Coefficient
  

[Thomas Klemm]

Here's a slightly improved program for the HP-42S:

Code:

00 { 52-Byte Prgm }
01 STO 00
02 R↓
03 STO "A"
04 GTO 01
05▸LBL 00
06 RCL "A"
07 TRANS
08 ×
09 STO "_"
10 DIM?
11 +
12 LASTX
13 DIM "_"
14 RCL "_"
15 TRANS
16 STO "_"
17 DIM?
18 1
19 -
20 DIM "_"
21 RCL "_"
22 TRANS
23 RSUM
24▸LBL 01
25 DSE 00
26 GTO 00
27 END

It allows to use an arbitrary [ k×1 Matrix ] and an integer exponent \(n > 0\).

Example

Use matrix [ 1 2 3 ] and exponent \(n = 5\).

3 ENTER 1 NEWMAT
EDIT
1 →
2 →
3
EXIT

5

R/S

x: [ 11×1 Matrix ]

The result is:

\(
\begin{bmatrix}
1 \\
10 \\
55 \\
200 \\
530 \\
1052 \\
1590 \\
1800 \\
1485 \\
810 \\
243 \\
\end{bmatrix}
\)

Compare this to: (x^2+2x+3)^5

[Thomas Klemm]

To illustrate the program, here's an example of a single execution of the loop:

Code:

RCL "A"
TRANS
×

\(
\begin{bmatrix}
1 \\
4 \\
10 \\
12 \\
9 \\
\end{bmatrix}
\begin{bmatrix}
1 & 2 & 3
\end{bmatrix}
=
\begin{bmatrix}
1 & 2 & 3 \\
4 & 8 & 12 \\
10 & 20 & 30 \\
12 & 24 & 36 \\
9 & 18 & 27 \\
\end{bmatrix}
\)

Code:

STO "_"
DIM?
+
LASTX
DIM "_"

\(
\begin{bmatrix}
1 & 2 & 3 \\
4 & 8 & 12 \\
10 & 20 & 30 \\
12 & 24 & 36 \\
9 & 18 & 27 \\
0 & 0 & 0 \\
0 & 0 & 0 \\
0 & 0 & 0 \\
\end{bmatrix}
\)

Code:

RCL "_"
TRANS
STO "_"

\(
\begin{bmatrix}
1 & 4 & 10 & 12 & 9 & 0 & 0 & 0 \\
2 & 8 & 20 & 24 & 18 & 0 & 0 & 0 \\
3 & 12 & 30 & 36 & 27 & 0 & 0 & 0 \\
\end{bmatrix}
\)

Code:

DIM?
1
-
DIM "_"

\(
\begin{bmatrix}
1 & 4 & 10 & 12 & 9 & 0 & 0 \\
0 & 2 & 8 & 20 & 24 & 18 & 0 \\
0 & 0 & 3 & 12 & 30 & 36 & 27 \\
\end{bmatrix}
\)

Code:

RCL "_"
TRANS

\(
\begin{bmatrix}
1 & 0 & 0 \\
4 & 2 & 0 \\
10 & 8 & 3 \\
12 & 20 & 12 \\
9 & 24 & 30 \\
0 & 18 & 36 \\
0 & 0 & 27 \\
\end{bmatrix}
\)

Code:

RSUM

\(
\begin{bmatrix}
1 \\
6 \\
21 \\
44 \\
63 \\
54 \\
27 \\
\end{bmatrix}
\)

Compare this to: \((x^4+4x^3+10x^2+12x+9)(x^2+2x+3)\) \(= x^6 + 6 x^5 + 21 x^4 + 44 x^3 + 63 x^2 + 54 x + 27\)

[Albert Chan]

Instead of multiply to expand polynomial, we could also do division.

(x^2+x+1)^5 = (x^3-1)^5 / (x-1)^5

Numerator coefficients is simply 5th row of pascal triangle, with alternate sign:

x^15 - 5x^12 + 10x^9 - 10x^6 + 5x^3 - 1

Synthetic Division, (x-1) 5 times, we get the expanded polynomial needed.
With symmetry, we only need 2 terms. (for exponent >= 7, we don't even need last 2 columns)

Code:

1> 1   0   0   -5    0    0 
1> 1   1   1   -4   -4   -4
1> 1   2   3   -1   -5   -9
1> 1   3   6    5    0   -9
1> 1   4   10  15   15    6 
   1   5   15  30   45   51

P(sum >= 7) = (1+5+15+30)/3^5 = 17/81 ≈ 0.2099

Or, we skip synthetic divisions, and go straight for answer

P(sum >= 7) = (1+5+T₅+(Te₅-5))/3^5 = (1+T₅+Te₅)/3^5 = (1+comb(6,2)+comb(7,3))/3^5 = 17/81

https://en.wikipedia.org/wiki/Triangular_number
https://en.wikipedia.org/wiki/Tetrahedral_number

[Allen]

(04-29-2024 04:55 PM)Thomas Klemm Wrote:  \(
\begin{align}
(x^2+x+1)^5 &= x^{10} + 5 x^9 + 15 x^8 + 30 x^7 + 45 x^6 + 51 x^5 + 45 x^4 + 30 x^3 + 15 x^2 + 5 x + 1 \\
\\
p &= \frac{1 + 5 + 15 + 30}{3^5} \\
&= \frac{51}{243} \\
&= \frac{17}{81} \\
&\approx 0.20988
\end{align}
\)

@dm If you like this kind of wizardy with counting, which Mr. Klemm and others here are simply world-class at doing on any platform, may I suggest A book called "Generatingfunctionology" by Wilf ( HTML page with terms)?

Related (unevaluated) Youtube lectures.

[Dave Britten]

Here's a short Python program I wrote a long time ago for generating probability tables:

https://gist.github.com/dave-britten/c4b...65ce54497b

Roughly half of it is just processing the user input! Enter one die at a time, and enter a blank line to terminate input.

EDIT: It's so old, in fact, that I just realized it's for Python 2! Here's a quick update to Python 3:

Code:

import re

dice = [] #Input dice as polynomial coefficients    

print("Enter face values as non-negative integers, separated by spaces, commas, or semicolons. (e.g. \"1,2;3 4 5 6\")")
while True:
    rawdie = [] #List of face values
    qty = 1 #Quantity of a single die spec

    inp = input("Die " + str(len(dice) + 1) + "? ")
    if inp == "":
        break
    #rawdie = [int(face) for face in inp.split(" ")]
    rawdie = [int(face.strip()) for face in re.split(" |,|;", inp) if face.strip()]
    qty = input("Quantity? [1] ")

    if qty == "":
        qty = 1
    else:
        qty = int(qty)

    die = [0] * (max(rawdie) + 1)
    for i in rawdie:
        die[i] += 1

    dice += [die for i in range(qty)]

    print(str(rawdie) + " x " + str(qty))

lastresult = dice[0]

for current in dice[1:]:
    result = [0] * (len(lastresult) + len(current) - 1)
    for x, v1 in enumerate(lastresult):
        if v1 == 0:
            continue
        for y, v2 in enumerate(current):
            result[x + y] += v1 * v2
    lastresult = result

skipzero = True
p = 0.0
fcum = 0.0
for i in lastresult:
    p += i

print("Permutations: {:.0f}".format(p))
print("{:<6} {:<12} {:>10} {:>10} {:>10}".format("Total","Freq","Prob", "p>=", "p<="))
for i, f in enumerate(lastresult):
    fcum += f
    if f == 0 and skipzero:
        continue
    skipzero = False
    print("{:<6} {:<12} {:>10.3%} {:>10.3%} {:>10.3%}".format(i, f, f / p, (p - fcum + f) / p, fcum / p))

[Thomas Klemm]

I took the liberty to refactor your Python program:

Code:

import re


def input_dice_as_polynomial_coefficients():
    dice = []

    print(
        'Enter face values as non-negative integers, separated by spaces, commas, or semicolons. (e.g. "1,2;3 4 5 6")'
    )
    while inp := input("Die " + str(len(dice) + 1) + "? "):
        # List of face values
        raw_die = [
            int(face.strip()) for face in re.split(r"[ ,;]", inp) if face.strip()
        ]
        # Quantity of a single die spec
        raw_qty = input("Quantity? [1] ")
        qty = 1 if raw_qty == "" else int(raw_qty)
        die = coefficients(raw_die)
        dice += [die] * qty
        print(f"{raw_die} x {qty}")
    return dice


def coefficients(raw_die):
    die = [0] * (max(raw_die) + 1)
    for i in raw_die:
        die[i] += 1
    return die


def calculate_frequencies(dice):
    previous, *rest = dice

    for current in rest:
        result = [0] * (len(previous) + len(current) - 1)
        for i, u in enumerate(previous):
            if u == 0:
                continue
            for j, v in enumerate(current):
                result[i + j] += u * v
        previous = result
    return previous


def print_permutations(frequencies):
    skip_zero = True
    p = sum(frequencies)
    f_cum = 0.0

    print(f"Permutations: {p:.0f}")
    print(f"{'Total':<6} {'Freq':<12} {'Prob':>10} {'p>=':>10} {'p<=':>10}")
    for i, f in enumerate(frequencies):
        f_cum += f
        if f == 0 and skip_zero:
            continue
        skip_zero = False
        print(
            f"{i:<6} {f:<12} {f / p:>10.3%} {(p - f_cum + f) / p:>10.3%} {f_cum / p:>10.3%}"
        )


if __name__ == "__main__":
    dice = input_dice_as_polynomial_coefficients()
    frequencies = calculate_frequencies(dice)
    print_permutations(frequencies)

• extract method to break the program into smaller pieces
• this allows to write unit-tests for some of the functions
• use the __name__ == "__main__" check
• this allows to import the file as module without running it
• use f-strings
• use the walrus-operator :=
• rename some variables
• use black to format the code
  

Examples

Code:

Enter face values as non-negative integers, separated by spaces, commas, or semicolons. (e.g. "1,2;3 4 5 6")
Die 1? 0 1 2
Quantity? [1] 5
[0, 1, 2] x 5
Die 6? 
Permutations: 243
Total  Freq               Prob        p>=        p<=
0      1                0.412%   100.000%     0.412%
1      5                2.058%    99.588%     2.469%
2      15               6.173%    97.531%     8.642%
3      30              12.346%    91.358%    20.988%
4      45              18.519%    79.012%    39.506%
5      51              20.988%    60.494%    60.494%
6      45              18.519%    39.506%    79.012%
7      30              12.346%    20.988%    91.358%
8      15               6.173%     8.642%    97.531%
9      5                2.058%     2.469%    99.588%
10     1                0.412%     0.412%   100.000%

Code:

Enter face values as non-negative integers, separated by spaces, commas, or semicolons. (e.g. "1,2;3 4 5 6")
Die 1? 0 0 0 1 1 2
Quantity? [1] 5
[0, 0, 0, 1, 1, 2] x 5
Die 6? 
Permutations: 7776
Total  Freq               Prob        p>=        p<=
0      243              3.125%   100.000%     3.125%
1      810             10.417%    96.875%    13.542%
2      1485            19.097%    86.458%    32.639%
3      1800            23.148%    67.361%    55.787%
4      1590            20.448%    44.213%    76.235%
5      1052            13.529%    23.765%    89.763%
6      530              6.816%    10.237%    96.579%
7      200              2.572%     3.421%    99.151%
8      55               0.707%     0.849%    99.859%
9      10               0.129%     0.141%    99.987%
10     1                0.013%     0.013%   100.000%

[Dave Britten]

(04-30-2024 05:11 PM)Thomas Klemm Wrote:  I took the liberty to refactor your Python program:

Code:

import re


def input_dice_as_polynomial_coefficients():
    dice = []

    print(
        'Enter face values as non-negative integers, separated by spaces, commas, or semicolons. (e.g. "1,2;3 4 5 6")'
    )
    while inp := input("Die " + str(len(dice) + 1) + "? "):
        # List of face values
        raw_die = [
            int(face.strip()) for face in re.split(r"[ ,;]", inp) if face.strip()
        ]
        # Quantity of a single die spec
        raw_qty = input("Quantity? [1] ")
        qty = 1 if raw_qty == "" else int(raw_qty)
        die = coefficients(raw_die)
        dice += [die] * qty
        print(f"{raw_die} x {qty}")
    return dice


def coefficients(raw_die):
    die = [0] * (max(raw_die) + 1)
    for i in raw_die:
        die[i] += 1
    return die


def calculate_frequencies(dice):
    previous, *rest = dice

    for current in rest:
        result = [0] * (len(previous) + len(current) - 1)
        for i, u in enumerate(previous):
            if u == 0:
                continue
            for j, v in enumerate(current):
                result[i + j] += u * v
        previous = result
    return previous


def print_permutations(frequencies):
    skip_zero = True
    p = sum(frequencies)
    f_cum = 0.0

    print(f"Permutations: {p:.0f}")
    print(f"{'Total':<6} {'Freq':<12} {'Prob':>10} {'p>=':>10} {'p<=':>10}")
    for i, f in enumerate(frequencies):
        f_cum += f
        if f == 0 and skip_zero:
            continue
        skip_zero = False
        print(
            f"{i:<6} {f:<12} {f / p:>10.3%} {(p - f_cum + f) / p:>10.3%} {f_cum / p:>10.3%}"
        )


if __name__ == "__main__":
    dice = input_dice_as_polynomial_coefficients()
    frequencies = calculate_frequencies(dice)
    print_permutations(frequencies)

Cool! There's definitely plenty of room to fancy it up, especially with the 9 or so years of Python feature additions since I wrote it.

I also did a Turbo Pascal 5.5 port if anybody needs to do this on a 200LX.

[Thomas Klemm]

(04-30-2024 05:11 PM)Thomas Klemm Wrote:  

• this allows to write unit-tests for some of the functions
  

Here's an example of parametrized tests using pytest:

Code:

import pytest
from dice import calculate_frequencies, coefficients


@pytest.mark.parametrize(
    "raw_die, expected",
    [
        ([0, 1, 2], [1, 1, 1]),
        ([0, 0, 0, 1, 1, 2], [3, 2, 1]),
    ],
)
def test_coefficients(raw_die, expected):
    assert coefficients(raw_die) == expected


@pytest.mark.parametrize(
    "dice, expected",
    [
        ([[1, 1, 1]] * 5, [1, 5, 15, 30, 45, 51, 45, 30, 15, 5, 1]),
        ([[3, 2, 1]] * 5, [243, 810, 1485, 1800, 1590, 1052, 530, 200, 55, 10, 1]),
    ],
)
def test_calculate_frequencies(dice, expected):
    assert calculate_frequencies(dice) == expected

To run the tests, pip install pytest into the virtual environment and execute:

pytest test_dice.py

Code:

================================================================================​============================================ test session starts ===========================================================================​==================================================
platform darwin -- Python 3.12.2, pytest-8.2.0, pluggy-1.5.0
rootdir: /Users/whatever
collected 4 items

test_dice.py ....                                                                                                                                                                                                                                                      [100%]

================================================================================​============================================= 4 passed in 0.01s ===========================================================================​===================================================

[Thomas Klemm]

(04-30-2024 12:00 PM)Dave Britten Wrote:  Here's a short Python program I wrote a long time ago for generating probability tables:

Look what I found: Calculating odds of rolling a given total with a specified number of dice
Cf. Post #14

[Dave Britten]

(04-30-2024 06:25 PM)Thomas Klemm Wrote:  

(04-30-2024 12:00 PM)Dave Britten Wrote:  Here's a short Python program I wrote a long time ago for generating probability tables:

Look what I found: Calculating odds of rolling a given total with a specified number of dice
Cf. Post #14

Ah ha! There's that thread I wasn't coming up with the right combination of search terms to find.