# 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)

• When you say efficient, are you concerned with execution time or having clean code? – SuperBiasedMan Sep 30 '15 at 18:04
• Execution time. – Wrddot Sep 30 '15 at 18:09
• Why would you expect the percentages to be exactly equal ever? That's not how probability works... – Barry Sep 30 '15 at 18:15
• They can be equal, particularly since these are rounded, but it will certainly take a long time. – SuperBiasedMan Sep 30 '15 at 18:17

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

• Now you only have to resolve the hard-coded 7, 13, and 36 ;) What if we rolled three dice? Or what if we used twenty-sided dice? – mkrieger1 Sep 30 '15 at 21:32
• @mkrieger1 Ask and ye shall receive. – Barry Sep 30 '15 at 21:56

For performance, your huge bottleneck is that you make a new list as a copy of results every single time. This is expensive and wasteful. Instead, you could just update the number that you've just run a diceroll for. It's pretty simple, you can index emp_percent like you index results, and then run your calculation just on that index. So you're using the same calculation but only applying it to a single number and leaving the rest of the list unchanged.

while True:
r = random.randrange(1,7) + random.randrange(1,7)
count += 1
results[r] += 1
emp_percent[r] = round(results[r] / count, 3)

if emp_percent == expected:
break


Some other notes about your initialisation. results doesn't need to have floats, since you're just incrementing it by 1 each time. Later, when you copy this for emp_percent and divide by count as a float, you'll get the result as a float anyway. I also think you don't need the comment explaining it. results is a clear name. (same with count). percent_results is a clearer name, it matches with results and doesn't contain a confusing shortened word.

results = [0] * 13
percent_results = [0] * 13
count = 0


Also since you have the parameters in advance, you could just use while percent_results != expected instead of using break.