Hands-on Tutorials
Monte Carlo Integration in Python over Univariate and Multivariate Functions
How to use Monte Carlo methods to approximate integration of complex functions
Monte Carlo integration is a basic Monte Carlo method for numerically estimating the integration of a function f(x). We will discuss here the theory along with examples in Python.
Theory
Suppose we want to solve the integration of f(x) over a domain D.
In the case of a univariate function (i.e. with one variable), the domain is simply one-dimensional and the integration is from a to b.
We can rearrange some terms and express the above equation as,
In other words, the integration is equivalent to finding the expected value of g(x), where we have defined g(x)=f(x)/p(x) over a domain. We can approximate this by sampling x from the domain for N times, which means,
If we are sampling from a uniform distribution, and let’s say the function is univariate, this means the probability of drawing any given x is simply p(x)=1/(b-a). If we substitute this into the above approximation expression, we see that,
This is effectively calculating the mean value of f(x) over the interval a to b and multiplying by the length of the interval. In other words, we are finding the area of a rectangle with width = interval width and height = expected value of f(x).
This also works with any dimensions. In case of a univariate function, the domain is simply a line (i.e. b-a). For a bivariate function, the domain is the area. Generally,
This means we can approximate the integration by multiplying the domain volume by the expected value of the function over the domain.
Example | Univariate
As an example, in Python we can perform the following to approximate the integration of f(x)=x² from -2 to 2.
We get the following results,
where the Monte Carlo approximation is very close to the analytical solution. Visually, the integration of f(x)=x² from -2 to 2 is shown below in blue. The approximation is the rectangle highlighted in red.
Example | Multivariate
We can also perform integration for multivariate functions. The procedure is the same as before. However, instead of sampling over a line (from a to b), we now need to sample over a higher-dimensional domain. For simplicity, we will illustrate the integration of a multivariate function over a domain with the same a and b for each variable. This means in a function with two variables (x1 and x2), the domain is square-shaped; and for a function with three variables, cube shaped.
The results
Example | Multivariate integrated over other domains
The domain over which the integration is performed can be more complicated and difficult to sample from and to calculate its volume. We can for example integrate our bivariate function over a circular domain instead of a square domain. The idea is nonetheless the same — with uniform sampling, we wish to sample over the domain, and approximate the integration via the product of the domain volume and expected value of the function over the domain.
Let’s use the same bivariate function f(x)=10-x1²-x2² and integrate over a unit circle. Uniform sampling over exactly the unit circle is harder than just sampling over a square region (that covers the unit circle). From this we can 1) calculate the area of the domain as the product of the area of the sampled square by the proportion of sampled points inside the domain, 2) the expectation as the mean f(x) of the sampled points inside the domain. Below is a visualization of the sampling (sampled points inside the unit circle shown in green),
In Python, this looks like,
With results as,
As a last example, we can also integrate over multivariate probability distributions. Let’s integrate the following multivariate normal distribution over the unit circle,
where
We define these in Python and perform our Monte Carlo integration,
As confirmation, we use the pmvnEll function that can integrate multivariate normal distributions over ellipsoids (circles included).
library(shotGroups)
pmvnEll(r=1, sigma=rbind(c(1,0.8,0.5), c(0.8,1,0.8), c(0.5,0.8,1)),
mu=c(0,0.5,1), e=diag(3), x0=c(0,0,0))With the results closely matching each other,
Additional comments
Obviously the examples provided are simple and some have analytical solutions and/or Python/R packages for specific cases. But they are useful to get a grasp of the mechanics behind Monte Carlo integration. Clearly the Monte Carlo method described is readily generalizable to more complicated functions with no closed form solutions. In addition, there are many more optimized ways to perform sampling (e.g. stratified sampling, importance sampling, etc) and readers are encouraged to read more into those topics if interested.
Additional resources
https://towardsdatascience.com/monte-carlo-integration-in-python-a71a209d277e
Originally published at https://boyangzhao.github.io on January 13, 2021.
