The probability of getting a certain amount after rolling n cubes. Ruby

Question

What is the best solution to find the probability of rolling a sum with n cubes? I solve it by finding

• mean.
• standard deviation.
• z_score for numbers below x
• z_score for numbers above x
• conversion as to probabilities
• subtracting one from the other

This is what I have 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) 

For me, when I choose x = 200 , the probability is 0. Is this the right solution?

Answer 1

Here is my final version.

• Changed offset amounts in get_z_score to x-0.5 and x+0.5 respectively, for a more accurate result.
• Added return 0 if x < n || x > n * sides return 0 if x < n || x > n * sides to cover cases where the sum is less than the number of cubes and higher than the number of dice times the number of sides.
• Added test tests with the 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 

Answer 2

Adding the result of two independent probability distributions coincides with the folding of two distributions. If the distributions are discrete, then this is a discrete convolution.

So, if one stamp is presented as:

 probs_1d6 = Array.new(6) { Rational(1,6) } 

Then 2d6 can be calculated as follows:

 probs_2d6 = [] probs_1d6.each_with_index do |prob_a,i_a| probs_1d6.each_with_index do |prob_b,i_b| probs_2d6[i_a + i_b] = ( probs_2d6[i_a + i_b] || 0 ) + prob_a * prob_b end end probs_2d6 # => [(1/36), (1/18), (1/12), (1/9), (5/36), (1/6), # (5/36), (1/9), (1/12), (1/18), (1/36)] 

Although this is an n-square for the sides of the cubes, a logical combination can reduce this, making it generally less flexible for more complex settings. The best part about this approach is that you can continue to add more cubes and make other more exotic combinations. For example, to get 4d6, you can drill two results for 2d6. Using smart solutions allows you to solve floating point problems.

I missed one detail, you need to save the initial offset (+1 for a normal six-sided matrix) and add it together to find out what the probabilities are connected with.

I made a more complex version of this logic in gem games_dice in a floating point, and not in Rational, which can handle several other combinations of dice.

Here's a basic rewrite of your method using the above approach in a naive way (just combining the effects of the cube one at a time):

 def probability_of_sum(x, n, sides=6) return 0 if x < n single_die_probs = Array.new(sides) { Rational(1,sides) } combined_probs = [1] # Represents chance of getting 0 when rolling 0 dice :-) # This is not the most efficient way to do this, but easier to understand n.times do start_probs = combined_probs combined_probs = [] start_probs.each_with_index do |prob_a,i_a| single_die_probs.each_with_index do |prob_b,i_b| combined_probs[i_a + i_b] = ( combined_probs[i_a + i_b] || 0 ) + prob_a * prob_b end end end combined_probs[ x - n ] || 0 end puts probability_of_sum(400, 100).to_f # => 0.0003172139126369326 

Note that the method actually calculates the full probability distribution from 100 to 600, so you only need to call it once and save the array (plus offset +100) once, and you can do other useful things, such as the probability of getting more than a certain number . All with great accuracy due to the use of Rational numbers in Ruby.

Because in your situation you have only one type of cubes, we can avoid using Rational to the end, working only with integers (essentially counting combined values) and divide by the total number of combinations (sides to the power of the number of rolls). This is much faster and returns values ​​for 100 dice per second:

 def probability_of_sum(x, n, sides=6) return 0 if x < n combined_probs = [1] # Represents chance of getting 0 when rolling 0 dice :-) n.times do start_probs = combined_probs combined_probs = [] start_probs.each_with_index do |prob_a,i_a| sides.times do |i_b| combined_probs[i_a + i_b] = ( combined_probs[i_a + i_b] || 0 ) + prob_a end end end Rational( combined_probs[ x - n ] || 0, sides ** n ) end 

Answer 3

There is an exact solution including a variable sum of binomial coefficients. I wrote this in several places (on Quora and MSE ), and you can find it elsewhere, although there are some erroneous versions. Be careful if you implement this, you may need to cancel binomial coefficients that are much larger than the final result, and if you use floating point arithmetic, you may lose too much accuracy.

Neil Slater's suggestion of using dynamic programming to compute convolution is a good one. It is slower than adding binomial coefficients, but reliable enough. You can speed it up in several ways, for example, by squaring and using the fast Fourier transform, but many people will find that they will be redundant.

To fix your method, you should use a (simple) continuity correction in the normal approximation and confine yourself to a context in which you have enough bones, and you estimate far enough from the maximum and minimum that you expect from the normal approximation to be good, or in absolute, either in a relative sense. The correction for continuity is that you replace counter n with an interval from n-1/2 to n + 1/2.

An exact calculation of the number of rolling methods of only 200 is 7745954278770349845682110174816333221135826585518841002760, so the probability is divided by 6 ^ 100, which is about 1.188563 x 10 ^ -20.

A normal approximation with 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 turned off 35 times, so it is not large in relative terms. A normal approximation gives a better relative approximation closer to the center.

Answer 4

I also decided that, using the method of Monte Carlo , and the results are relatively close.

 # x - sum of points to find probability for # n - number of dice # trials - number of trials def monte_carlo(x, n, trials=10000) pos = 0 trials.times do sum = n.times.inject(0) { |sum| sum + rand(1..6) } pos += 1 if sum == x end pos / trials.to_f end puts monte_carlo(300, 100, 30000)