+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Calculators (and very old HP Computers) (/forum-3.html)
+--- Forum: General Forum (/forum-4.html)
+--- Thread: Dice probabilities (/thread-21671.html)
Pages: 1 2
Dice probabilities - dm319 - 04-29-2024 04:23 PM
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.
RE: Dice probabilities - Thomas Klemm - 04-29-2024 04:55 PM
\(
\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}
\)
RE: Dice probabilities - dm319 - 04-29-2024 05:49 PM
That looks good to me, but can you explain for someone who is not so bright?
RE: Dice probabilities - Thomas Klemm - 04-29-2024 06:07 PM
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.
RE: Dice probabilities - dm319 - 04-29-2024 07:14 PM
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?
RE: Dice probabilities - dm319 - 04-29-2024 07:22 PM
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%
RE: Dice probabilities - Albert Chan - 04-29-2024 09:22 PM
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
RE: Dice probabilities - Gil - 04-29-2024 10:36 PM
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
\>>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
\>>RE: Dice probabilities - Thomas Klemm - 04-30-2024 12:27 AM
(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 ENDInitialisation
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?
RE: Dice probabilities - Thomas Klemm - 04-30-2024 12:40 AM
Also the following links might be of interest:
RE: Dice probabilities - Thomas Klemm - 04-30-2024 06:43 AM
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 ENDIt 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
RE: Dice probabilities - Thomas Klemm - 04-30-2024 08:06 AM
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\)
RE: Dice probabilities - Albert Chan - 04-30-2024 10:59 AM
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 51P(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+T5+(Te5-5))/3^5 = (1+T5+Te5)/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
RE: Dice probabilities - Allen - 04-30-2024 11:19 AM
(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.
RE: Dice probabilities - Dave Britten - 04-30-2024 12:00 PM
Here's a short Python program I wrote a long time ago for generating probability tables:
https://gist.github.com/dave-britten/c4b54c879d65ce54497b
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))RE: Dice probabilities - Thomas Klemm - 04-30-2024 05:11 PM
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%RE: Dice probabilities - Dave Britten - 04-30-2024 05:58 PM
(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.
RE: Dice probabilities - Thomas Klemm - 04-30-2024 06:04 PM
(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) == expectedTo 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 ==============================================================================================================================RE: Dice probabilities - Thomas Klemm - 04-30-2024 06:25 PM
(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
RE: Dice probabilities - Dave Britten - 04-30-2024 07:24 PM
(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.
