Using Monte Carlo method to find a specific probability of a certain dice output in Python 3.2

Question

I'm in a statistics class and we are constantly being given "dice problems" with various constrains. This is a probability problem, where I need to evaluate the probability of an event by using the Monte Carlo method. I know I could integrate my way through it but I'd like to write a programme that allows me to simply modify constraints in terms of how many dice I have, how many rolls I made and how many sides do those dice have.

This is one of the problems I'm working on. Simpler, cause I don't want to get carried away with my code writing "abilities".

Assume that a dice have 9 sides each. Estimate the probability that when you roll 5 dice at least 3 will have the same value.

This is the general template for the problems we are given: Roll an X number of n-­‐sided dice, each of which has sides numbered from 1 to n. Estimate the probability that we will get 3 or more dice with the same outcome.

I want to write a function for the sample problem, lets say, which takes as input an integer n, which is the number of faces in each die, and computes the probability that 3 or more dice have the same value. My biggest problem is the constraint "at least 3 out of 5". I've looked through other similar problems on Stackoverflow, but none of them really touch base with me. How would you write the code for constraint? I am using Python 3.2.

class Die(object):
  def __init__(self, sides = 9):
    self.sides = sides

  def roll(self):
    return randint(1, self.sides)

I'm stuck here. Any input is helpful, thanks!

Answer 1

I don't think I would use a class here. you just generate the dice roll and then check to see if the dice roll should be counted or not. In this case, I'd use a Counter to do the counting just to make the code a little more clean:

from collections import Counter
from random import randint

def roll(ndice,nsides=9):
    return [randint(1,nsides) for _ in range(ndice)]

def count_it():
    c = Counter(roll(5))
    return c.most_common(1)[0][1] >= 3

ntries = 100000
print (sum(1 for _ in range(ntries) if count_it())/ntries)

It looks to me like you have about a 10% chance. The trick with Monte-Carlo is determining whether you've converged or not. You can do this a few times with different numbers of ntries. The bigger your make ntries, the smaller the spread will be in your output. Eventually, when the spread is small enough, you say you've converged on the solution with some certainty.

Answer 2

Just use your Die class a large number of times:

# roll a lot of dice!
myDie = Die(9) # 9 sides

roll_counts = {side:0 for side in range(1, myDie.sides + 1)} 

numRolls = int(1e6)
for x in xrange(numRolls):
    roll_counts[myDie.roll()] += 1

Then analyze your distribution as needed:

for side in sorted(roll_counts):
    side_pct = float(roll_counts[side]) / numRolls * 100
    print 'side {} comprised {}% of all rolls'.format(side, side_pct)

EDIT: I'm aware this solution doesn't solve your homework problem, but hopefully it gives you the tools you need to roll and count your dice. It's likely you'll need to do the above for multiple die at a time and have a way of comparing their equality at each roll.