Motivating Monte Carlo Methods

In this section, I’m going to do my best to convince you that Monte Carlo methods can greatly improve the quality of your life. But first, a warm-up exercise to get your mathematical juices flowing.

Consider the following function:

g(x) = \frac{6 \pi x^4 \exp(-x^2) (\log(x) + 3^x)}{(x+6)^{1/2}}.

Now, compute the following definite integral of this function:

\int\limits_{0}^{1} \frac{6 \pi x^4 \exp(-x^2) (\log(x) + 3^x)}{(x+6)^{1/2}} dx .

Go. Do it. Open a calculus textbook if you need to. Google “solving complicated definite integrals”. Don’t be lazy. I promise it’ll be worth the effort.

Did you figure it out? Me neither.

I have no idea how in the hell to solve this problem, nor am I particularly interested in diving into my old calculus textbooks to try to do so.

What if I were to tell you that we could solve this problem with just 3 simple lines of code? Too good to be true? Nope – I’m telling the truth here. What method will I be using? You guessed it – a simple Monte Carlo method.

Monte Carlo to the rescue

As I mentioned in the previous section, Monte Carlo methods leverage randomness in order to solve problems. The way we incorporate randomness towards solving problems is by performing random sampling.

Let’s tackle the integral posed at the beginning of this section. Let

X_1, ..., X_n \overset{iid} \sim U(0,1), 

denote a sample of size n taken from the uniform distribution on (0,1) (in statistics and probability, the notation ‘iid’ stands for independent and identically distributed. In other words, each sample is drawn independently of the other samples, and all samples are drawn from the same distribution). In layman’s terms, just pick n random numbers from 0 to 1.

Next, I want to emphasize a very important point. Even though the integral

g(x) = \frac{6 \pi x^4 \exp(-x^2) (\log(x) + 3^x)}{(x+6)^{1/2}}.

looks very complicated, it is trivial to plug values of x into this equation. For example, if x = 0.22, then

g(0.22) = \frac{6 \pi 0.22^4 \exp(-0.22^2) (\log(0.22) + 3^{0.22})}{(0.22+6)^{1/2}}.

Namely, we can easily plug in the n samples we drew from the uniform distribution into this equation to get values

g(X_1), g(X_2), ..., g(X_n).

Before moving ahead, let’s recall what the law of large numbers tells us. Given a large enough sample  X_1, ..., X_n , the sample mean \bar{X} will be a good approximation to the population mean E[X] ; in other words,

\frac{1}{n} \sum_{i=1}^{n} X_i \approx  E[X].

Similarly, the law of large numbers implies that

\frac{1}{n} \sum_{i=1}^{n} g(X_i) \approx  E[g(X)],

Don’t let the g() terms throw you off. We are simply stating that the mean of a bunch of numbers (left side of the approximation symbol) will be a good approximation to the population mean of those numbers (the right side of the approximation symbol).

By the definition of expectation, we immediately see that

E[g(X)] = \int_{0}^{1}g(x)1dx 
is exactly the integral we are trying to approximate (I put the 1 inside of the integral to remind you that the pdf of a uniform(0,1) random variable is 1)! Therefore we have found a simple scheme to approximate this integral!

Monte Carlo method to approximate the nasty integral

  1. Generate X_1, ..., X_n \overset{iid} \sim U(0,1) with large n
  2. Compute g(X_i) for each X_{i}
  3. Compute the mean of the g(X_i) ‘s. The mean will give us an approximation to the integral we seek.

As I promised, the corresponding code is very simple:

nasty_function <- function(x){6*pi * (x^4)*exp(-x^2)*(log(x) + 3^x)/((x+6)^(0.5))}
u <- runif(n = 10000, min = 0, max = 1)

This gives me the value 1.56836. We can use R to find the exact (well, up to a very small error term) solution to this problem:

integrate(nasty_function, lower = 0, upper = 1)

The exact answer, it turns out, is about 1.589. Our simple Monte Carlo method did a pretty good job! Of course, we can increase n to get better and better approximations.

What happens if we repeat this process many times – that is, we take a random sample of size n, compute its mean, take another random sample of size n, compute its mean, and on and on. Could you predict what will happen?

Let's code it up and see for ourselves.

results <- c()
for(i in 1:1000){
u <- runif(n, min = 0, max = 1)
results[i] <- mean(nasty_function(u))

hist(results, prob = TRUE, main = "Histogram of 1000 approximations of nasty function")


We see the familiar gaussian distribution centered around the true value of the integral. Of course, this shouldn’t surprise you – it is an instance of the Central Limit Theorem in action.

One more example



I hope I got your appetite wet for Monte Carlo!