EP24040 - Computational Physics
Kirit Makwana
Department of Physics, IIT Hyderabad
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP24040 - Computational Physics 1 / 36
Monte Carlo Methods
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP24040 - Computational Physics 2 / 36
Motivation
The integral computation by the Trapezoidal or Simpson rule can run
into trouble due to a couple of issues
The number of computations scales as N n , where N is the number of
points in each dimension and n is the number of dimensions
Although physical space has only 3 dimensions, in many cases we
need many dimensional integrals
For ex. in statistical mechanics each particle has its own dimension
Secondly, the boundary conditions can be very complicated if they
depend on each other
For ex., if the first index boundary depends on the remaining N − 1
dimensions, then the second dimension boundary depends on the
remaining N − 2 dimensions and so on
Another problem is model-fitting, where a model can have hundreds
of paramters, and we need to scan over different values of each
parameter
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP24040 - Computational Physics 3 / 36
Monte-Carlo
The idea of Monte Carlo integration is to randomly sample the space
of integration using just a few points to get the integral
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP24040 - Computational Physics 4 / 36
Monte Carlo
A circle of radius a is embedded in a square of side 2a
Randomly put points in the square, with a uniform probability of
landing anywhere inside the square
If there are p% of points that fall inside the circle, then
Area of circle = p ∗ 0.01 ∗ (2a)2 (1)
The trick here is that the points should be uniformly distributed inside
i.e. the points have equal probability of landing anywhere inside the
square
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP24040 - Computational Physics 5 / 36
Random numbers
Random numbers do not have any deterministic behavior, but they do
have a distribution function
For ex., height of individuals cannot be predicted individually, but it
follows a Gaussian distribution with some mean and standard
deviation
For generating true random numbers, you need an actual physical
process that is random, for ex., taking input from a Geiger counter of
radioactice decay, or input from a noise fluctuations in a signal, etc.
In practice, random numbers can be generated by pseudo-random
number generators
[Link]
how-to-generate-random-numbers-in-python/
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP24040 - Computational Physics 6 / 36
Pseudo random number generator(PRNG)
Pseudo-random numbers are generated by complex algorithms that
are deterministic, but extremely complicated such that they produce
numbers that appear random for all practical purposes
Typically these algorithms start with a seed number, then apply a
complex set of mathematical operations on this number to produce a
new number
This process is repeated with the new number as a seed, thus
generating a sequence of numbers that are pseudo-random
They are pseudo-random, because the algorithm is deterministic, and
the same seed produces the same sequence of numbers
One advantage of this is reproducibility - we can get same results by
using same seed. This helps in testing and verification
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP24040 - Computational Physics 7 / 36
PRNG in python
[Link]
In python numpy we have
from numpy import random
x = random . r a n d ( 1 0 )
print (x)
x = [Link](100)
x=[Link](100, size=(5))
We can also set the seed
random . s e e d ( 5 1 6 6 )
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP24040 - Computational Physics 8 / 36
Uniform distribution
Random arrays
Uniform distribution histogram
To find the area of a unit circle x 2 + y 2 = 1, take a uniform random
distribution of x and y from -1 to 1
This can simply be generated by transforming the set of
x ∈ [0, 1) → 2x − 1
x = 0 maps to -1 and x = 1 maps to 1, while the distribution remains
uniform
To generate randomly distributed points in the square of side 2, take
two sets of random points x ∈ [−1, 1) and y ∈ [−1, 1) and then
combinations of (x, y )
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP24040 - Computational Physics 9 / 36
Uniform distribution in square box
The probability of a point lying in a region ∆x < 2 is ∆x/2
Similarly the probability of it lying in a region ∆y < 2 is ∆y /2
Then, the probability of a point lying in a small square of side ∆x is
P(A ∩ B) = P(A)P(B)
This gives ∆x 2 /4, and this is uniform everywhere inside the square
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP24040 - Computational Physics 10 / 36
Area of unit circle
#D e f i n e t h e number o f random p o i n t s t o c r e a t e , N
#C r e a t e 2 a r r a y s , x , y , w i t h N random p o i n t s from =1 t o 1
r s q = x **2+y ** 2
b = np . where ( r s q <=1.0 ,1 ,0)
f r a c t i o n = np . sum( b ) /N
area = f r a c t i o n *( area of square )
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP24040 - Computational Physics 11 / 36
Exercise
Calculate the volume of a 10D unit sphere
A 10D unit sphere is defined by the locus of points (x1 , x2 , x3 , ..., x10 )
such that
(x12 + x22 + x32 + ...x10
2
)≤1 (2)
This can be solved analytically to give
π5
V10 = (3)
120
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP24040 - Computational Physics 12 / 36
Integration
Rb
Now to calculate a f (x)dx
Take the maximum ymax and ymin of this function
Create a random distribution of numbers in the rectangle a < x < b
and ymin < y < ymax
Count the number of points inside the curve with y > 0, count the
number of points inside the curve with y < 0, subtract them and
divide by total number of points to get the area fraction
Then the integral will just be fraction multiplied by area of rectangle
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP24040 - Computational Physics 13 / 36
Exercise
Lets do the below integral with Monte Carlo method
Z 3π/2
cos(x)dx (4)
0
The integral is -1
Maximum value of 1 and minimum is -1, so take a rectangle of
(0, 3π/2) in x and (−1, 1) in y
We have to be careful of the positive and negative areas
The real power of this method is in higher dimensions
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP24040 - Computational Physics 14 / 36
Error in MC method
We can repeat this calculation M times with a different set of random
numbers
From this we will get a mean and a variance of the results
PM
˜
i=1 Ii
⟨I˜⟩ = (5)
M
M
1 X ˜
σ2 = (Ii − ⟨I˜⟩)2 (6)
M −1
i=1
The error in the calculation can be estimated as the standard
deviation σ
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP24040 - Computational Physics 15 / 36
Statistical basis
An integral can be viewed as an average over some probability
distribution function ρ(x)
Z b Z b
f (x)
I = f (x)dx = ρ(x) dx (7)
a a ρ(x)
f (x)
= (8)
ρ(x)
One way of calculating an average is to randomly sample the
distribution and calculate its average
So we take a set of N uniformly distributed random points x1 , ..., xN
N
1 X f (xi )
I = (9)
N ρ(xi )
i=1
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP24040 - Computational Physics 16 / 36
Statistical basis
For a uniform distribution then
ρ(x)∆x = (∆x)/(b − a) =⇒ ρ(x) = 1/(b − a)
This gives
N
(b − a) X
I = f (xi ) (10)
N
i=1
In a standard Riemann integral, we would take N intervals and then
∆x = (b − a)/N and we would get the exact same expression as
above , except the xi would be equi-spaced points in the intergration
domain
Monte-Carlo intergration samples random points instead of
equi-spaced points
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP24040 - Computational Physics 17 / 36
Area connection
The area approximation used in the MC integration is
m
I˜ = (b − a)(fmax ) (11)
N
Where m is the number of points that lie below the function. This
can be written as
x X N Ny
fmax (b − a) X
I˜ = H(xi , yj ) (12)
Nx Ny
i=i j=1
where
(
0 if yj < f (xi )
H(xi , yj ) = (13)
1 if yj > f (xi )
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP24040 - Computational Physics 18 / 36
Area connection
This can be written as
x N Ny
X
(b − a) X fmax
I˜ = H(xi , yj ) (14)
Nx Ny
i=1 j=1
x N
(b − a) X
≈ f (xi ) (15)
Nx
i=1
This is exactly the Monte Carlo integral
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP24040 - Computational Physics 19 / 36
Non-uniform distribution
It is not necessary that we have to use a uniform distribution
We can use other kinds of distributions that may be more appropriate
for the problem
How to obtain a distribution function of a form f (y )?
Lets say there is an invertible function y (x) which maps
x ∈ [0, 1]| → y ∈ [a, b]
Suppose we generate a set of points xi which are uniformly
distributed in [0, 1] and then map them to yi = y (xi )
What is the probability distribution function of yi ?
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP24040 - Computational Physics 20 / 36
Non-uniform distributions
Let the probability distribution be f (y )
The probability of x within an interval x → x + dx is dx
This will be equal to the probability of finding the mapped point y in
y → y + dy , which is f (y )dy
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP24040 - Computational Physics 21 / 36
Non-uniform distribution
This gives
dx
f (y )dy = dx =⇒ f (y ) = (16)
dy
This means that if we want to generate a distribution f (y ), then we
need to create a mapping y (x) such that
dy 1
= (17)
dx f (y )
Example, lets say we want to generate a distribution of random
numbers with f (y ) = sin(y ) with y ∈ [0, π/2]
R π/2
Check that 0 sin(y )dy = 1, which is a basic property of probability
distribution function
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP24040 - Computational Physics 22 / 36
Sine distribution function
Here we have to solve
dy 1
= (18)
dx sin(y )
=⇒ y = cos−1 (c − x) (19)
If we take c = 1 then x = 0 maps to y = 0 and x = 1 maps to
y = π/2
Make a histogram of this distribution and check
Now an integral can be calculated by sampling xi from an appropriate
distribution
R ρ(x) and then calculating the average
f (x)dx = ⟨f (x)/ρ(x)⟩
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP24040 - Computational Physics 23 / 36
Gaussian/Normal distribution
The normal (Gaussian) distribution is
1 (y − µ)2
1
f (y ) = √ exp − (20)
σ 2π 2 σ2
σ is the standard deviation and µ is the average
The special case of µ = 0 and σ = 1 is the standard normal
distribution N(0, 1)
It can be verified that the normalization condition is satisfied
Z ∞
f (y )dy = 1 (21)
−∞
To get the standard normal distribution, we need to solve the equation
y2
Z
1
dy √ exp − =x (22)
2π 2
However, there is no closed form expression for this integral
There are a couple of different ways to produce this distribution
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP24040 - Computational Physics 24 / 36
Box-Muller transform
The indefinite integral of the Gaussian function cannot be expressed
in any standard function form
As a result we cannot sample it by transforming a uniform distribution
Consider the Gaussian distribution with zero mean and unit standard
deviation
y2
1
f (y ) = √ exp − (23)
2π 2
Suppose we generate samples xi and yi from this Gaussian
What is the probability of the point (xi , yi ) to lie in the region
x → x + dx and y → y + dy
f (x)f (y )dxdy (24)
Convert this probability into polar co-ordinates
1
exp(−r 2 /2)rdrdθ (25)
2π
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP24040 - Computational Physics 25 / 36
Box-Muller transform
We can equivalently generate points in 2D with this probability
distribution by sampling the radial co-ordinate with a distribution
function
g (r ) = exp(−r 2 /2)r (26)
and the theta co-ordinate by the distribution function
dθ/2π (27)
The r distribution can be obtained by solving
2 /2
dx = g (r )dr =⇒ x = 1 − e −r (28)
Equivalently
p
r = − 2 ln(x) (29)
and the θ distribution can be obtained from a uniform distribution
between [0, 2π)
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP24040 - Computational Physics 26 / 36
Box-Muller Transform
The Box-Muller transform is implemented by generating two random
numbers x1 and x2 from standard uniform distribution
Transform to two new random variables
p
r = −2 ln(x1 ) (30)
θ = x2 × (2π) (31)
Then two new variables can be generated that are standard normal
distributed
y1 = r cos θ; y2 = r sin θ (32)
Non-standard Gaussian with finite µ and σ can also be generated
Exercise: Implement and verify the Box-Muller transform
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP24040 - Computational Physics 27 / 36
Random walk in one dimension
The process of diffusion can be thought of as a random walk of
particles, with each particle free to move independently
Consider a particle constrained to move in the x direction in fixed
steps of size unity
It can move in the positive or negative x direction with equal
probability
What is the probability W (m, N) that the particle reaches distance m
after N steps? (m < N)
To reach m after N steps, the particle must have taken some
(N + m)/2 steps in the forward direction, and (N − m)/2 steps in the
backward direction
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP24040 - Computational Physics 28 / 36
1D random walk
The number of ways of choosing (N + m)/2 steps out of N steps is
N N!
N+m = N−m N+m (33)
2 ( 2 )!( 2 )!
The probability of taking each particular ordering of steps is (1/2)N
Therefore, the probability function is
N
1 N!
W (m, N) = (34)
2 ( 2 )!( N+m
N−m
2 )!
If N is odd, m is odd, and similarly if N is even then m is even
Lets take the case of even N and m
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP24040 - Computational Physics 29 / 36
Average displacement
The
P rms displacement undergone by a particle after N steps will be
( Nm=−N,−N+2,.. m 2 W (m, N))1/2
Using the fact that
N
1
W (m, N) = × co-efficient of u (N−m)/2 in (1 + u)N (35)
2
We can write
N
X
N
(1 + u) = (coeff .u (N−m) )u N−m (36)
m=0
N
X
= (coeff .u (N−k)/2 )u (N−k)/2 (37)
k=−N,−N+2,...
Here we have replaced k = 2m − N. Take derivative w.r.t. u
N
N−1
X N − k (N−k−2)/2
N(1 + u) = (coeff .u (N−k)/2 ) u
2
k=−N,−N+2,...
(38)
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP24040 - Computational Physics 30 / 36
Average displacement
Take u = 1 gives
N N
X N Xk
N2N−1 = (coeff .u (N−k)/2 ) − (coeff .u (N−k)/2 )
2 2
k=−N,−N+2,.. k
(39)
From this we get
N
X
kW (k, N)2N = 0 (40)
k=−N
This gives the expected result that the average displacement is zero
Similarly, by taking the second derivative we can show that
This gives ⟨m2 ⟩ = N √
In the limit N ≫ m using Stirling’s approximation (n! ∼ 2πn(n/e)n )
r
2
W (m, N) = exp(−m2 /(2N)) (41)
Nπ
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP24040 - Computational Physics 31 / 36
Markov chain
If the step size is of length l so that x = ml and the particle
undergoes n steps per unit time then
x2
1
W (x, t) = exp − (42)
2(πDt)1/2 4Dt
where D = nl 2 /2
A random walk is an example of a Markov chain
Markov chains are processes in which the next step depends only on
the current step
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP24040 - Computational Physics 32 / 36
Metropolis algorithm
In many physical cases, we will not know the global probability
density, but only local probabilites like in random walk
In other cases, we may know the target distribution function, ϕ(x),
only upto a multiplicative constant, and there is no way to integrate
and normalize it
In such cases Markov chain simulations can help
The idea is to start from some point and make a random walk around
it, but trying to move into regions where ϕ(x) is larger
Pick a starting point x0
Take a random walk step x1 = x0 + ∆x
Define the acceptance ratio
ϕ(x1 )
a = min 1, (43)
ϕ(x0 )
Pick a random number between 0 and 1, and if it falls under the
acceptance ratio, move to x1 or stay at x0
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP24040 - Computational Physics 33 / 36
Bayes theorem
Bayes theorem relates to conditional probability of event A given
event B - P(A|B)
It states that
P(B|A)P(A) P(A|B) P(A)
P(A|B) = =⇒ = (44)
P(B) P(B|A) P(B)
Suppose we have a target probability distribution function f (x) and
we are at point a
We take a step to point b, assume f (b) > f (a), then for the
Metropolis algorithm P(b|a) = 1
And for the Metropolis algorithm P(a|b) = f (a)/f (b)
From Bayes theorem, this means
P(a) f (a)
P(a|b) = = (45)
P(b) f (b)
So we will get the target distribution function from this collection of
points a and b
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP24040 - Computational Physics 34 / 36
Examples
Create a probability distribution proportional to the following function
f (x) = − sin(x) ln(x) (46)
with x ∈ [0, 1]
Generate this distribution with the Metropolis algorithm
Example
If a salesman has to travel to four cities but can only travel to
neighboring cities
Ideally, the salesman would like to visit a city as much as it’s
population
The salesman can decide whether to go to the next city or stay in the
same city depending on their population ratios as per the metropolis
algorithm
This will give visits to a city proportional to it’s population
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP24040 - Computational Physics 35 / 36
Caveats
The step size needs to be chosen carefully, if its too narrow, you will
sample a lot of points in a small neighborhood
If the step size is too large, you may get stuck on some point for a
long time
If you select a bad starting point, then the initial distribution might
not be accurately portray the target distribution
For this reason, some amont of starting samples are rejected - this is
called as burn-in period
If the sampler is getting stuck at a point for a long time, then lag is
introduced. Lag means that samples are taken only after some n
number of iterations of the cycle
This reduces the correlation between samples
If the target distribution has distinct regions separated by regions of
very low probability, the Metropolis algorithm may not cross over into
the other regions
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP24040 - Computational Physics 36 / 36