Monte Carlo integration

I recently came across this question: "You have a function called random that randomly generates a number between 0 to 1. Use this to calculate pi." Worded differently, the question is asking you to estimate pi using random numbers. As you can read in the script for Life of Pi (one of my favourite movies):

Pi, the sixteenth letter of the Greek alphabet, which is also used in mathematics to represent the ratio of any circle's circumference to its diameter - an irrational number of infinite length, usually rounded to three digits, as 3.14.

I later learned that the method used to answer this question is called Monte Carlo integration. The idea is that you randomly generate a pair of numbers, x and y, which represents a point on the Cartesian coordinate system and check whether this point lies inside or outside of a circle drawn on the same coordinate system. If we use random numbers between 0 and 1, the circle will have a radius of 1 and will be a unit circle.

How do we know whether a point lies inside or outside of the unit circle? We can use the Pythagorean theorem, which you may recall as x^2 + y^2 = z^2 . The hypotenuse in our case is the radius of the circle, so therefore if \sqrt{x^2 + y^2} \leq 1 , then a point lies within the circle. By randomly generating lots of points we can calculate the area of the circle by counting the number of points inside the circle divided by the total number of points and multiplying this number by 4. We multiply by 4 because when we use points ranging from 0 to 1 and a unit circle, we are only calculating the area in one quadrant.

But why is the area of unit circle an estimation of pi? Again as you may recall the area of a circle is \pi r^2 and when the radius is 1, the area is equal to pi. Below is the Perl code (pi.pl) to approximate pi.

#!/usr/bin/env perl

use strict;
use warnings;

my $usage = "Usage: $0 <number> <random_seed>\n";
my $n = shift or die $usage;
my $seed = shift or die $usage;

srand($seed);
my $inside = 0;

for (my $i = 0; $i < $n; ++$i){
    my $x = rand(1);
    my $y = rand(1);
    if (sqrt($x**2 + $y**2) <= 1){
        ++$inside;
    }
}

printf("%.5f\n", $inside / $n * 4);

The approximation improves as we generate more random numbers.

for n in 1000 10000 100000 1000000; do ./pi.pl $n $n; done
3.16400
3.15400
3.13592
3.14102

I have other examples of using Monte Carlo methods to solve problems.




Creative Commons License
This work is licensed under a Creative Commons
Attribution 4.0 International License
.

Leave a Reply

Your email address will not be published. Required fields are marked *

This site uses Akismet to reduce spam. Learn how your comment data is processed.