Monday, January 5, 2009

Monte Carlo Simulation for Dummies

History: The idea behind Monte-Carlo simulations gained its name and its first major use in 1944, in the research work to develop the first atomic bomb. The scientists working on the Manhattan Project had intractably difficult equations to solve in order to calculate the probability with which a neutron from one fissioning Uranium atom would cause another to fission. The equations were complicated because they had to mirror the complicated geometry of the actual bomb, and the answer had to be right because, if the first test failed, it would be months before there was enough Uranium for another attempt.

They solved the problem with the realization that they could follow the trajectories of individual neutrons, one at a time, using teams of humans implementing the calculation with mechanical calculators. At each step, they could compute the probabilities that a neutron was absorbed, that it escaped from the bomb, or it started another fission reaction. They would pick random numbers, and, with the appropriate probabilities at each step, stop their simulated neutron or start new chains from the fission reaction.

The brilliant insight was that the simulated trajectories would have identical statistical properties to the real neutron trajectories, so that you could compute reliable answers for the important question, which was the probability that a neutron would cause another fission reaction. All you had to do was simulate enough trajectories.

 

I know for most of people that went above there heads. So lets try to look it in a simple way.

Monte Carlo algorithm to compute the value of tex2html_wrap_inline68478from a sequence of random numbers. Consider a square positioned in the x-y plane with its bottom left corner at the origin as shown in Figure. The area of the square is r2, where r is the length of its sides. A quarter circle is inscribed within the square. Its radius is r and its center is at the origin of x-y plane. The area of the quarter circle is π*r2/4.

  
Figure: Illustration of a Monte Carlo method for computing π.

Suppose we select a large number of points at random inside the square. Some fraction of these points will also lie inside the quarter circle. If the selected points are uniformly distributed, we expect the fraction of points in the quarter circle to be

f = (πr2/4) / r2 = π/4

Therefore by measuring f, we can compute π. Program shows how this can be done.

class Program

    {

        static void Main(string[] args)

        {

            int trials = 10000000;

            Console.WriteLine("The value generated after number of {0} trials is {1}",trials, Pi(trials));

            Console.ReadLine();

        }

 

        public static double Pi(int trials)

        {

            Random rnd = new Random();

            int hits = 0;

            for (int i = 0; i <>

            {

                double x = rnd.NextDouble();

                double y = rnd.NextDouble();

                if (x * x + y * y <>

                    ++hits;

            }

            return 4.0 * hits / trials;

        }

    }

Program: Monte Carlo program to compute π.

The Pi method uses the Random.NextDouble() defined to generate (x,y) pairs uniformly distributed on the unit square (r=1). Each point is tested to see if it falls inside the quarter circle. A given point is inside the circle when its distance from the origin, Root(x2+y2) is less than r. In this case since r=1, we simply test whether x2 + y2 <>

How well does Program work? When 1000 trials are conducted, 792 points are found to lie inside the circle. This gives the value of 3.168 for π, which is only 0.8% too large. When 108 trials are conducted, 78535956 points are found to lie inside the circle. In this case, we get π = 3.14143824 which is within 0.005% of the correct value!