rubymathprobabilitystandard-deviation

Probability of getting specific sum after rolling n dice. Ruby


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?


Solution

  • 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