Wiki

Monte Carlo Pi

Did you know you can estimate the value of pi by throwing random points on a square?

alt text

Now draw a square around that circle, like this:

alt text

The area of the square is 4 = (1+1)*(1+1), and the area of the circle is = 3.1459...

If you throw darts at the square randomly, the number that land in the circle divided by the total darts thrown will approximate pi! Neat, huh? We use a quarter of the square and then multiply the estimate by 4. This is important, because the fraction of (darts in circle/darts thrown) will always be less than or equal to 1. To get pi, I start by using 1/4 of the square, and so multiplying by 4 brings us "full circle".

[Note: This is also one simple way to perform rejection sampling over the distribution defined by the area under pi/2 radians on the unit circle.]

The approximation approaches the true value of pi as the number of darts approaches infinity (i.e. you throw enough darts to cover the entire circle), and assuming truly random throws. A successful run should produce something like:

alt text

Pulling the Repo

You can pull this repo with:

$ git clone https://github.com/entangledloops/monte_carlo_pi.git

See the readme for further details on running this script on your machine.

Python Code

#!python

# Author:      Stephen Dunn
# Date:        8/16/13
# File:        pi.py
# Description: Monte-Carlo Pi Estimator
# Notes:       1) exceeding ~pi to 4 decimal places with anything better than luck has low probability
#              2) yes, this could make better use of pygame functions (e.g. not hand-computing the semicircle)

import sys
import pygame
import math
import random

window = ""
width =  800
height = 600
resolution = 2000 # precision of line segments
points = 10000 # of points

def graph():

  x = y = ex = ey = 0
  dp = (math.pi/2)/resolution

  # draw a quarter-circle
  for i in range(0,resolution):
    x = ex
    y = ey
    ex = ex+dp
    ey = ey+dp

    pygame.draw.aaline(window, (128,255,128), \
        ( math.cos(x)*width, math.sin(y)*height ), \
        ( math.cos(ex)*width, math.sin(ey)*height ), 3)

  # estimate pi from random points
  hits = 0
  for i in range(0,points):
    px = random.random()
    py = random.random()
    c = math.sqrt( (px*px) + (py*py) ) 
    if (c <= 1.0): hits += 1
    px *= width
    py *= height
    pygame.draw.circle( window, (148,0,211), (int(px),int(py)), 1, 1 )

  api = 4.0*(float(hits)/float(points)) # approximate pi
  font = pygame.font.SysFont('Arial', 18)
  pi_str = u'\u03c0' + ' ' + u'\u2248' + " " + "{0:.5f}".format(api)
  pi_str.encode('utf-8')
  print pi_str
  text = font.render(pi_str, True, (72,209,204))
  window.blit( text, (width/4,height/4) )

#----------------------------------
random.seed()
pygame.init()

window = pygame.display.set_mode( (width, height) )

print "graphing..."
graph()
print "done"

pygame.display.flip()
while True:
  for event in pygame.event.get():
    if event.type == pygame.QUIT:
      sys.exit(0)