Probability of getting specific sum after rolling n dice. Ruby

Question

What is the best solution for finding probability of rolling sum with n dice? I'm solving it by finding

1. mean.
2. standard deviation.
3. z_score for the numbers below x
4. z_score for the numbers above x
5. converting both to probabilities
6. subtracting one from the other

This is what I've done so far.

# sides - number of sides on one die
def get_mean(sides)
  (1..sides).inject(:+) / sides.to_f
end

def get_variance(sides)
  mean_of_squares = ((1..sides).inject {|sum, side| sum + side ** 2}) / sides.to_f
  square_mean = get_mean(sides) ** 2

  mean_of_squares - square_mean
end

def get_sigma(variance)
  variance ** 0.5
end

# x - the number of points in question
def get_z_score(x, mean, sigma)
  (x - mean) / sigma.to_f
end

# Converts z_score to probability
def z_to_probability(z)
  return 0 if z < -6.5
  return 1 if z > 6.5

  fact_k = 1
  sum = 0
  term = 1
  k = 0

  loop_stop = Math.exp(-23)
  while term.abs > loop_stop do
    term = 0.3989422804 * ((-1)**k) * (z**k) / (2*k+1) / (2**k) * (z**(k+1)) / fact_k
    sum += term
    k += 1
    fact_k *= k
  end

  sum += 0.5
  1 - sum
end

# Calculate probability of getting 'х' total points by rolling 'n' dice with 'sides' number of sides.
def probability_of_sum(x, n, sides=6)

  mean = n * get_mean(sides)
  variance = get_variance(sides)
  sigma = get_sigma(n * variance)

  # Rolling below the sum
  z1 = get_z_score(x, mean, sigma)
  prob_1 = z_to_probability(z1)

  # Rolling above the sum
  z2 = get_z_score(x+1, mean, sigma)
  prob_2 = z_to_probability(z2)

  prob_1 - prob_2
end

# Run probability for 100 dice
puts probability_of_sum(400, 100)

What concerns me is, when I pick x = 200, the probability is 0. Is this correct solution?

Answer 1

There is an exact solution involving an alternating sum of binomial coefficients. I have written it out in a few places (on Quora and MSE), and you can find it elsewhere although there are some flawed versions. Be careful that if you implement that, you may need to cancel binomial coefficients that are much larger than the final result, and if you use floating point arithmetic you might lose too much precision.

Neil Slater's suggestion to use dynamic programming to compute the convolution is a good one. It is slower than the summation of binomial coefficients, but reasonably robust. You can speed it up in a few ways, such as using exponentiation by squaring, and by using the Fast Fourier Transform, but many people will find those to be overkill.

To fix your method, you should use the (simple) continuity correction to the normal approximation, and restrict to the context where you have enough dice and you are evaluating far enough from the maximum and minimum that you expect a normal approximation will be good, either in an absolute or relative sense. The continuity correction is that you replace the count of n with the interval from n-1/2 to n+1/2.

The exact count of the number of ways to roll a total of 200 is 7745954278770349845682110174816333221135826585518841002760, so the probability is that divided by 6^100, which is about 1.18563 x 10^-20.

The normal approximation with the simple continuity correction is Phi((200.5-350)/sqrt(3500/12))-Phi((199.5-350)/sqrt(3500/12)) = 4.2 x 10^-19. This is accurate in absolute terms, since it is very close to 0, but it is off by a factor of 35 so it's not great in relative terms. The normal approximation gives a better relative approximation closer to the center.