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]\).

Gaussian funciton

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.

Convergence

Source code