The Monte Carlo Method
Written by Mike James   
Article Index
The Monte Carlo Method
Simulations
Finding Pi
Working out the matrix product

As easy as pie

pie1

Hitting the target – at random

 

OK, finding the area might not impress you but what about calculating one of the fundamental mathematical quantities using randomness. Pi is a well known target for programmers trying to make their name by computing it to a million more decimal places but what about making it a real target for a change.

If you set up a circle with radius one then its area is exactly Pi and so finding the area using random numbers gives you an estimate of pie – I’m tempted to say as easy as…

The idea of finding Pi using random numbers was invented in the eighteenth century by a French naturalist called Georges-Louis Leclerc, Comte de Buffon. He asked what was the probability of a needle thrown at random on a floor made up of planks with a width twice its length landing on a crack between planks. This very linear question has the answer 1/Pi and so you can count the number of times it lands on a crack and use this to estimate 1/Pi and by simple arithmetic Pi itself.

It was actually tried out in 1901 and yielded 355/113 for the value of Pi which isn’t bad considering all they did was toss a needle onto the floor and count the number of times it hit the crack!

Let's implement the algorithm in JavaScript. To make things easier we can place a circle with a radius 0.5 in the middle of the unit square. This circle has area 1/4 Pi and so counting the number of hits as a ratio of the total will estimate 1/4 Pi. 

The program is fairly easy. First we set up a loop that will repeat the random "shot" at the target:

var num = 1000;
var total = 0;
var hit = 0;
for (var i = 1; i <= num; i++) {

Next we generate two random numbers x and y and work out how for this point is from the center of the circle at 0.5,0.5:

var x = Math.random();
var y = Math.random();
var r = Math.sqrt(Math.pow(x - 0.5, 2)
          + Math.pow(y - 0.5, 2));

If the point is closer than 0.5 then we have a hit:

  total++;
  if (r < 0.5) hit++;
 Text1.value = hit / total * 4;
}

If you run the program using Firefox for a range of values you get something like:

num     Pi
10      3.2
100     3.16
1000    3.1
10000   3.1512
100000  3.14
1000000 3.143372

As you can see convergence isn't fast but it generating random numbers is fairly cheap.

The complete program as an HTML page is:

<!DOCTYPE html>
<html lang="en">
 <head>
  <meta charset="utf-8" />
  <title>Random Pi</title>
 </head>
 <body>
  <input type="text" id="Text1"/>
  <script>
   var num = 1000000;
   var total = 0;
   var hit = 0;
   for (var i = 1; i <= num; i++) {
    var x = Math.random();
    var y = Math.random();
    var r = Math.sqrt(Math.pow(x - 0.5, 2)
        + Math.pow(y - 0.5, 2));
    total++;
    if (r < 0.5) hit++;
    Text1.value = hit / total * 4;
   }
  </script>
 </body>
</html>
 

More Than Just Areas

You might be convinced that you can work out areas, volumes and even Pi by random numbers but where next?

The answer is that you can estimate the results of just about any numerical computation using randomness in much the same way. The only problem is that explaining how it works would get us ever deeper into mathematics so a simple example will have to do.

You can use the Monte Carlo method to solve linear equations like Ax=b where b is a known vector and A is a known matrix.

Usually this problem is solved by inverting the matrix or a similar numerical method but when this is large finding the inverse is an expensive problem. Again randomness comes to the rescue and it is possible to estimate the x vector in much the same way as the needle dropping estimated Pi.

The actual steps to get to the solution are complicated but what about just working out the matrix product y=Ab, where A and b are known.

This isn’t such a big problem but it is a step on the way to solving Ax=b and at first it doesn't  seem to have anything at all to do with random numbers.

See if you can work out how to do it first.



 
 

   
RSS feed of all content
I Programmer - full contents
Copyright © 2014 i-programmer.info. All Rights Reserved.
Joomla! is Free Software released under the GNU/GPL License.