Monte-Carlo simulation

The expression "Monte-Carlo simulation" refers to a large number of techniques that have the same approach to solving complex problems :

• The problems are usually completely deterministic, and could theoretically be solved by exact mathematical calculations or at least by numerical approximation methods.
• But in practice, the mathematical solutions are intractable, and/or numerical approximation methods take an impractically long time to reach a result.

This kind of situation is very common in all fields of  mathematics, physics, chemistry, biology, economy, finance, sociology etc...The problems then encountered involve many variables, functions or differential equations.

It turns out that it is sometimes possible to solve such a problem numerically by introducing a random sampling scheme of the domain of the variables.

A method resorting to such an approach will be generically called a "Monte-Carlo simulation".

1-dimensional integration

We here address one of the simplest applications of Monte-Carlo simulation.

How do we calculate the integral between a and b of the function represented in the illustration below ?

If no primitive function of f(x) is known, the integral cannot be calculated analytically. But if the function f(x) can be easily calculated at any point of the integration region, we can obtain a good approximation of the integral by resorting to numerical methods.

There are many such methods. The simplest one (and also the least effective !) is to divide up the range of x into a series of N adjacent rectangles (see lower image of the above illustration). The height of each rectangle is adjusted to the value of f(x) for x taken as the mid-point of the base of the rectangle. The (algebraic) sum of the areas of the rectangles is taken as an approximation of the (algebraic) area under the curve representing f(x), that is, its integral I :

If f(x) is well-behaved, and if the rectangles are narrow enough, then the calculated area will be a close approximation of the value of the integral.

-----

Note that this approximation of the integral is, within a geometric normalization factor, the sum of the values of f(x) taken at evenly spaced points in the integration region.

Multi-dimensional integration

Now what if we try to use this technique for functions of more than just one variable ? We then run into two new difficulties :

1. The integration region is not defined by just a pair of numbers as before. Instead, it is a closed (hyper-)surface whose shape may be very complicated even for simple problems (see interactive animation) and impossible to describe analytically.
    
2. The second difficulty is even more serious. Let's get back to the 1-dimensional integral, and suppose we decided to use 100 rectangles (a small number indeed). We then turn to a 100-dimensional problem (not so high by today's standards), and decide to keep the same resolution as before on each of the axes. We now need to define 100¹⁰⁰  hyper-rectangles, a number out of reach of even the fastest computers.
   This ubiquitous problem is encountered in any local technique (either deterministic or probabilistic), and is called the "curse of dimensionality".
   Note that the number of rectangles is chosen mostly by considering the rate of change of the function, and that we therefore cannot arbitrarily decide to reduce the number of rectangles down to a manageable level, lest we lose so much of the detailed description of the function that further calculations will be grossly off target.

Monte-Carlo integration

Now comes an idea that is both simple and effective.

We noted that an integral is just the sum of the values of the function taken at evenly spaced points of the integration region. Now replace "evenly" by "randomly uniformly distributed". You have now your first Monte-Carlo method, so simple that it is often affectionately called the "crude" or "brute force" Monte-Carlo integration method.

So, in its simplest version, a Monte-Carlo calculation of the integral :

would go as follows :

• Use a n-dimensional uniform distribution to draw N points p_i  (i = 1, 2, ..., N)  inside the integration region (whose volume is V).
•  Add all the f(p_i), multiply this sum by V/N.
• The result is an estimate of the value of the integral I.
   

Advantages of the Monte-Carlo integration method

Crude as it is, the simple Monte-Carlo integration method makes the two aforementioned problems less critical :

1. We do not need to know the mathematical form of the boundaries of the integration region anymore, but only its volume, which can be estimated if we have a simple rule that tells us whether a given point is "inside" or "outside" this integration region, a much simpler problem (see interactive animation).
    
2. But the main interest of the Monte-Carlo integration method lies in its behavior in high-dimension spaces. This is a difficult question, but whose answer can be expressed in simple words : take any deterministic integration method that relies on calculating the values of the function to be integrated on the N vertices of a regular multi-dimensional grid (again, there are many such methods). Then, as we consider higher and higher dimensions, there will ultimately come a dimension d when the Monte-Carlo integration method will be more efficient than the considered deterministic method. What this means is that, for a given value of N, the estimate of the integral produced by the MC method will be, more often than not, closer to the true value of the integral than the value produced by the deterministic method.

-----

What is the price to pay for these advantages ?

Because of the random nature of the sampling of f(x), the number obtained at the end of a MC run is the realization of a random variable. In other words, if you run several MC calculations of a given integral, you'll never get twice the same value. These values follow a probability distribution with a definite variance. So whereas a deterministic integration method provides an approximation of the true value of an integral, the MC integration method provide an estimate of this value. The difference, although subtle, is important : the error of an approximation can often be bounded upwards, whereas there is usually no upper limit to the error of an estimation. But by making the number of random points large enough, the probability that this error exceeds a given value e can be made arbitrarily small (this is just a consequence of the Weak Law of Large Numbers). Of course, a consequence is that Monte-Carlo simulations runs are always very long.

Interactive Animation

 We illustrate here the concept of Monte-Carlo integration on two problems.

• Estimation of the area of an irregular and user-defined region in the [0, 1]² square. Although the problem appears to be very simple, no analytical solution exists, and even this elementary simulation provides a useful practical way of estimating this area.
• Calculation of p. This is achieved by estimating the area of a disc inside the [0, 1]² square. Needless to say, this approach is not to be considered as a practical way of calculating p, but rather as a pedagogical illustration of the concept of "Monte-Carlo integration method".
  

____________________________________________________________

Related readings