Tuesday, January 06, 2015

Monte Carlo simulations to determine the value of PI

Consider a circle with radius R. The area of the circle is PI x R^2. The area of the square that just embeds our circle (as in image below) is (2R)^2.


One takes a print of the image, sticks it on a cardboard and throws a dart at it. The chance the dart will fall within the circle is the ratio of area of the circle to the area of the square.

Chance of dart landing within the circle =  Area of circle 
                                                                    Area of square

                                                                 =  PI x R^2 
                                                                       (2R)^2

                                                                 =  PI x R^2 
                                                                      4 x R^2

4 x (Chance of dart landing within the circle) = PI

Whether the dart is within the circle or beyond the circle is simply determined by the distance at which it is from the centre of the figure (as illustrated).



To determine the chance that a dart thrown at the figure will land within the circle, one throws the dart several times and notes the distance of its landing from the origin. 

The ratio  number of times the dart lands at a distance less than equal to the radius of the circle 
                                                              total number of throws

gives us the probability (i.e. chance) of the dart landing within the circle in total number of throws.

There is no bias expected. Hence, the probability of the dart falling within the quarter circle shown in the red square above is the same as that of the dart falling within the entire circle.

Programatically one now needs to randomly select (x, y) to represent a point at which the dart lands, determine if it is at a distance greater than R or not. If it was a unit circle (i.e. R = 1) the comparisons are a little easier since there is no need to take a root of (x^2 + y^2). 

Below is a program that gives us the probability of selecting a point (x, y) such that it falls in the quarter circle:


#include <iostream>
#include <ctime>

void usage()
{
std::cout << "Usage: find_pi <number of tries>" << std::endl;
return;
}

int main(int argc, char* argv[])
{
if (argc != 2)
{
usage();
return -1;
}

int numTries = atoi(argv[1]);

int total = 0;
int inCircle = 0;

srand(static_cast<unsigned>(time(0)));

for (int total = 0; total < numTries; total++)
{
double x = static_cast<double>(rand()) / RAND_MAX;
double y = static_cast<double>(rand()) / RAND_MAX;
double z = x*x + y*y;

if ((z - 1.0) <= DBL_EPSILON)
{
++inCircle;
}
}

std::cout << (double)inCircle * 4.0 / numTries << std::endl;

return 0;
}

No comments: