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_{i1}$ and $x_i^{\ast} \in [x_{i1}, 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_{i1}$ and $x_i^{\ast} \in [x_{i1}, 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.
1 

1  $ ./pi 
The $\pi$ computed approximately using Monte Carlo integration is $3.14175$.
References
Monte Carlo Integration