Ramblings of a math and CS teacher

December 11, 2005

Python Monte Carlo Simulator (Very Simplistic!)

Filed under: Uncategorized — danschellenberg @ 7:14 pm

In teaching Math A30 here in Saskatchewan, there is a segment of the curriculum devoted to teaching students the basics of Data Analysis. Classic Monte Carlo simulations are touched upon, and students are given the opportunity to experience Monte Carlo simulations, whether through tossing dice, drawing cards, or any other “equally likely options” based event.

While some Monte Carlo simulation extensions can be made through the use of advanced calculators, I thought that a more interesting way for my students (and myself!) to experience simplistic Monte Carlo simulations would be to create a simple Python computer program that would automate the task for us. After a few hours, I came up with something that I think is quite useable. The result of one run of the program is seen below (formatting lost in HTML conversion, sorry!):

Trial # | Die Toss | # Tosses
1 | 1 2 6 6 3 5 6 6 2 4 | 10
2 | 6 3 1 5 3 1 5 1 4 1 5 6 1 6 2 | 15
3 | 3 6 3 5 6 1 3 4 1 4 5 3 5 5 3 2 | 16
4 | 6 1 3 1 1 2 4 5 | 8
5 | 6 1 2 2 2 5 2 4 6 2 1 2 4 5 2 2 1 2 6 5 4 4 1 5 3 | 25
6 | 1 5 1 1 2 4 4 2 4 2 5 5 5 1 2 4 6 5 1 5 5 1 4 6 6 4 3 | 27
7 | 3 6 1 2 2 4 1 1 5 | 9
8 | 4 5 3 1 4 2 6 | 7
9 | 3 1 3 2 5 5 6 6 6 4 | 10
10 | 1 4 4 1 6 2 5 6 2 3 | 10

Total tosses: 137
Number of trials: 10
Average tosses per trial: 13.7

This is a clone of the example given in Burt Thiessen’s Math A30 textbook, wherein he gives the example of 6 different hockey cards being placed in cereal boxes, and poses the question “How many boxes of cereal will you have to purchase until you have collected the entire set of six hockey cards?”.

Based on the output of the script from above, students would guess that on average, you would have to buy 14 boxes of cereal to be sure you would have all 6 hockey cards. However, taking it further, we could change the “trials” variable in the program to be 1000 instead of 10, and find that the average we came up with this time is 14.7. We therefore alter our original guess, now believing we would have to buy 15 boxes on average. We could continue increasing the number of trials until we are satisfied that we have come up with an acceptable conclusion.

A fairly straightforward exercise, but I think the students will appreciate this program, especially after having spent a few minutes attempting to do just a few trials on their own…

If you are interested in the Python source code, see below, or check out the plain text version.

Source code:

# Monte Carlo Simulation
# Dan Schellenberg      Dec 11, 2005

# set low and high values
low = 1
high = 6
trials = 10

################################################################################
# This section of code is used to format the output into nice columns
# Code from http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/302380

import re

LEFT = ‘left’
RIGHT = ‘right’
CENTER = ‘center’

def align(text, width=70, alignment=LEFT):
    ”’ Align the “text” using the given alignment, padding to the given
    width. Strip off any existing whitespace on the side being aligned to
    and pad with spaces (‘
‘) on the opposite side.

    Code from http://www.faqts.com/knowledge_base/view.phtml/aid/4476
    ”’
    if alignment == CENTER:
        text = text.strip()
        space = width – len(text)
        return ‘ ‘*(space/2) + text + ‘ ‘*(space/2 + space%2)
    elif alignment == RIGHT:
        text = text.rstrip()
        space = width – len(text)
        return ‘ ‘*space + text
    else:
        text = text.lstrip()
        space = width – len(text)
        return text + ‘ ‘*space

class FormatColumns:
    ”’Format some columns of text with constraints on the widths of the
    columns and the alignment of the text inside the columns.
    ”’

    def __init__(self, columns, contents, spacer=‘ | ‘, retain_newlines=True):
        ”’
        ”columns”   is a list of tuples (width in chars, alignment) where
                    alignment is one of LEFT, CENTER or RIGHT.
        ”contents”  is a list of chunks of text to format into each of the
                    columns.
        ”’

        assert len(columns) == len(contents),
            ‘columns and contents must be same length’
        self.columns = columns
        self.num_columns = len(columns)
        self.contents = contents
        self.spacer = spacer
        self.retain_newlines = retain_newlines
        self.positions = [0]*self.num_columns

    def format_line(self, wsre=re.compile(r’s+’)):
        ”’ Fill up a single row with data from the contents.
        ”’

        l = []
        data = False
        for i, (width, alignment) in enumerate(self.columns):
            content = self.contents[i]
            col =
            while self.positions[i] < len(content):
                word = content[self.positions[i]]
                # if we hit a newline, honor it
                if ‘n’ in word:
                    # chomp
                    self.positions[i] += 1
                    if self.retain_newlines:
                        break
                    word = word.strip()

                # make sure this word fits
                if col and len(word) + len(col) > width:
                    break

                # no whitespace at start-of-line
                if wsre.match(word) and not col:
                    # chomp
                    self.positions[i] += 1
                    continue

                col += word
                # chomp
                self.positions[i] += 1
            if col:
                data = True
            col = align(col, width, alignment)
            l.append(col)

        if data:
            return self.spacer.join(l).rstrip()
        # don’t return a blank line
        return

    def format(self, splitre=re.compile(r‘(n|rn|r|[ t]|S+)’)):
        # split the text into words, spaces/tabs and newlines
        for i, content in enumerate(self.contents):
            self.contents[i] = splitre.findall(content)

        # now process line by line
        l = []
        line = self.format_line()
        while line:
            l.append(line)
            line = self.format_line()
        return ‘n’.join(l)

    def __str__(self):
        return self.format()

def wrap(text, width=75, alignment=LEFT):
    return FormatColumns(((width, alignment),), [text])

################################################################################
# This is the actual Monte Carlo stuff…

from random import randint

def rollAll( low, high, trials ):
# Must “roll” each element in the list at least once…
    
    # print header
    print FormatColumns(((8, CENTER), (58, CENTER), (8, CENTER)), ["Trial #", "Die Toss", "# Tosses"])
    
    grand_total = 0
    grand_tosses = 0    
    
    i = 0
    while i < trials:
    #for each trial, each number must be rolled at least once
        arraySize = range(low, high+1)
        allChosen = []
    
        #populate list with all False elements
        for element in arraySize:
            allChosen.append(False)
        
        #initialize variables
        allDone = False
        tosses = 0
        total = 0
        numData = “”
        
        while allDone == False:
            num = randint(low, high)
            numData = numData + str(num) + ” “
            tosses = tosses + 1
            total = total + num
            
            # set the “this number rolled” flag to true
            allChosen[num - 1] = True
            
            # check if all “numbers” have been “rolled”
            for element in allChosen:
                if element == False:
                    allDone = False
                    break
                else:
                    allDone = True

        tossData = str(tosses)
        grand_total = grand_total + total
        grand_tosses = grand_tosses + tosses
        print FormatColumns(((8, CENTER), (58, LEFT), (8, CENTER)), [str(i+1), numData, tossData])
        i = i + 1
    
    #print results
    print
    print “Total tosses:”, grand_tosses
    print “Number of trials:”, trials
    print “Average tosses per trial:”, float(grand_tosses)/trials
    
    return 1

rollAll(low, high, trials)

Blog at WordPress.com.