Home Course Concepts About

This notebook is an element of the risk-engineering.org courseware. It can be distributed under the terms of the Creative Commons Attribution-ShareAlike licence.

This notebook contains an introduction to the use of NumPy for numerical integration using the Monte Carlo method.

Monte Carlo methods for numerical integration

This notebook contains an illustration of the use of Monte Carlo methods for numerical integration. See the associated course materials for an introduction to the use of stochastic simulation methods and to download this content as a Jupyter/Python notebook. Numerical integration is often used to evaluate risk measures in the finance industry.

First integral

Task: resolve the integral

$\int_1^5 x^2 dx$

Analytical method

We start by resolving this integral using the standard analytical method, assisted by the SymPy symbolic mathematics library.

In [1]:
import sympy
from sympy.interactive import printing

x = sympy.Symbol("x")
sympy.integrate(x**2, (x, 1, 5))
$\displaystyle \frac{124}{3}$
In [2]:
# note that _ references the value of the previous cell
$\displaystyle 41.3333333333333$

Numerical method

The expression we wish to integrate is very simple here, so calculating its integral analytically is easy (even without resorting to Python’s symbolic mathematics package!). In many cases, however, analytical approaches to integration are not feasible:

  • the expression we wish to integrate is very complicated, possibly without a closed analytical integral
  • it is only known in “black box” form (for instance a software module): we can evaluate it at various points but don’t know the corresponding equation

In such situations, stochastic simulation (“Monte Carlo”) methods allow us to generate an approximation of the integral, simply by evaluating the expression a large number of times at randomly selected points in the input space and counting the proportion that are less than the integrand at that point. The larger the number of simulations we run, the better the approximation (and note that computers are very very fast today, so on simple problems the simulation will run in a blink of an eye!).

In [3]:
import numpy
In [4]:
N = 100_000
accum = 0
for i in range(N):
    x = numpy.random.uniform(1, 5)
    accum += x**2
measure = 5 - 1
measure * accum/float(N)
$\displaystyle 41.2241963497441$

Second integral

Task: resolve the integral $\int_3^2 x^2 + 4 x sin(x)$.

Analytical solution

In [5]:
x = sympy.Symbol("x")
float(sympy.integrate(x**2 + 4 * x * sympy.sin(x), (x, 2, 3)))
$\displaystyle 11.8113589250983$

Stochastic solution

In [6]:
N = 100_000
accum = 0 
for i in range(N):
    x = numpy.random.uniform(2, 3)
    accum += x**2 + 4*x*numpy.sin(x)
measure = 3 - 2
measure * accum/float(N)
$\displaystyle 11.8102144834769$

Exercise: now undertake the same comparison of analytical and stochastic simulation methods to evaluate the integral

$$\int_1^3 e^{x^2}$$

(Hint: the result should be approximately 1464.2)

A 2D integral

Now we move to a 2D integral, $\int\int cos(x^4) + 3y^2 dx dy$ over $x ∈ [4,6]$ and $y ∈ [0, 1]$.

Monte Carlo integration methods become increasingly attractive when compared to other numerical integration methods (such as quadrature) as the number of dimensions increases.

Let’s examine the analytical solution, again thanks to sympy.

In [7]:
x = sympy.Symbol("x")
y = sympy.Symbol("y")
float(sympy.integrate(sympy.cos(x**4) + 3 * y**2, (x, 4, 6), (y, 0, 1)))
$\displaystyle 2.00505508674967$

Now compare with Monte Carlo estimation:

In [8]:
N = 100_000
accum = 0
for i in range(N):
    x = numpy.random.uniform(4, 6)
    y = numpy.random.uniform(0, 1)
    accum += numpy.cos(x**4) + 3 * y * y 
measure = 2 * 1
measure * accum/float(N)
$\displaystyle 1.9946300281134$
In [ ]: