Simulations

Let's say that we wanted to know the probability that if we roll a pair of dice, the sum will be exactly 7. One way to do this (the way we'll do it) is to run a large number of simulations. In each, we'll simulate the roll of two dice. If we then divide the number of simulations that result in a sum of 7 by the total number of simulations,  we'll have an estimate of our probability.

How many simulations should we run? Thousands for use. Perhaps hundreds of thousands or even millions. How do we simulate the roll of a pair of dice? By use of the randint function is the random module. We'll simply have it generate random integers in the range from 1 to 6.

Simulations of this kind, ones in which we generate a large number of random values to estimate a certain probability, are called Monte Carlo simulations.

Here's my code.

import random

def sum_to_7(num_simulations):

    sums_of_7 = 0

    for _ in range(num_simulations):

        die_1 = random.randint(1, 6)

        die_2 = random.randint(1, 6)

        if die_1 + die_2 == 7:

            sums_of_7 += 1

    return sums_of_7 / num_simulations

A few comments:

1. The number of simulations is a parameter. This allows us to experiment with how the number of simulations effects the estimated probability.

2. We loop of course, once for each simulation. In each simulation, we "roll two dice"; that is, we generate two random integers in the range 1 - 6.

3. If in one simulation, the sum of our two dice is 7, we increment sums_of_7.

4. When the for loop is over, we return the number of successes divided by the total number of simulations.

Here's my output for a few inputs. As I said, the input is the number of simulations run.

>>> sum_to_7(10000)

0.1711

>>> sum_to_7(100000)

0.16636

>>> sum_to_7(1000000)

0.166292

Note that we should not expect that any of these outputs are precisely the probability that we'll get a sum of 7. But on the assumption that Python's randint function really is random and given that we did run a large number of simulations, it is very likely to be very close. (The exact probability is 1/6, and our estimates are indeed very close to that.)

Simulation One: Generalize

Generalize the above simulation in this way: accept the number of dice and the sum that counts as a success as inputs. (Above, the success sum was 7.) Call the function dice_sum. It's first line might be:

def dice_sum(num_dice, success_sum, num_simulations):

Here you'll want a for loop inside a for loop. The outer will begin with something like:

for _ in range(num_simulations):

The inner loop will simulate a roll of the given number of dice and so will iterate once for each die. This means it will begin with something like:

for die in range(num_dice):

A few test cases:

>>> dice_sum(2, 7, 10000)

0.1677

>>> dice_sum(3, 12, 10000)

0.1127

Simulation Two: 6's

Use Monte Carlo simulations to estimate the probability of these events:

1. Six fair dice are tossed and at least one 6 appears.

2. Twelve fair dice are tossed and at least two 6's appear.

3. Eighteen fair dice are tossed and at least three 6's appear.

As above, each function should accept the number of simulations that will be run as an input.

Simulation Three: Birthdays

In a group of two people chosen at random, it's of course quite unlikely that the two have the same birthday. In a group of three people chosen at random, it's a little less unlikely at least two will have the same birthday. As we increase the size of the group, the probability that at least two have the same birthday increases.

A question naturally arises: What's the smallest group size that will give us a probability of at least 0.5 that at least two have the same birthday?

Use a Monte Carlo simulation to estimate this probability.

A few suggestions:

1. As above, you'll have a loop that begins as for _ in range(num_simulations).

2. You'll also have an inner while loop. What will it do? Add birthdays to a list of birthdays until a birthday appears twice. That means you'll initialize a list-type variable (perhaps birthdays) to the empty list before you begin the while loop; and in the while loop, you'll add in each new birthday to the list of birthdays. The while loop will terminate when the most recent birthday is already in the list.

3. If birthdays is a list of birthdays and birthday is a particular birthday, birthdays.append(birthday) will add that birthday to the list.

4. We can test whether an item is in a list with a Boolean expression of the form an_item in a_list .

I'll let you think about how to use random.randint() to pick birthdays at random. There's a simple way to do that.