Monte Carlo Approximation of Pi

Date: 11/24/2001 at 11:47:22
From: Sharmila V.
Subject: Monte Carlo approximation of pi

Using the circle-in-square model, find an algorithm to determine the 
value of pi. You can assume the existence of two functions - 

   uniform (a,b) - returns a uniformly distributed random variable of 
     type real between a and b.

   radis (x,y) - given the coordinates x and y, the function returns 
     the distance from the centre of the circle/square at which the   
     coordinates lie.

How many trials will it take to achieve a 90% confidence interval for 
pi? or a 95% confidence interval?

Thanks!

Date: 11/25/2001 at 08:39:29
From: Doctor Jubal
Subject: Re: Monte Carlo approximation of pi

Hi Sharmila,

Algorithms like these, which use random numbers to approximate 
deterministic outcomes, are called Monte Carlo methods, after the 
famous casino city where random numbers produce a deterministic 
outcome (the house wins the gambler's money).

The unit circle has an area of pi * r^2 = pi * 1^2 = pi, while the 
square it is inscribed in has side length 2, so its area is 4.

One way of of approximating pi is to pick a point at random inside the 
square. With probability pi/4, that point also lies inside the unit 
circle, because the unit circle has pi/4 the area of the square.

So this gives a way to approximate pi/4. Pick a large number of points 
at random, and then approximate

           number of points inside unit circle
  pi/4 ~ --------------------------------------
                 total number of points

Of course, you can find pi itself by multiplying by four at the end.

So to implement this scheme, you need a way of randomly picking a 
point inside the square and a way of testing whether that point is 
inside the unit circle.

To pick a random point inside the square, we need to know where the 
square is. Since the unit circle is centered at the origin and has 
radius one, the square it is inscribed in has edges at x = +1, 1 and 
y = +-1. So if we call uniform(-1,1) twice, once for the x-coordinate 
and once for the y-coordinate, we'll have a random point in the 
square.

The unit circle is the set of points that are exactly 1 unit away from 
the origin. Anything less than 1 unit away is inside it. Anything more 
than 1 unit away is outside it. So once you have the x and y 
coordinates, that coordinate pair lies inside the unit circle if 
radius(x,y) < 1.

Repeat a number of times, recording how many of the points lie inside 
the circle along with the total number of points, and you have an 
answer.

How good an answer depends on the number of points. For this, we'll 
need some statistics 

The usual measure of a random process's variability is the variance 
s^2.  s^2 is defined

  s^2 = <r^2> - <r>^2

where <r^2> is the expectation of the random variable r squared, and 
<r>^2 is the square of the random variable expectation (in other 
words, the square of the mean.

The random variable r in this case is the true/false property of 
whether or not our randomly chosen point is inside the unit circle.  
Let's assign the property of being inside the circle a value of 1, and 
being outside the circle a value of 0. With probability pi/4, the 
point is inside the circle, so the mean (expectation) of r is

  <r> = (pi/4)(1) + (1-pi/4)(0) = pi/4

Technically, we already knew this, since it's the whole basis for the 
algorithm, but it's nice to know we're on the right track.

The expectation of r^2 is also

  <r^2> = (pi/4)(1^2) + (1-pi/4)(0^2) = pi/4

so the variance s^2 is

  s^2 = <r^2> - <r>^2 = pi/4 - (pi/4)^2 = (pi/4)(1 - pi/4)

For our problem, the variance itself is only useful because we can use 
it to derive a couple of other more useful quantities, namely the 
standard deviation s, which is just the square root of the variance:

  s = sqrt[(pi/4)(1-pi/4)]

and the relative error, which is the standard deviation divided by the 
mean

  s/m = sqrt[(pi/4)(1-pi/4)] / (pi/4) = sqrt(4/pi - 1) = .5227232...

The relative error measures how big the "average" error is compared to 
the correct answer. In this case, the "average" error is more than 
half the size of the correct answer, which is downright terrible! But 
this is for only a single trial, which will give us one of only two 
answers, pi = 4 or pi = 0, so it's understandable why it's so bad.

The relative error decreases as the square root of the number of 
trials, that is (N is the number of trials)

  s/m = (s/m for one trial)/sqrt(N)
               4/pi - 1
  s/m = sqrt( ----------- )
                   N

So to find the 90% and 95% confidence intervals, you need to find the 
values of N that decrease the relative error to .1 and .05, 
respectively.

Does this help?  Write back if you'd like to talk about this some
more, or if you have any other questions.

- Doctor Jubal, The Math Forum
  http://mathforum.org/dr.math/   

Date: 11/25/2001 at 10:28:04
From: Sharmila Varshney
Subject: Re: Monte Carlo approximation of pi

Thanks for your quick reply to my question. The only thing is that it 
is a unit circle in a unit square. I hope this does not change the 
approach. I will test it out with the area of square equal to 2 
instead of 4.

I am using your forum for the first time and am very impressed. I 
think you all are doing a wonderful job!. I did not think that in this 
time and age people actually did this kind of thing. Thanks a million.

Sharmila

Date: 11/25/2001 at 10:43:04
From: Doctor Jubal
Subject: Re: Monte Carlo approximation of pi

Hi Sharmila,

That approach won't work, because the unit circle (which has area pi) 
won't fit inside a unit square (which has area 1, not 2). If you want 
to use a unit square, what you could do is take the arrangement of 
circle and square I described in my first response and consider only 
the part of it that lies in the first quadrant.

That arrangement gives a unit square (area 1), with one fourth of the 
unit circle (area pi/4) inside it. The ratio of area occupied by 
circle to area occupied by square is still pi/4, because each of the 
fourth quarters of the circle-square combination is completely 
identical except for being rotated 90 degrees with respect to its 
neighbors.

To do this, the only modification you'd need to make to the approach 
I suggested to first time is to call uniform(0,1) rather than 
uniform(-1,1). This gives you random points that are guaranteed to be 
in the first quadrant. Because of symmetry, the test to see if the 
point is inside the unit circle still works, and because the 
probability of being in the unit circle is still pi/4, the error 
analysis still holds.

Does this help?  Write back if you'd like to talk about this some
more, or if you have any other questions.

- Doctor Jubal, The Math Forum
  http://mathforum.org/dr.math/   

Date: 11/25/2001 at 10:50:14
From: Sharmila Varshney
Subject: Re: Monte Carlo approximation of pi

I was just realizing the problem when I received your second email. 
Now I get the whole picture. Thanks so much! You have no idea how much 
of a help you have been.

Would the approach be different if I had to consider a circle with 
diameter 1 in a unit square?

Sharmila

Date: 11/27/2001 at 09:59:27
From: Doctor Jubal
Subject: Re: Monte Carlo approximation of pi

Hi Sharmila,

Absolutely not. The shapes and the ratios of area are all the same.  
The only thing to change is that since the circle now has radius 1/2 
and the square extends from -1/2 < x or y < 1/2, you want to generate 
random points with uniform(-1/2, 1/2) and then test if they're inside 
the circle by testing if radius(x,y) < 1/2.

You can do this algorithm with any size square you want just by 
adjusting the constants. Most of the time, though, it's done using 
1/4 of a circle with radius 1, simply because the library random 
number generating routines on many platforms generate random numbers 
with a uniform(0,1) distribution.

- Doctor Jubal, The Math Forum
  http://mathforum.org/dr.math/