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

Here is my final version.

1. Changed sum offsets in get_z_score to x-0.5 and x+0.5 respectively for more precise result.
2. Added return 0 if x < n || x > n * sides to cover cases, where sum is less then the number of dice and higher then number of dice multiplied by number of sides.
3. Added Benchmark tests with result

Main Functions

# 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)
  return 0 if x < n || x > n * sides

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

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

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

  prob_1 - prob_2
end

Benchmark

require 'benchmark'

Benchmark.bm do |x|
  x.report { @prob = probability_of_sum(350, 100).to_f }
  puts "\tWith x = 350 and n = 100:"
  puts "\tProbability: #{@prob}"
end
puts

Benchmark.bm do |x|
  x.report { @prob = probability_of_sum(400, 100).to_f }
  puts "\tWith x = 400 and n = 100:"
  puts "\tProbability: #{@prob}"
end
puts

Benchmark.bm do |x|
  x.report { @prob = probability_of_sum(1000, 300).to_f }
  puts "\tWith x = 1000 and n = 300:"
  puts "\tProbability: #{@prob}"
end

Result

       user     system      total        real
   0.000000   0.000000   0.000000 (  0.000049)
    With x = 350 and n = 100:
    Probability: 0.023356331366255034

       user     system      total        real
   0.000000   0.000000   0.000000 (  0.000049)
    With x = 400 and n = 100:
    Probability: 0.00032186531055478085

       user     system      total        real
   0.000000   0.000000   0.000000 (  0.000032)
    With x = 1000 and n = 300:
    Probability: 0.003232390001131513