cheating_royal_riddle.R

# "Cheating royal riddle" from https://www.youtube.com/watch?v=hk9c7sJ08Bg

# Die designations.
x <- c(2, 7, 7, 12, 12, 17)
y <- c(3, 8, 8, 13, 13, 18)

# Tabulate the outcomes to get counts per face per die. Note that tabulate
# discards non-positive values, but this is fine since all of the values of
# interest are positive.
tx <- tabulate(x)
ty <- tabulate(y)

# Normalize to get a probability distribution per roll outcome. Note that the
# first element is the probability of rolling a 1 (not a zero).
tx <- tx/sum(tx)
ty <- ty/sum(ty)

# Convolve the two dice to get the probabilities _per roll_. Note that this uses
# fft, so the results are inexact (in particular, values of zero may be _near_
# zero in the output). Because we're adding 2 vectors each with a min index of
# 1, the value of the first index is actually for a result of 2.
roll <- convolve(tx, rev(ty), type="o")
min_value <- 2

num_rolls <- 20

# We now have the distribution of outcomes for a _single_ roll. Roll an
# additional 19 times to get the final distribution.
result <- roll
for (i in seq(2, num_rolls)) {
    result <- convolve(result, rev(roll), type="o")
}

# Value of the first index due to convolving with a 2-based initial index.
min_value <- 2 * num_rolls

df <- tibble(score=seq_along(result) + min_value - 1,
             p=result,
             cdf=cumsum(p))

# View the PDF
df %>% ggplot(aes(score, p)) + geom_col(width=1)
# View the CDF
df %>% ggplot(aes(score, cdf)) + geom_col(width=1)