Integrating a function, put simply, allows one to find the area under that function over a certain domain. Expanding this into higher dimensions, the integral of a 3-D function allows one to find the volume under that function. This is all fine and dandy if the function is easily integrable, but what if some function is impossible to integrate?

Should a function be impossible to integrate (also stated as, has no elementary anti-derivatives), then there are numerical solutions, basically ways to guess really well what the true solution is. Monte Carlo Integration is one of these numerical solutions.

1. Find some region that you want to integrate over. Call the region \(R\). This region should have as many dimensions as there are inputs to the function you want to integrate. For example, for \(f(x, y, z, w)\), \(R\) should be in 4 dimensions.

2. Find the volume/area of \(R\). Denote the volume/area as \(V\)

3. Choose many points (the more the merrier!) that are uniformly distributed within \(R\). Let the set of points be \(X\), where the number of points in \(X\) is \(N\).

4. Find the sum of \(f\) evaluated at all points in \(X\), \(\sum{f(x)}, x \in X\). Denote this sum as \(S\).

5. Let \(Q_{N}=SV/N\).

6. As \(N\) increases, \(Q_{N}\) approaches the value of the integral of \(f\) over \(R\).

This method is stochastic, as the values in \(X\) are randomly selected. For very complicated reasons beyond this presentation, deterministic numerical integration methods error increases exponentially with the number of dimensions, whereas the error of Monte Carlo Integration is independent of the number of dimensions of the function. The downside of Monte Carlo Integration is that it is computationally expensive, but it tends to be worth it in high dimension problems.

The standard toy problem to explain Monte Carlo Integration is to find the volume of (read, the integral of) the function \[ f(x, y) = \begin{cases} 1 & \text{if $\sqrt{x^{2}+y^{2}} < 1$},\\ 0 & \text{if $\sqrt{x^{2}+y^{2}} \geq 1$}, \end{cases} \] over \(-1 < x < 1\) and \(-1 < y < 1\). You may have noticed that \(f(x)\) is \(0\) everywhere except within the unit circle, where its value is \(1\). Using the formula for the volume of a cylinder, we can then tell that the volume is \(V = (base)(height) = (\pi r^{2})(1) = \pi\). Therefore, it will be easy to check whether or not the Monte Carlo Integration is working correctly!

MonteCarloIntegrateToyProblem = function(N)
{
  # Step 2
  V = (1 - (-1))*(1 - (-1))
  # Step 3
  XValues = runif(N, min = -1, max = 1)
  YValues = runif(N, min = -1, max = 1)
  # Step 4
  S = 0
  for (n in 1:N)
  {S = S + circle(XValues[n], YValues[n])}
  # Step 5
  Q_N = S*V/N
  return(Q_N)
}

Now, lets look at a more realistic case where one would use Monte Carlo Integration. For this example, we shall use a multi-dimensional function with no elementary antiderivative. Looking at a table of functions with no elementary antiderivative, we shall use \[ f(x, y, z) = \frac{e^{(x^{2})}\sin{(y^{2})}}{1+e^{z}}, \] Integrating over \(-1 < x < 1\), \(-1 < y < 1\), and \(-1 < z < 1\), as this is impossible to integrate and is in several dimensions! Now, lets put it into R:

func = function(x, y, z)
{
  return((exp(x^2)*sin(y^2))/(1+exp(z)))
}

Well there you have it, the integral of \(f(x, y, z)\) is about 1.816793, or in mathematical notation, \[ \int_{-1}^{1}\int_{-1}^{1}\int_{-1}^{1}[\frac{e^{(x^{2})}\sin{(y^{2})}}{1+e^{z}}]\,dx\,dy\,dz\,\approx1.816793 \] Now, some of you may be wondering: how often does one really need to numerically integrate high-dimensional functions? After all, we do live in a mere three dimensions. However, in many cases, the inputs of a function may be any slew of variables that are unrelated to the space we live in. Therefore, it is not uncommon to have to numerically integrate a function with far more than three inputs.

To have an intuitive look at the function, and to show off a cool plot, lets view \(f(x, y, z)\) in 3-D:

In conclusion, Monte Carlo Integration is a very useful numerical integration technique that particularly shines in high-dimensional problems. This is because the method is inherently stochastic. This requires the generation of random numbers, which is computationally expesive, and thus the method tends to be outshone by other methods in low dimensions.