Introduction to Monte Carlo Integration Techniques

Overview

Teaching: 30 min
Exercises: 5 min

Questions

• How to integrate a function using a monte carlo method?

• What is the acceptance-rejection method?

Objectives

• Understand basics of Monte Carlo integration technique

(Pseudo) Random Numbers:

Random numbers are a sequence of numbers \(r_1, r_2, ...r_N\) such that there are no correlations between \(r_i\) and \(r_j\), where \(i,j \in N\). However, computationally, it is not possible to define a sequence of truly random numbers. One may define pseudorandom numbers on a computer using for example the following equation:

such that \(r_{1}\) is a number provided by the user and referred to as “seed”. The following python code shows how one may generate their own random numbers (using the linear congruent algorithm):

import math


def genrandom(ri, a, b, M):
    ri1 = (a * ri + b) % M
    return ri1

Using this function, one can generate a sequence random numbers. One notices that the random numbers tends to repeat themselves after a sequence of length M. For instance, the following program shows the sequence of numbers generated by the random number generator:

r, a, b, M = 1, 1, 1, 5
for i in range(M * 3):
    r = genrandom(r, 1, 1, 5)
    print(r)

Python comes alongwith modules, such as random, numpy, which may be used to generate random numbers. For instance, a uniform sequence of random numbers, that is, random numbers between [0, 1] may be generated by iterating over the following base program:

import random

random.random()

Monte Carlo Inegration Technique

In order to use the Monte Carlo Integration technique, we use the following observation that an integral \(\int_{a}^{b} f(x) dx\) can be written as \(\frac{b-a}{N} \Sigma_{i=1}^{N} f(x_i)\), where \(x_i\) belong to the phase space over which the integration is carried out. The following is the python implementation of the technique where we integrate the \(f(x) = 2x\) with b = 1, a = 0.

import random


def func(x):
    return 2 * x


# evaluate f(x_i) at each point and add
Npoint = 100
fi_sum = 0
for i in range(Npoint):  # we use 100 points
    fi_sum += func(random.random())

# integration result
b = 1
a = 0
N = Npoint
result = (b - a) / N * fi_sum
print(result)

Notice that the obtained result may fluctuate depending on how many number of points are used to integrate.

Acceptance-Rejection Technique

A common Monte Carlo technique in integrating a function is the acceptance-rejection method which proceeds as follows:

1. Generate a random number, x uniformly distributed between [xmin, xmax].
2. A second random number, u, is then generated between [0, fmax].

3. Accept x, if u < f(x), otherwise reject x.

   Python based examples to be added for acceptance-rejection method.

Key Points