In the early 1600s, Galileo was asked to explain the fact that, although the number of triples of integers from 1 to 6 with sum 9 is the same as the number of such triples with sum 10, when three dice are rolled, a 9 seemed to come up less often than a 10—supposedly in the experience of gamblers.
Write a program to simulate the roll of three dice a large number of times and keep track of the proportion of times that the sum is 9 and the proportion of times it is 10.
Can you conclude from your simulations that the gamblers were correct?
I decided to figure this out using probability logic as well as a simulation.
First, I was a bit confused by the definition of the “triples of integers from 1 to 6”. But I eventually decided that they meant combination with replacement.
Let’s make a data frame in R representing all combinations with replacement for the integers from 1 to 6.
combinations <- data.frame(First = rep(c(1,2,3,4,5,6),times=c(21,15,10,6,3,1)),
Second = c(rep(c(1,2,3,4,5,6),times=c(6,5,4,3,2,1)),rep(c(2,3,4,5,6),times=c(5,4,3,2,1)),rep(c(3,4,5,6),times=c(4,3,2,1)),rep(c(4,5,6),times=c(3,2,1)),rep(c(5,6),times=c(2,1)),6),
Third = rep(NA,times=56),
stringsAsFactors=FALSE)
values_to_fill_in_for_third <- c(1:6,2:6,3:6,4:6,5:6,6)
for(i in 1:6)
{
if(i > 1){values_to_fill_in_for_third <- values_to_fill_in_for_third[(9 - i):length(values_to_fill_in_for_third)]}
combinations$Third[combinations$First == i] <- values_to_fill_in_for_third
}
combinations## First Second Third
## 1 1 1 1
## 2 1 1 2
## 3 1 1 3
## 4 1 1 4
## 5 1 1 5
## 6 1 1 6
## 7 1 2 2
## 8 1 2 3
## 9 1 2 4
## 10 1 2 5
## 11 1 2 6
## 12 1 3 3
## 13 1 3 4
## 14 1 3 5
## 15 1 3 6
## 16 1 4 4
## 17 1 4 5
## 18 1 4 6
## 19 1 5 5
## 20 1 5 6
## 21 1 6 6
## 22 2 2 2
## 23 2 2 3
## 24 2 2 4
## 25 2 2 5
## 26 2 2 6
## 27 2 3 3
## 28 2 3 4
## 29 2 3 5
## 30 2 3 6
## 31 2 4 4
## 32 2 4 5
## 33 2 4 6
## 34 2 5 5
## 35 2 5 6
## 36 2 6 6
## 37 3 3 3
## 38 3 3 4
## 39 3 3 5
## 40 3 3 6
## 41 3 4 4
## 42 3 4 5
## 43 3 4 6
## 44 3 5 5
## 45 3 5 6
## 46 3 6 6
## 47 4 4 4
## 48 4 4 5
## 49 4 4 6
## 50 4 5 5
## 51 4 5 6
## 52 4 6 6
## 53 5 5 5
## 54 5 5 6
## 55 5 6 6
## 56 6 6 6Now, which of these combinations have a sum of 9 vs. a sum of 10?
print("Sum of 9:")## [1] "Sum of 9:"combinations[rowSums(combinations) == 9,]## First Second Third
## 11 1 2 6
## 14 1 3 5
## 16 1 4 4
## 25 2 2 5
## 28 2 3 4
## 37 3 3 3print("Sum of 10:")## [1] "Sum of 10:"combinations[rowSums(combinations) == 10,]## First Second Third
## 15 1 3 6
## 17 1 4 5
## 26 2 2 6
## 29 2 3 5
## 31 2 4 4
## 38 3 3 4Compare to the sums from all possible permutations of the integers 1 to 6 given 3 slots.
sum_per_permutation <- rep(NA,times=216)
permutations <- vector("list",length=216)
counter = 0
for(i in 1:6)
{
for(j in 1:6)
{
for(k in 1:6)
{
counter = counter + 1
sum_per_permutation[counter] <- i + j + k
permutations[[counter]] <- c(i,j,k)
}
}
}
table(sum_per_permutation)## sum_per_permutation
## 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
## 1 3 6 10 15 21 25 27 27 25 21 15 10 6 3 1We find that if we use permutations instead of combinations, a sum of 10 has a probability of 27/216, while a sum of 9 has a probability of 25/216.
Permutations are a more accurate mathematical representation of rolling dice than combination with replacement.
The major difference is the sixth combination for each sum.
For the combination 3,3,3, there is only one way to roll it.
Versus the combination 3,3,4, which could be rolled as:
All three of these are equally likely, increasing the chance of getting a sum of 10 relative to getting a sum of 9.
Let’s run a simulation of rolling 3 dice a bunch of times.
The maximum standard error we will get 95% of the time is roughly equal to 1/sqrt(n).
Here, this means we need to set n quite large.
For a standard error no larger than +/- 1/216, we would need 46,656 trials.
And having a standard error of a full 1/216 still seems quite large.
On modern computers, running this kind of simulation 1 million times is not that crazy.
So let’s do that. This will give a standard error of around 0.2 per 216, which seems a lot better.
set.seed(1392)
sum_of_dice <- rep(NA,times=1e6)
for(i in 1:1e6)
{
sum_of_dice[i] <- sum(sample(1:6,3,replace=TRUE))
}
print("Mathemetical probabilities as proportions:")## [1] "Mathemetical probabilities as proportions:"round(table(sum_per_permutation)/216,digits=4)## sum_per_permutation
## 3 4 5 6 7 8 9 10 11 12
## 0.0046 0.0139 0.0278 0.0463 0.0694 0.0972 0.1157 0.1250 0.1250 0.1157
## 13 14 15 16 17 18
## 0.0972 0.0694 0.0463 0.0278 0.0139 0.0046print("Probabilities per sum:")## [1] "Probabilities per sum:"round(table(sum_of_dice)/1e6,digits=4)## sum_of_dice
## 3 4 5 6 7 8 9 10 11 12
## 0.0047 0.0140 0.0278 0.0463 0.0695 0.0967 0.1159 0.1255 0.1245 0.1156
## 13 14 15 16 17 18
## 0.0973 0.0694 0.0465 0.0278 0.0139 0.0047The probabilities of getting a sum of 9 or 10 are 0.1159 and 0.1255 respectively according to our simulation (or 25.03/216 and 27.11/216 respectively).
This is very very close to the mathematically predicted probabilities of 0.1157 and 0.1250 (or 25/216 and 27/216 respectively).