CODEWITHSUNDEEP
pythonprobabilitydice
Ginger Leprechaun
Ginger Leprechaun

Reputation: 13

How to calculate the probability of the sum of n rolled m sided dice of being a certain value or greater

I'm trying to write some Python code which will calculate the probability of the sum of n rolled m sided fair dice of being a certain value or greater. The inputs are the sum of all dice rolls, the number of rolls made and the number of sides on the dice. The output is the percentage chance of the sum of dice rolls being that value or greater.

I'm basing my calculations on the equations I found in this paper and scaling it for any sided die being rolled any number of times: https://digitalscholarship.unlv.edu/cgi/viewcontent.cgi?article=1025&context=grrj

I've made some code which "works" but is incredibly slow for when the dice has lots of faces so is only useful for dice with 20 faces or less.

import numpy as np

def probability_calculator(roll_total, num_of_rolls, dice_faces):       

    if num_of_rolls == 0:
        probability = 100    
    elif num_of_rolls == 1:
        probability = (dice_faces + 1 - roll_total) * 100/dice_faces    
    else:
        inverse = (dice_faces ** num_of_rolls) ** -1
        side_list = np.linspace(1, dice_faces, dice_faces)
        expanded_list = np.zeros(dice_faces * num_of_rolls)         
        stacked_side_list = side_list   

        for i in range(num_of_rolls - 1):
            stacked_side_list = np.vstack((stacked_side_list, side_list))
            
        index_array = np.zeros(num_of_rolls, dtype=int)
            
        while True:                
            value = 0

            for i in range(num_of_rolls):                    
                value = value + stacked_side_list[i][index_array[i]]
                
           expanded_list[int(value) - 1] += 1

           if sum(index_array) == (dice_faces - 1) * num_of_rolls:
               break

           for i in range(num_of_rolls):  
               if index_array[i] == dice_faces - 1:
                   index_array[i] = 0          
               else:
                    index_array[i] += 1                        
                    break

       probability = inverse * sum(expanded_list[roll_total - 1:]) * 100
  
   return probability

As you can see, this is incredibly inefficient code which if you were to roll only four 100 sided die you would have to iterate through the while loop 100^4 = 100,000,000 times.....

I am pretty certain there is some mathematical equation which can simplify this code and make it run many orders of magnitude faster but maths isn't my favourite subject and I don't know of any equation or Python function that could help.

Upvotes: 1

Views: 1170

Answers (1)

no comment
no comment

Reputation: 10463

Had a quick look at the paper and got scared by the formulas. Implemented a presumably different way simply counting how often the different totals can be reached:

from functools import lru_cache

@lru_cache(None)
def sum_freq(total, rolls, faces):
        if not rolls:
            return not total
        return sum(sum_freq(total - die, rolls - 1, faces)
                   for die in range(1, faces + 1))

def probability_calculator(roll_total, num_of_rolls, dice_faces):
    return sum_freq(roll_total, num_of_rolls, dice_faces) / dice_faces**num_of_rolls

Demo for your "four 100 sided die":

prob_314 = probability_calculator(314, 4, 100)
prob_any = sum(probability_calculator(total, 4, 100)
               for total in range(1, 401))
print(f'{prob_314:%}')
print(f'{prob_any:%}')

Output:

0.113564%
100.000000%

Output for ten 100 sided dice:

0.050065%
100.000000%

Output for 42 such 100 sided dice:

0.000000%
100.000000%

Try it online!

Upvotes: 1

Related Questions