Let's say you know the equation of a circle (r² = x² + y²), but want to calculate the area of a circle. And, uh, you forgot how to do that. Happens to all of us, honest.

Fortunately, we've got a little python program that uses the monte-carlo method to estimate the value of pi! Let's get going!

What we're going to do here is create a square and inscribe a circle inside it. Then we're going to start picking random points within the square -- points that are inside the circle will be red. The area outside the circle will be blue.

Since we know the area of a square (A = L²), how many points fell inside the circle, and how many fell outside, we can use that information to figure out the area of the circle, and from that, a value for pi.

First, you'll need to download and install VPython, a python module for easily displaying simple 3D graphics. Once you've got that all done with, you're ready to go.

I was going to write an explanation of what the code is doing and why, but the code with comments speaks for itself, I believe:

from visual import *
import random

# Written by Gary Arnold http://www.garyarnold.com
# Licensed under Creative Commons Attribution-Share Alike 3.0 Unported License
# http://creativecommons.org/licenses/by-sa/3.0/
# You may copy, distribute, transmit, and adapt this code
# You must attribute the original author (as the original author)
# You may distribute the resulting code only under the same, similar or a
#   compatible license

# Approximates PI by comparing the ratio of the area of a square with the area
#  of a circle inscribed within in.

# Radius of the circle we're considering
#  We're using 5 here to show there's no trickery in using 1 or anything
radius = 5.0

# To cut down on the math, we'll pre-compute R^2
radiusSquared = radius * radius

# Paint the "outside" area blue
box(pos=(0,0,0), length=2.0 * radius, height=2.0 * radius, \
    width=0, color=color.blue)

# Number of points inside the circle
numInside = 0

# Total number of random points we've considered
total = 0

# Text to display our calculations
ratioLabel = label(pos=(0,-radius - 0.2,0), \
    text='ratio = 0 / 0 = 0.0\nX = 4.0 * ratio = 0.0' )

# Just keep going forever...
while 1:
    # Collection of points that are inside the circle
    insidePoints = []

    # Process 1000 random points
    for p in xrange(1000):
        # Create random X,Y on the range [-radius,radius)
        x = random.uniform(-radius,radius)
        y = random.uniform(-radius,radius)

        # Check if X,Y coordinates are within the circle
        #  (Note that we don't bother with a square root...in this case,
        #   it makes no difference because we're just comparing)
        if ( x * x ) + ( y * y ) <= radiusSquared:
            # Point is inside circle -- add it to point collection and
            #  increment number of points inside circle
            insidePoints.append( (x, y) )
            numInside += 1.

    # We just processed 1000 more points, so increase the total
    total += 1000.

    # Create a new convex object using the points inside the circle
    #  (This is faster than adding more and more points to an existing object)
    convex(pos=insidePoints, color=color.red)

    # Calculate area of square
    # Side length of square:
    #   L = 2.0 * radius
    # Area of square:
    #   A = L * L
    #   A = 2.0 * radius * 2.0 * radius
    # Area of square is 4.0 * radius * radius
    areaOfSquare = 4.0 * radius * radius

    # Because we're talking about the area of a planar object, the area of
    #   the circle must be something of the form:
    # Some factor "X" * radius * radius
    # For a square, X = 4.0
    # Let's figure out what X is for a circle
    # The area of the circle is obviously smaller than the area of the square
    #  (Since we can see the circle is contained inside the square)
    # Let's use the ratio of (points inside) to (total points) to tell us
    #  what "X" is for a circle, like so:
    # Area of circle = ratio * areaOfSquare
    # Area of circle = X * radius * radius
    # So:
    # X * radius * radius = ratio * areaOfSquare
    # X * radius * radius = ratio * 4.0 * radius * radius
    # X                   = ratio * 4.0

    ratio =  numInside / total
    X = ratio * 4.0
    ratioLabel.text = 'ratio = %d / %d = %1.3f\nX = 4.0 * ratio = %1.3f' \
        % (numInside, total, ratio, X)

There you have it! Run this code by saving it to disk (there's a handy link at the bottom of this posting) and double-clicking the file or typing "python pi.py" in a console. As the circle fills in, you'll see the ratio slowly converge on 3.141, the first few digits of pi.

Plenty of other math and probability experiments are easily possible with python -- how many times do you need to roll a die before you roll a six, on average? How common is it to deal a straight flush in poker? What are fractals?

Maybe I'll try some of those out in a future post...

This blog entry is licensed under a Creative Commons Attribution-Share Alike 3.0 Unported License