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/
Search the Dr. Math Library:
Ask Dr. MathTM
© 1994- The Math Forum at NCTM. All rights reserved.