March 14 is Pi Day

March 14 is Pi Day, so let’s use a Monte Carlo method to estimate π.

See:

One Monte Carlo method for computing π is to draw a circle inscribed in a square, and randomly place dots in the square. The ratio of dots inside the circle to the total number of dots will approximately equal π/4.

– quoted (approximately) from the first Wikipedia article listed above.

The code below uses that type of Monte Carlo method to estimate π, when executed in Processing Python Mode. Points are marked in the sketch as they are chosen randomly. This implementation has obvious limitations, with one being the finite number of pixels in the sketch. More effective methods would use genuine floating point mathematics instead of relying on an image. This is offered here primarily as a celebratory demo. :fireworks:

# A Monte Carlo method for computing π 
# is to draw a circle inscribed in a square,
# and randomly place dots in the square.
# The ratio of dots inside the circle
# to the total number of dots
# will approximately equal π/4.
# - quoted from Wikipedia: Pi
# https://en.wikipedia.org/wiki/Pi

hits = 0.0 # current number of random points within the circle
import random
def setup():
    size(512, 552)
    
    # draw a blue circle inscribed in a red square
    rectMode(CENTER)
    ellipseMode(CENTER)
    textAlign(CENTER, CENTER)
    textSize(32)
    background(255)
    noSmooth()
    stroke(255, 0, 0)
    fill(255, 0, 0)
    rect(255, 255, 512, 512)
    stroke(0, 0, 255)
    fill(0, 0, 255)
    ellipse(255, 255, 512, 512)
    
def draw():
    global hits
    # choose a random pixel within bounds of the square
    x = random.randint(0, 511)
    y = random.randint(0, 511)
    pix = get(x, y)
    if blue(pix) == 255:
        # the pixel is blue so the point is in the circle
        hits += 1
        # mark point; maintain blue component
        stroke(0, 255, 255)
        point(x, y)
    else:
        # the pixel is not blue so the point is outside the circle
        # mark point; maintain red component
        stroke(255, 255, 0)
        point(x, y)
    # report the current estimate every 60 frames
    if frameCount % 60 == 0:
        stroke(255)
        fill(255)
        rect(width / 2, height - 20, width, 40)
        fill(0)
        # estimated pi is 4.0 * ratio of hits to total number of points
        # frameCount equals the total number of points
        text(str(4.0 * (hits / frameCount)), width / 2, height - 20)


Image capture of π estimation in progress.

Edited on March 14, 2021 to make minor adjustments to the Python code and to update the image.

6 Likes