Finding the value of Pi

If you're interested in statistics or probability, then you'll probably love the Count Bayesie blog. In one of its recent posts gave a bunch of examples of Monte Carlo simulation. While that sounds complex its really just fancy maths talk for generating random numbers, and using them to test something. If we test something in a bunch of random ways, then the results will be representitive of the thing we're testing.

It's a pretty fundamental idea when it comes to solving difficult maths using a computer, and represents a fundamental shift from the way we used to do things "in the old days". While once upon a time we might have tried to solve a complex problem by converting it into a single equation which we could evaluate to give us an answer, once we have a computer its often quicker/easier to just let it do the work and run millions of simple tests rather than one bigger one.

Lets take a basic example: Imagine a square, entered at the origin which is 2 units on each side. It has an area of 4. Now imagine a circle of radius one. It sits in the middle of the square and just touches each side. It has an area of Pi r^2, or (since r=1) simply pi.

If we pick a random point inside the square it might fall inside the circle or not. What's the probability of that? Well the points are spread over an area of 4, but we only catch those that and in the cenral area of size Pi so we catch Pi/4 of them.

So if we generate n points, then we should catch nPi/4 of them in the circle as "hits". Jiggle that around and we see that Pi=4h/n. What's the point of that? It means you can calculate Pi by doing the experiment and dropping pins randomly onto a piece of paper. However this is actually harder do do than it sounds, and takes ages so this stuff really only took off once we got computers. Here's the code in Sniff:

make hits number
make iterations number
make x number
make y number
when start
.set iterations to 10000
.set hits to 0
.repeat iterations
..set x to (pick random -1000000 to 1000000)/1000000
..set y to (pick random -1000000 to 1000000)/1000000
..if x*x+y*y<1
...change hits by 1

.say [hits/iterations * 4 ]

This prints out 3.1444 which isn't too bad. The great thing is that the more times we go around the loop the better it gets. 10 million iterations gives 3.14141 which is getting pretty close (to 3.14159265...). Trying to go past 10,000,000 iterations runs into a problem - regular computers need special code to handle such big numbers correctly, and Sniff isn't set up for it, so this is about as accurate as we can get with this simple approach. Going the other way, if we want the answer faster 100 iterations gives an answer of 3.08 (bonus question: plot the error against iterations, to see how quickly we get the right answer).

This approach can be applied to all sorts of problems (including lots in computer graphics used for movies), but its generally used for "hard" problems that mathematicians can't solve any other way, so some of the examples can get a bit tricky. However the basic idea is simple enough its worth trying out.