dice roll math with large n (>100)

Question

I promise this is not just another dice rolling homework problem. I implemented a function to calculate the probability of obtaining less than a sum s when rolling n m-sided dice. My function works for small values of n but I am finding weird results for large values of n. See the attached plot. Anyone have insight into what is going on?

My probability function

Implemented from this math stack exchange

probability <- function(s, m, n) {

  i <- 0:((s-1-n) / m)
  m^(-n) * sum((-1)^i * choose(n, i) * choose(s - 1 - i * m, n))

}

Starts breaking with ~ n > 80

n <- 1:90 # number of dice
m <- 6 # number of sides
s <- floor(mean(1:m)*n) # sum of faces
p <- mapply(probability, s = s, m = m, n = n)
plot(n, p, main = paste("probability of rolling less than floor(", mean(1:m),"* n) as sum of n rolls of a", m, "sided dice"))

Answer 1

As mentioned in the comments on the origianl question, the problem is that the probability function is asking R to calculate really huge numbers (choose(80,40) = 1.075072e+23) and we are hitting the numerical precision limits of R.

An alternative approach that doesn't involve huge numbers but instead uses lots of numbers is to run monte carlo simulations. This generates a distribution of dice roll sums and compare the observed sum to the distribution. It will take longer to run, but is much easier to do and will not have the numerical precision problems.

mc <- Vectorize(function(s, m, n, reps = 10000) {
  x <- replicate(reps, sum(sample(m, n, replace = TRUE)))
  ecdf(x)(s-1)
})



n <- 1:90 # number of dice
m <- 6 # number of sides
s <- floor(mean(1:m)*n) # sum of faces
analytic_prob <- mapply(probability, s = s, m = m, n = n)
mc_prob <- mapply(mc, s = s, m = m, n = n)


plot(n, analytic_prob, main = paste("probability of rolling less than floor(", mean(1:m),"* n) as sum of n rolls of a", m, "sided dice"),
     sub = "monte carlo in red")
points(n, mc_prob, col = "red")