Saturday, July 26, 2008

Estimating the value of pi

Here's something I learned recently and found quite interesting. A neat way to calculate the value of pi using estimation is to take a unit circle inside a square of side 2 units and take random shots on the surface of the square.

unit circle in square








Now, the ratio of the area between the circle and square can be given by,

pi*r^2/L^2 = hits/shots




Thus, we can estimate the value of pi as,

pi=4*hits/shots




Now, I wrote a little script in Python which does just that. I also put in the value of pi to 10 decimal places and wrote out the error produced. Initially, I didn't take a circle of unit radius but that just produced too much error:



import math;
import random;
import sys;


#calculate the value of pi using probability
maxtries = float(sys.argv[1]);


radius = 1.0;
side = 2.0 * radius;
tries = 0.0;
hits = 0.0;
pi = 3.1415926536;


while (tries < maxtries):
tries = tries + 1;
#get a point within the square


rx = random.random()*side;
ry = random.random()*side;


#distance from center of circle
dist = math.sqrt(math.pow(rx - radius, 2.0) + math.pow(ry - radius, 2.0));
if (dist <= radius): hits = hits + 1;


#estimated value of pi
piest = (hits*4/tries);
print "Radius: %d, Tries: %d, Hits: %d, pi: %1.10f, Error:%1.10f\n" % (radius, tries, hits, piest, piest - pi);





My output looks something like:

~$ python pi.py 10.0
Radius: 1, Tries: 10, Hits: 6, pi: 2.4000000000, Error:-0.7415926536

~$ python pi.py 100.0
Radius: 1, Tries: 100, Hits: 83, pi: 3.3200000000, Error:0.1784073464

~$ python pi.py 1000.0
Radius: 1, Tries: 1000, Hits: 790, pi: 3.1600000000, Error:0.0184073464

~$ python pi.py 10000.0
Radius: 1, Tries: 10000, Hits: 7822, pi: 3.1288000000, Error:-0.0127926536

~$ python pi.py 100000.0
Radius: 1, Tries: 100000, Hits: 78591, pi: 3.1436400000, Error:0.0020473464

~$ python pi.py 1000000.0
Radius: 1, Tries: 1000000, Hits: 786368, pi: 3.1454720000, Error:0.0038793464

~$ python pi.py 10000000.0
Radius: 1, Tries: 10000000, Hits: 7853556, pi: 3.1414224000, Error:-0.0001702536

~$ python pi.py 100000000.0
Radius: 1, Tries: 100000000, Hits: 78539696, pi: 3.1415878400, Error:-0.0000048136

~$ python pi.py 1000000000.0
Radius: 1, Tries: 1000000000, Hits: 785382068, pi: 3.1415282720, Error:-0.0000643816

The last observation is the most interesting. It ran for a few hours on my machine but did absolutely nothing to improve the error. In fact, the error actually increased!

6 comments:

Sharad said...

Cool stuff, I did know about this pi approximation. But great to see someone go and write a program.

On your observation, I guess its just one experiment. I could toss a coin a 100 times and get 45/55 heads tails. And I could toss it a 1000 times and get 420/580. Odd to say this for a mathematical concept - luck of the draw :) But hey, I wasn't there eagerly waiting for it to finish.

By the way (sorry, perfectionist speaking here), when you copy/paste-ed the program from your editor you forgot the 'tries' increment.

Cool stuff, thanks for posting. Well worth the read.

Unknown said...

Thanks! fixed it! I am glad that no one is plagiarizing my code. :)

I actually tried the last observation a couple of times before giving up.

I thought this would be a binomial distribution with p = pi/4 with a probability of getting k successes increasing with the value of n (tries). Thus I expected the error to decrease with more observations.

I'll work the math out sometime and let you know what results I get.

Sharad said...

A very basic understanding of Python is needed to plagiarize your code (as the indentation has been lost). Guess what I did :)

It is binomial with p = pi/4. IMHO, the probability of k successes is meaningless in this matter.

(Conjecture follows). As we start running into the limits of floating point arithmetic, the rounding that the computer is performing starts becoming significant. Ipso facto, I am pretty sure it can be proved that increasing our tries beyond a point will not help improve the error.

Unknown said...

Wrote it in Java?

btw, based on your conjecture this would mean that decreasing the rounding that needs to be done should increase the accuracy of the estimation of pi. Instead of multiplying by 4 should the error decrease if we multiplied by 4 billion?

Sharad said...

No, Python is probably the better language. I just meant that I copied and ran it :)

I did think about your Q when I responded earlier - would it help if we worked with larger numbers? Guess an analogous question would be: will a computer provide a more accurate value for, lets say, 1/983745 or for 1000000/983745. Too much pizza inside me, can't think about it right now. My gut-feel is that we would change the significance of the rounding error by the very amount that we would reduce it's magnitude.

Anonymous said...

Who knows where to download XRumer 5.0 Palladium?
Help, please. All recommend this program to effectively advertise on the Internet, this is the best program!