# Integration using the Monte Carlo method

In a previous article we used the Monte Carlo method to estimate the value of $$\pi$$. Here, we will use the Monte Carlo method in one of its most common applications: numeric integration. The idea behind Monte Carlo integration is the same we applied to estimate $$\pi$$. The estimation of an integral in a region with volumne $$V$$ is done by estimating the fraction of random points that fall below the function multiplied by $$V$$.

# When is Monte Carlo integration used

The convergence of Monte Carlo integration is of order $$O(n^{1/2})$$ and is independent of the dimensionality. Hence, Monte Carlo integration usually beats numerical intergration for high dimensional integration since numerical integration (quadrature) converges as $$O(nd)$$m where $$d$$ is the number of dimensions.

Even for low dimensional problems, Monte Carlo integration may have an advantage when the volume to be integrated is concentrated in a very small region and we can use information from the distribution to draw samples more often in the region of importance.

# The case study

Typically the method would be applied to calculate an integral of a function without a closed-form analytical solution. However, in order to have an analytical solution as benchmark we will use a function with a known definite integral. Consider the Gaussian function:

$$f(x) = a e^{-\frac{(x-b)^2}{2c^2}}$$

where $$a = 0.4$$, $$b = 0$$, and $$c = 1$$. Let's calculate the definite integral between $$[-4 4]$$.

## Analytical solution

We can use the Python sympy package to calculate the analytical solution:

from sympy import symbols, integrate, exp

x = symbols('x')
func = 0.4 * exp(-(x-0)**2 / 2*1**2)
expr = integrate(func, (x, -4, 4))
integral_f_gauss_analyt = expr.evalf()


The analytical solution is: $$1.00258779942818$$.

## Deterministic numeric integration

As another comparisson, lets calculate the integral numerically using the trapezoid rule.

import numpy as np

# Create the Gaussian function
def f_gauss(x, a=0.4, b=0, c=1):
"""Gaussian function
Has closed-form analytical solutions for the integral
x : independent variable
a, b, c : function parameters
"""
return a * np.exp(-(x-b)**2 / 2*c**2)

# Define domain
INT_RANGE = (-4, 4)

# Define discretization steps
N = 1000

# Calculate the integral using the trapezoid rule
x = np.linspace(INT_RANGE[0], INT_RANGE[1], N)
y = f_gauss(x)
integral_f_gauss_trapz = np.trapz(y, x=x)


The result of the trapezoid rule is: $$1.0025877936915661$$.

## Monte Carlo integration

Finally, let's calculate the same integral using the Monte Carlo method. This is a non-deterministic approach, every run will result in a slightly diferent estimation.

# Define the number of random points
N = 1000000

# Create the random points
x_points = [np.random.uniform(-4, 4) for _ in range(N)]
# Calculate the integrand
f_x = [f_gauss(x) for x in x_points]
# Calculate the volume of the domain
V = INT_RANGE[1] - INT_RANGE[0]
# Calculate the integral
integral_f_gauss_montecarlo = sum(f_x) * V / N


The result of the Monte Carlo method is: $$1.0022406995836632$$.

# Monte Carlo method convergence

Here is how the Monte Carlo method converges to the final estimate using $$n$$ random points.

## Source code

• Jupyter Notebook available here