Probability that the sum of 10d4 beat the sum of 7d6

######################################################################
# Calculate the probability that the sum of 10d4 beat the sum of 7d6
# using convolutions to the PMFs and, hence, the joint PMF
######################################################################

#' Function to convolute two discrete probability distributions with
#' support on 0, 1, ... I.e. we compute the PMF of Z = X + Y.
#'
#' @param fX - the PMF of X given as a named vector, where the names
#' represent the support X_min, ..., 0, 1, ..., X_max
#' @param fY - the PMF of Y given as a named vector, where the names
#' represent the support Y_min, ..., 0, 1, ..., Y_max
#'
#' @return A named vector containing the PMF of Z = X+Y.
#'
convolute <- function(fX, fY) {
  fZ <- rep(0, 1 + (length(fX) - 1) + (length(fY) - 1))
  names(fZ) <- as.character(seq(
    min(as.numeric(names(fX))) + min(as.numeric(names(fY))),
    max(as.numeric(names(fX))) + max(as.numeric(names(fY)))
  ))
  for (i in names(fX)) {
    for (k in names(fY)) {
      j <- as.numeric(i) + as.numeric(k)
      fZ[as.character(j)] <- fZ[as.character(j)] + fX[i] * fY[k]
    }
  }
  fZ
}

# Function to compute the PMF of X = \sum_{i=1}^n X_i when the
# X_i are iid RV with PMF `pmf_one`.
convolute_n_times <- function(n, pmf_one) {
  i <- 1
  pmf <- pmf_one
  while (i<n) {
    i <- i+1
    pmf <- convolute(pmf, pmf_one)
  }
  return(pmf)
}

# Compute PMF for the sum of 10d4, i.e. X=\sum_{i=1}^n X_i with X_i being
# iid U(\{1,2,3,4}) RVs.
pmf_d4 <- c(0,rep(1/4, 4)) %>% setNames(0:4)
pmf_10d4 <- convolute_n_times(n=10, pmf_d4)
plot(pmf_10d4, type="h", xlab="x", ylab=expression(f[X](x)))
# Calculate expectation of X
10 * sum(0:4 * pmf_d4)
sum(0:40 * pmf_10d4)

# Compute PMF for the sum of 7d6, i.e. Y=\sum_{i=1}^n Y_i with Y_i being
# iid U(\{1,2,3,4,5,6}) RVs.
pmf_d6 <- c(0,rep(1/6, 6)) %>% setNames(0:6)
pmf_7d6 <- convolute_n_times(n=7, pmf_d6)
plot(pmf_7d6, type="h", xlab="y", ylab=expression(f[Y](y)))
# Expectations
7*sum(0:6 * pmf_d6)
sum(0:(7*6) * pmf_7d6)

# Calculate the joint probability of X and Y (which are independent)
pmf_XY <- outer(pmf_10d4, pmf_7d6, "*")
image(x=as.numeric(names(pmf_10d4)), y=as.numeric(names(pmf_7d6)), pmf_XY, xlab="x", ylab="y")
# Rows correspond to value of X, columns correspond to columns of Y
rownames(pmf_XY)
colnames(pmf_XY)
# Indicator function I(x>y)
d4_win <- (row(pmf_XY)-1) > (col(pmf_XY)-1)
# P(X>Y) = E(I(X>Y))
# i.e. P(X>Y) = \sum_{x=0}^{40} \sum_{y=0}^{42} I(x>y) f_{X}(x) f_{Y}(y).
sum(d4_win * pmf_XY)

# Alternative:
# https://topps.diku.dk/torbenm/troll.msp
# # Variation of Euler problem 205 with dice numbers yielding a prob of 1/2.
# count ((sum 10d4) > (sum 7d6))