HP Forums

Pages: 1 2

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.

\(
\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}
\)

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

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.

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?

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%

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

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

\>>

(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?

Also the following links might be of interest:

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

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\)

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

(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.

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

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%

(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. Wink

(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 ===========================================================================​===================================================

(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

(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. Big Grin

Pages: 1 2

dm319

Thomas Klemm

dm319

Thomas Klemm

dm319

dm319

Albert Chan

Gil

Thomas Klemm

Thomas Klemm

Thomas Klemm

Thomas Klemm

Albert Chan

Allen

Dave Britten

Thomas Klemm

Dave Britten

Thomas Klemm

Thomas Klemm

Dave Britten