First, I’ll import dplyr for analysis, charting and piping. I will also set a random seed so my results are consistent upon knitting.

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

set.seed(905)

Here, I’ll generate a for loop, setting i as an iterative number of trials. The sample function allows me to generate a random series of 3 values from 1 to 6, simulating my dice rolls. I set the argument replace equal to TRUE, since the rolls are independent of one another.

#initialize vector of dice roll sums
roll_sums <-c()

trials <- 50000
#run the simulation predetermined number of times

for (i in 1:trials) {
  dice_rolls <- sample(1:6, size=3, replace=TRUE)
  roll_sums <- c(roll_sums, sum(dice_rolls))
}

count_of_9s <- sum(roll_sums == 9)
count_of_10s <- sum(roll_sums == 10)

print(paste("Count of 9s:",count_of_9s, "a probability of ", count_of_9s/trials))

## [1] "Count of 9s: 5846 a probability of  0.11692"

print(paste("Count of 10s:",count_of_10s, "a probability of ", count_of_10s/trials))

## [1] "Count of 10s: 6100 a probability of  0.122"

Let’s visualize the full distribution:

roll_sums %>%
  table() %>%
  barplot()

As we can see, there were more dice rolls resulting in a sum of 10 than of 9 (and actually, 11 seems to have had the highest probability of all).

Now, I want to address the premise that the number of possible 3-die rolls adding up to 9 or 10 are the same. I can make a list of all possible rolls, then narrow down to those adding up to 9 or 10 respectively.

# Generate all possible outcomes for three rolls
possible_outcomes <- expand.grid(roll1 = 1:6, 
                                 roll2 = 1:6, 
                                 roll3 = 1:6)

# Calculate sum
possible_outcomes$sum = rowSums(possible_outcomes)

count_9sum <- possible_outcomes %>%
  filter(sum == 9) %>%
  nrow()

count_10sum <- possible_outcomes %>%
  filter(sum == 10) %>%
  nrow()

count_9sum / nrow(possible_outcomes)

## [1] 0.1157407

count_10sum / nrow(possible_outcomes)

## [1] 0.125

Indeed, of all possible 3-die rolls, more of them add up to 10 than 9. Perhaps 17th century gamblers aren’t the most trustworthy sources of information!