This notebook is an element of the 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 provide an introduction to the use of stochastic simulation methods.

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')
integral = sympy.integrate(x**2)
integral.subs(x, 5) - integral.subs(x, 1)
In [2]:
float(integral.subs(x, 5) - integral.subs(x, 1))

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)

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)