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. The associated lecture slides on risk-engineering.org provide an introduction to the use of stochastic simulation methods.

$\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
printing.init_printing(use_latex='mathjax')

x = sympy.Symbol('x')
integral = sympy.integrate(x**2)
integral.subs(x, 5) - integral.subs(x, 1)

Out[1]:
$$\frac{124}{3}$$
In [2]:
float(integral.subs(x, 5) - integral.subs(x, 1))

Out[2]:
$$41.333333333333336$$

## 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 = 100000
count = 0
for i in range(N):
x = numpy.random.uniform(1, 5)
y = numpy.random.uniform(0, 5**2)
if y < x**2:
count += 1
area = 4 * 5**2
area * count / float(N)

Out[4]:
$$41.413$$

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)