Probability of getting exactly 2 6's out of 4 dice rolls

two_sixes.R

library(tidyverse)

ROLLS <- 4
SIMS <- 10000
  
#calculate
ncombos <- choose(ROLLS,2)
two_sixes <- 1/6 * 1/6 * 5/6 * 5/6

print(sprintf("choose(4, 2) * 25/(6^4) = %f", ncombos * two_sixes))

# Verify with Simulation

# simulate 10000 instances of 4 dice rolls
rolls <- sample(1:6, ROLLS * SIMS, replace=TRUE)
data <- data.frame(matrix(rolls, nrow = SIMS, ncol = ROLLS))
names(data) <- c("Roll1", "Roll2", "Roll3", "Roll4")
head(data)

# mutate it to be a 1 for each 6, and 0 for any other number
just_sixes <- sapply(rolls, FUN = function(x){ ifelse(x == 6, 1, 0)})
data <- data.frame(matrix(just_sixes, nrow = SIMS, ncol = ROLLS))
names(data) <- c("Roll1", "Roll2", "Roll3", "Roll4")
head(data)

# get row sums (total number of 6 rolls for each instance)
results <- apply(data, 1, sum)

# get the percent of instances that have exactly 2 target rolls
print(sprintf("simulation found %f%% of rolls with 2 sixes", mean(results == 2)))