Pi day is March 14th - get it? (3.14) I am a big fan of Pi. Here is my first post to celebrate the awesomeness of Pi (I know this is early, but I was too excited to wait).
How can you determine Pi?
Oh sure, tons of high schools do the classic experiment. Measure the circumference and diameter of as many round things as possible. Plot diameter vs. circumference. The slope will be Pi. Really, this is a great lab to do for all sorts of ages. The key thing is that students can see what Pi really means. I am not going to talk about this lab, I am going to do some thing cooler.
What if I had a 1 meter by 1 meter square taped out in the grass near the physics building and shot ping pong balls off the roof at this square? Ping pong balls are very difficult to aim, so you could easily imagine that they would fall in a random pattern in the square. I don't really want to set this up, so I am going to do it with a program instead. Here is a program that will generate random points in a 1 x 1 square.
This is what the output would look like:
Now, here is the key. What if I calculate the distance of each of these point from the origin? This would simply be:
Doing this, I can separate all the points into those that are more than 1 unit away from the origin and those that are less than 1 unit away. In this plot, (of 1000 points) the red dots are more than one unit away and the blue are 1 unit or less.
Note that this plot does not have exactly the same horizontal and vertical scale - no idea why it came out like that. However, this is enough dots that maybe you can see a pattern. The blue dots are filling up 1/4th of a circle. The full area of this circle would be:
Since I am a physicist, I have trouble leaving the area as unitless. What if I look at the ratio of blue dots to total dots? This should be the same as the ratio of the area of the quarter circle to the whole square, or:
Doing this for my last random number run, I get pi = 3.04. What happens as I increase the number of random dots? This a plot of the estimation of Pi for numbers of dots starting at 1000 up to 100,000 dots. I did this 5 times because each time is different.
So, you can see that as the number of random points increases, the estimated value of Pi is a little less spread out and closer to 3.14-something.
P.S. Thanks to Dave for this idea. You know who you are.
I also made a Scratch program version of this calculation. (Scratch is a graphical programming language made at MIT - you know, for kids).
IIRC this could be thought of as a Monte Carlo method for evaluating Pi.
Yes - good point, I guess I should have mentioned that.
Also reminiscent of one of the first documented usages of a Monte Carlo method (from http://en.wikipedia.org/wiki/Asaph_Hall):
On June 5, 1872 Hall submitted an article entitled "On an Experimental Determination of Pi" to the journal Messenger of Mathematics. The article appeared in the 1873 edition of the journal, volume 2, pages 113-114. In this article Hall reported the results of an experiment in random sampling that Hall had persuaded his friend, Captain O.C. Fox, to perform when Fox was recuperating from a wound received at the Second Battle of Bull Run. The experiment involved repetitively throwing at random a fine steel wire onto a plane wooden surface ruled with equidistant parallel lines. Pi was computed as 2ml/an where m is the number of trials, l is the length of the steel wire, a is the distance between parallel lines, and n was the number of intersections. This paper is a very early documented use of random sampling (which Nicholas Metropolis would name the Monte Carlo method during the Manhattan Project of World War II) in scientific inquiry.
Dude! Someone else that uses matplotlib! Awesome!
But Hall cheated. He stopped when the result was very close to the known value of pi. One more trial and the result would have gone a long way from pi in one direction or another. While the procedure is valid the number of trials to be confident of getting a good result makes it infeasable.
I decided to add a little more pie to your pi estimates. There is a screen shot and other details (including the R code) here: How to make pi... using R.
Note that this plot does not have exactly the same horizontal and vertical scale - no idea why it came out like that.
In order to square the axes, you need to insert the following code after you plot the points:
ax = gca()
ax.set_aspect(1.0 / ax.get_data_ratio())
Thanks! That worked.
Hey Rhett. Super cool post.
What about the rest of the code in Python. Can we see it?
Here is the code (sloppy code warning).
Let's have some fun with Pi. ;)