Monte Carlo Integration

Introduction

It is very common that we want to compute an integral for some $f(x)$ defined in $\Omega$.

$$
I = \int_{\Omega}^{} f(x) dx
$$

However, in some circumstances, the analytical result of the expected value might not be derived easily if $f(x)$ is very complex.

In this blog post, I would like to discuss how to compute an integral approximately using Monte Carlo integration.

Prerequisites

Law of Large Numbers

Law of large numbers states that for random variable $X$ with $\mathbb{E}[X] = \mu$ and $\mathbb{Var}[X] = \sigma^2$ and we also have a sequence of independent and identically distributed (i.i.d.) random variables $X_1, X_2, \cdots, X_n$, the sample average random variable

$$
\overline{X}_n = \frac{1}{n} \sum_{i=1}^{n} X_i
$$

converges to $\mu$

$$
\overline{X}_{n} \rightarrow \mu, \text{ when } n \rightarrow \infty
$$

This means that for any random variable, if the number of random samples is sufficiently large, the average of random samples from the i.i.d random variables will be sure to be the population mean $\mu$.

This is quite similar to the central limit theorem, except that the law of large numbers did not tell the distribution of sample average random variable will become a Gaussian distribution as $n \rightarrow \infty$.

Therefore, to estimate $\mathbb{E}[X]$, we prepare $n$ random samples $x_1, x_2, \cdots, x_n$ from i.i.d random variables $X_1, X_2, \cdots, X_n$, and $n$ is sufficiently large.

$$
\mathbb{E}[X] \approx \frac{1}{n} \sum_{i=1}^{n} x_i
$$

Riemann Integral and Riemann Sum

Riemann Integral

Riemann integral is actually the definition of the definite integral.

$$
\int_{\Omega}^{} f(x) dx = \lim_{||\Delta x|| \rightarrow 0}^{} \sum_{i=1}^{n} f(x_i^{\ast}) \Delta x_i
$$

where $\Delta x_i = x_i - x_{i-1}$ and $x_i^{\ast} \in [x_{i-1}, x_i]$.

Riemann Sum

Following Riemann integral, Riemann sum provides a way to compute the definite integral approximately.

$$
\int_{\Omega}^{} f(x) dx \approx \sum_{i=1}^{n} f(x_i^{\ast}) \Delta x_i
$$

where $\Delta x_i = x_i - x_{i-1}$ and $x_i^{\ast} \in [x_{i-1}, x_i]$.

As a special case of Riemann sum, suppose the volume of $\Omega$ is $V$, where

$$
V = \int_{\Omega}^{} dx
$$

We evenly divided the area $\Omega$ into $n$ pieces to get all the $x_i$, then $\Delta x_i = \frac{V}{n}$. We further set $x_i^{\ast} = x_i$. We have the Riemann sum

$$
\begin{align}
\int_{\Omega}^{} f(x) dx &\approx \sum_{i=1}^{n} f(x_i) \frac{V}{n} \\
&= \frac{V}{n} \sum_{i=1}^{n} f(x_i) \\
\end{align}
$$

The advantage of Monte Carlo integration over Riemann sum is that for some $\Omega$ it is easier to do random sampling in $\Omega$ than evenly diving $\Omega$.

Monte Carlo Integration

Monte Carlo integration is somewhat similar to Riemann sum. However, instead of evenly dividing the area $\Omega$ into $n$ pieces to get all the $x_i$, we randomly sample $x_i$ from $\Omega$. The expression of the formula remains almost the same as the special case of Riemann sum we discussed above though.

$$
\begin{align}
\int_{\Omega}^{} f(x) dx &\approx \frac{V}{n} \sum_{i=1}^{n} f(x_i) \\
\end{align}
$$

where

$$
V = \int_{\Omega}^{} dx
$$

and $x_i$ is a random sample of $X_i$ which is an independent and identically distributed (i.i.d.) random variables defined in $\Omega$ that follows uniform distribution.

Proof

By definition, we have the expected value for an arbitrary function.

$$
\begin{align}
\mathbb{E}[f(X)] &= \int_{\Omega}^{} p(x) f(x) dx \\
\end{align}
$$

$$
V = \int_{\Omega}^{} dx
$$

We set $p(x) = \frac{1}{V}$

$$
\begin{align}
\mathbb{E}[f(X)] &= \int_{\Omega}^{} p(x) f(x) dx \\
&= \int_{\Omega}^{} \frac{1}{V} f(x) dx \\
&= \frac{1}{V} \int_{\Omega}^{} f(x) dx \\
\end{align}
$$

Therefore,

$$
\int_{\Omega}^{} f(x) dx = V \mathbb{E}[f(X)]
$$

According to the law of large numbers, since $f(X)$ is a random variable, $\mathbb{E}[f(X)]$ could be computed approximately by sampling. That is to say, we have $n$ i.i.d random variables, $f(X_1), f(X_2), \cdots, f(X_n)$, and we randomly sample $f(x_1), f(x_2), \cdots, f(x_n)$ for $f(X_1), f(X_2), \cdots, f(X_n)$, respectively. To do such sampling for $f(X_i)$, we sample $x_i$ for $X_i$ uniformly from $\Omega$, following by computing $f(x_i)$.

$$
\begin{align}
\int_{\Omega}^{} f(x) dx &= V \mathbb{E}[f(X)] \\
&\approx \frac{V}{n} \sum_{i=1}^{n} f(x_i) \\
\end{align}
$$

This concludes the proof.

Examples

Here we would take a look at the applications of Monte Carlo integration.

Compute $\pi$ Approximately

To be honest, I have totally forgotten how $\pi$ is formally computed. However, with some prior knowledge about $\pi$ and circle, we could compute $\pi$ approximately.

We define a circle whose center is located at $(0,0)$ with a radius of $r = 1$. Formally, $x^2 + y^2 = 1$. With our prior knowledge, we know the area of this circle is

$$
\begin{align}
S &= \pi r^2 \\
&= \pi
\end{align}
$$

Therefore, to compute $\pi$ approximately, we just have to compute the area of the circle approximately.

We define $\Omega = [-1, 1] \times [-1, 1]$ and a function $f(x, y)$ defined in $\Omega$.

$$
\begin{align}
f(x, y) &=
\begin{cases}
1 & \text{if $x^2 + y^2 \leq 1$}\\
0 & \text{else}\\
\end{cases}
\end{align}
$$

$$
\begin{align}
I &= \int_{\Omega}^{} f(x, y) dx dy \\
&= S \times 1 \\
&= S \\
&= \pi \\
\end{align}
$$

Note that the integral $I$ is actually the volume of a cylinder which equals to $\pi$, instead of the area of the circle.

Therefore,

$$
\pi = \int_{\Omega}^{} f(x, y) dx dy
$$

To compute the integral approximately, we could either use Riemann sum or Monte Carlo integration. Here we use Monte Carlo integration to compute the integral.

$$
\begin{align}
V &= \int_{\Omega}^{} dx dy \\
&= 4
\end{align}
$$

Therefore,

$$
\begin{align}
\int_{\Omega}^{} f(x) dx &\approx \frac{V}{n} \sum_{i=1}^{n} f(x_i) \\
&= \frac{4}{n} \sum_{i=1}^{n} f(x_i) \\
\end{align}
$$

We could then write a compute program to help us finish the following sampling tasks.

pi.cpp

#include <random>
#include <iostream>

double f(double x, double y)
{
    double val = x * x + y * y;
    if (val <= 1.0)
    {
        return 1.0;
    }
    else
    {
        return 0.0;
    }
}

int main()
{
    // This will be used to obtain a seed for the random number engine
    std::random_device rd{};
    // Standard mersenne_twister_engine seeded with rd()
    // std::mt19937 gen{rd()};
    // Fixed random seed for reproducibility
    std::mt19937 gen{0};
    std::uniform_real_distribution<> disX{-1.0, 1.0};
    std::uniform_real_distribution<> disY{-1.0, 1.0};
    const int n{10000000};
    double x{0};
    double y{0};
    double sum{0};
    double pi{0};
    double V{4};
    for (int i = 0; i < n; i ++)
    {
        x = disX(gen);
        y = disY(gen);
        sum += f(x, y);
    }
    pi = V / n * sum;

    std::cout << "Pi: " << pi << std::endl;
}

$ ./pi 
Pi: 3.14175

The $\pi$ computed approximately using Monte Carlo integration is $3.14175$.

References

• Monte Carlo Integration