Simulation of 2 dice rolls

I wrote a program that records how many times 2 fair dice need to be rolled to match the probabilities for each result that we should expect.

I think it works but I'm wondering if there's a more resource-friendly way to solve this problem.

import random

expected = [0.0, 0.0, 0.028, 0.056, 0.083, 
         0.111, 0.139, 0.167, 0.139, 0.111,
         0.083, 0.056, 0.028]

results = [0.0] * 13  # store our empirical results here

emp_percent = [0.0] * 13  # results / by count

count = 0.0  # how many times have we rolled the dice? 

while True:
    r = random.randrange(1,7) + random.randrange(1,7)  # roll our die
    count += 1 
    results[r] += 1
    emp_percent = results[:]

    for i in range(len(emp_percent)):
        emp_percent[i] /= count
        emp_percent[i] = round(emp_percent[i], 3)

    if emp_percent == expected:
        break

print(count)
print(emp_percent)

I question the validity of this exercise. I wouldn't really expect the percentages to line up for a long, long, long time. Possibly never. Running a few times I got anywhere from 290k to nearly 8 million. That's a lot of die rolls, and a pretty wide range.

That said, we can certainly improve your code. First, there's random.randint() which is more directly what you want. Let's just put that in a function:

def die_roll():
    return random.randint(1, 6)

Now we don't need that weird 7.

Next, results should be ints. After all, we're iterating upward. We only need to convert them to floats when we do the calculation. We can also use itertools.count() for the loop:

results = [0] * 13
for count in itertools.count(start=1):
    results[die_roll() + die_roll()] += 1

Ok, now that we have our results, we need to compare against expected. Rather than doing a list copy, division, round on every element - let's go one at a time. After all, if the percentage for 2 doesn't line up, why even calculate the percentage for 3?

The percentage for any one roll is:

round(1. * r / count, 3)

for r in results. And the expectation is e in expected. We want to iterate over both lists simultaneously... "zipping" through them. For that, there's itertools.izip():

if all(round(1. * r/count, 3) == e
       for r, e in itertools.izip(results, expected)):
    break

Lastly, don't hardcode your percentages! These have well defined values. We have a pyramid (0, 1, 2, 3, 4, 5, 6, 5, 4, 3, 2, 1) that we just multiply by 1/36. So that's:

def prob(x):
    count = x-1 if x <= 7 else 13-x
    return round(count/36., 3)

expected = [0.0] + [prob(x) for x in range(1, 13)]

So the full solution is:

import random
import itertools

def prob(x):
    count = x-1 if x <= 7 else 13-x
    return round(count/36., 3)

expected = [0.0] + [prob(x) for x in range(1, 13)]

results = [0] * 13

def die_roll():
    return random.randint(1, 6)

for count in itertools.count(start=1):
    results[die_roll() + die_roll()] += 1

    if all(round(1. * r/count, 3) == e
           for r, e in itertools.izip(results, expected)):
        break

print count

Dropping Floats

Floating point arithmetic is always slow. We could do a bit better by keeping everything integral:

import random
import itertools

def prob(x):
    count = x-1 if x <= 7 else 13-x
    return 1000*count/36

expected = [0] + [prob(x) for x in range(1, 13)]
results = [0] * 13

def die_roll():
    return random.randint(1, 6)

for count in itertools.count(start=1):
    results[die_roll() + die_roll()] += 1

    if all(1000 * r/count == e
           for r, e in itertools.izip(results, expected)):
        break

print count

Generalize arbitrarily

Sometimes, the people want an over-engineered solution. So let's give them an over-engineered solution! We can use itertools.product() to help us determine the odds of particular die rolls. Everything else is straightforward:

def roll_until_match(sides=6, num_dice=2):
    c = collections.Counter(sum(p) for p in
                            itertools.product(xrange(1, sides+1), repeat=num_dice))

    expected = [1000 * c.get(i, 0) / (sides ** num_dice)
                for i in xrange(sides * num_dice + 1)] 

    results = [0] * (sides * num_dice + 1)

    def die_rolls():
        return sum(random.randint(1, sides) for _ in xrange(num_dice))

    for count in itertools.count(start=1):
        results[die_rolls()] += 1

        if all(1000 * r/count == e
               for r, e in itertools.izip(results, expected)):
            break

    print count
    print results

StackExchange Code Review Q#106167, answer score: 6