## Search in the document preview

First Lecture: Nasir M M 1 June 23, 2012

**Room A-114,
**

**Department of Physics & Applied Mathematics,
**

**Pakistan Institute of Engineering & Applied
Sciences, PIEAS
**

**P.O. Nilore, Islamabad, 45650
**

**email: [email protected]
**

Dr. Nasir M Mirza

*Lecture One:
*

*Introduction
*

Monte Carlo Methods

Docsity.com

First Lecture: Nasir M M 2 June 23, 2012

**Monte Carlo Methods
**

• In gambling often dice, roulette, cards, etc. are

used to generate random events.

• The random number generators on the computer

are more efficient (appr. 108 random numbers / s)

During the WW II John Neumann and Stanislaw Ulam

worked on the research of behavior of neutrons and

their penetration through different materials.

There appeared an important problem, how to

determine the percentage of neutrons in one given

group, which goes through eg. watter tank.

To solve this problem they used the technique of

roulette and called it Monte Carlo. It is name of a city full

of Casinos in Monaco.

Docsity.com

First Lecture: Nasir M M 3 June 23, 2012

**CONTENTS OF COURSE
**

• What is Monte Carlo Method

• Random Numbers

• A Bit of Probability Theory

• Sampling Random Variables

• Monte Carlo integration

• Importance Sampling

• Gibbs Sampling

• Markov Chain Monte Carlo

• Variance Reduction Techniques

Recommended Text:

• M. H. Kalos, P.A. Whitlock, Monte Carlo Methods,

Vol. I, John Wiley and Sons, NY.

**Monte Carlo (MC) Method**

Docsity.com

First Lecture: Nasir M M 4 June 23, 2012

**Monte Carlo methods** are a widely used class of

computational algorithms for simulating the behavior

of various physical and mathematical systems.

They are stochastic, that is nondeterministic as

they use random numbers as opposed to

deterministic algorithms.

It is repetitive algorithms and the large number

of calculations are involved.

Monte Carlo is suited to calculation using a

computer,

A **Monte Carlo algorithm** is a numerical method used to

find solutions to mathematical problems (which may have

many variables) that cannot easily be solved, for example,

by integral calculus, or other numerical methods.

Its efficiency relative to other numerical methods increases

when the dimension of the problem increases.

**Monte Carlo (MC) Method**

Docsity.com

First Lecture: Nasir M M 5 June 23, 2012

**When is MC helpful? **
We have all the probabilities we need for the
computation.
But actually doing the computation is too hard.
For example:

• Integration in many dimensions

• Optimization in many dimensions

• Simulation

**Monte Carlo (MC) Method**

Some define the Monte Carlo method as a numerical stochastic method for statistical simulation which utilizes sequences of random numbers to perform the simulation.

The Monte Carlo method is also playing games of chance where behavior and outcomes can be used to study some natural processes.

**Stochastic Processes?**

These are sequences of states whose evolution is determined by random events. These can be generated by random numbers in computers.

Docsity.com

First Lecture: Nasir M M 6 June 23, 2012

Modelling light transport in multi-layered tissues (MCML)

Graphics, particularly for ray tracing;

Monte Carlo methods in finance

Reliability Engineering

In simulated annealing for protein structure prediction

In semiconductor device research,

Environmental science,

Monte Carlo molecular modeling

Search And Rescue and Counter-Pollution.

In computer science (Las Vegas algorithm etc)

Modelling the movement of impurity atoms (or ions) in

plasmas in existing and tokamaks

In experimental particle physics, for designing detectors,

Nuclear and particle physics codes using the Monte Carlo

method

Many Other Applications.

**Some Applications of MC Methods**

Docsity.com

First Lecture: Nasir M M 7 June 23, 2012

What is the p ?

We usually consider it to be 3.14, but that’s not

precise. Ludolf von Ceulen enumerated for 35

decimal places:

= 3.141592653589793238462643383279

How to compute the circumference and the area of

the circle?

Everyone knows (or can find it in a clever

mathematical book), that the length L and the area A

of the circle are

*rL *p2
2*rA *p

Can we estimate ? YES, we can!!! :-)

The French naturalist Comte d´e Buffon showed that

the probability that a needle of length **l** thrown

randomly onto a grid of parallel lines with distance

d( > **l**) apart intersects a line is

*d
*

*l
p
*p

2

Docsity.com

First Lecture: Nasir M M 8 June 23, 2012

**Advantages of MC Methods:
**

• Direct simulation: PDF vs. conventional differential equations

• Discontinuous, multi-value, multi-dimension functions

• Doesn't require significant memory or computation time

• Easily to modified for more complex problem, always providing a solution.

**Disadvantages:
**

• Convergence problem

• Robust – subtle effect can be lost

**The aim of Monte Carlo method
**

**Problem 1:** Generating samples {**x**j} from a given target

density *P*(**x**).

**Problem 2:** Estimating expectations of functions.

**Methods:** Importance sampling and Rejection sampling

Docsity.com

First Lecture: Nasir M M 9 June 23, 2012

Let us find area under the curve using Monte Carlo simulation procedure. The shaded area is quarter of the circle shown in Figure.

122 *yx
*

0.0 0.2 0.4 0.6 0.8 1.0 1.2 0.0

0.2

0.4

0.6

0.8

1.0

1.2

*x
2
+ y
*

*2
= 1
*

**O C
**

*Y
i
*

*X
i
*

*
**BA
**

**y(x)
**

* x
*

If we throw series of darts at random on a rectangular board having dimensions of one unit one unit, what is the probability P that the darts will hit the shaded area? The probability is given by the following ratio:

*throwndartsofnumberTotal
*

*areashadedthehittingdartsofNumber
P *

** Monte Carlo method: Area Under the Curve
**

Docsity.com

First Lecture: Nasir M M 10 June 23, 2012

0.0 0.2 0.4 0.6 0.8 1.0 1.2 0.0

0.2

0.4

0.6

0.8

1.0

1.2

*x
2
+ y
*

*2
= 1
*

**O C
**

*Y
i
*

*X
i
*

*
**BA
**

**y(x)
**

* x
*

It is related to the shaded area by following:

*OABCAreaTotalPareashadedThe *

** ii yY **

** ii yY **

hit or a miss can be easily decided by the following if

it is a miss and

if it is a hit

** Monte Carlo method: Area Under the Curve
**

Docsity.com

First Lecture: Nasir M M 11 June 23, 2012

*Step 1:** Get two random numbers for Xi
*

*and Yi representing the abscissa and
*

*ordinate of impact point of the dart such
*

*that Xi & Yi both reside in [0, 1] domain.
*

*Step 2:** Then compare the Yi with the
*

*value of yi corresponding to Xi. If the
*

*conditions are satisfied we register a hit.
*

*Step 3:** Then calculate the probability, P
*

*using above equation.
*

*Step 4:** Repeat N times the step 1 through
*

*3.
*

*Step 5:** Then compute the area under the
*

*curve. The total area of the circle with
*

*radius one and center at origin, is four
*

*times this area. *

** Monte Carlo method: Area Under the Curve
**

Docsity.com

First Lecture: Nasir M M 12 June 23, 2012

**N = 100;
**

**hit = 0;
**

**mis = 0;
**

**xs = 1;
**

**ys = 0;
**

**rand('state', 0) % this sets the generator to an initial state
**

**nframes = 10; % number of frames for movie
**

**delta = 1.0/N
**

**for i = 1: nframes
**

**for j = 1: N
**

** xs = xs - delta;
**

** ys = abs(sqrt(1- (xs*xs)));
**

** plot(xs, ys, '--')
**

** axis([0 1 0 1])
**

** hold on
**

**end
**

**for k = 1: N
**

** x = rand;
**

** y = rand;
**

** if((x*x + y*y)<= 1 )
**

** hit = hit + 1;
**

** plot(x, y, 'r:o')
**

** hold on
**

** else
**

** mis = mis + 1;
**

** plot(x, y, '*')
**

** end
**

**end
**

**m(k)= getframe;
**

**end
**

**movie(m, 30) % play movie 20 times.
**

**probability = hit/N; area = probability**

**MATLAB PROGRAM: Area Under the Curve
**

Docsity.com

First Lecture: Nasir M M 13 June 23, 2012

Fig. The hitting darts within the circle are shown

as small circles and the star symbol is used for

darts falling outside.

** Monte Carlo method: Area Under the Curve
**

Docsity.com

First Lecture: Nasir M M 14 June 23, 2012

**% Program name; monte_area_2.m
**

**% Program to find area under curve by
**

**Monte Carlo method
**

**% as a function of number of histories
**

**N = 100; delta = 1.0/N;
**

**rand('state', 0) % initialize the
**

**generator to state zero
**

**probability(1) = 1; area(1)= 3.; xn(1)=
**

**N;
**

**exact_area(1)=3.14156
**

**for i = 2: 240
**

** hit = 0; mis = 0;
**

** for k = 1: N
**

** x = rand; y = rand;
**

** if((x*x + y*y)<= 1 ), hit = hit + 1;
**

** else
**

** mis = mis + 1;
**

** end
**

** end
**

** probability(i) = hit/N;
**

** area(i) = 4*probability(i);
**

** N = N + 100; xn(i) = N;
**

** exact_area(i) = 3.14156;
**

**end
**

** plot(xn, area, xn, exact_area)**

**MATLAB PROGRAM: Area Under the Curve
**

Docsity.com

First Lecture: Nasir M M 15 June 23, 2012

**Figure 5.3 Effect of number of histories on an
estimate of area within a circle using Monte Carlo
**

**method.**

** Monte Carlo method: Area Under the Curve
**

Docsity.com

First Lecture: Nasir M M 16 June 23, 2012

**Random Number Generation
**

John von Neumann suggested the *mid-square
*method. His idea was to take the square of the
preceding random number and extract the middle
digits.

The answer is that they are not really random, but
only seem so, and are in fact referred to as *pseudo-
random *or *quasi-random; *still we call them *random,
*with the appropriate reservation.

Von Neumann's method was also slow and awkward for statistical analysis; in addition the sequences tend to cyclic, and once a zero was encountered the sequence stopped.

For instance, let us use a four-digit numbers such as 5232, we square it, obtain 27373824; Then the next number consists of the middle four digits (namely, 3738 ) and the procedure is repeated.

This raises a logical question: how can such sequences, defined in a completely deterministic way, be random?

Docsity.com

First Lecture: Nasir M M 17 June 23, 2012

The random numbers generated by this or any other method are "good" ones if

• they are uniformly distributed,

• statistically independent, and

• reproducible.

*MmodTS *

For instance, consider *27 = 7 mod 5*.

Here,* S *is 27*; T *is7 with a difference of 20.

It is evenly divisible by 5.

We may also say that the quotient *(S – T )/M* is an
integer i.e. (27 – 7)/5 = 4.

(5.2)

**Random Number Generation
**

A good method is, moreover, necessarily fast and requires minimum memory capacity.

The difference of two numbers *S *and *T* is evenly

divisible by an integer *M*, then *S* is termed as

**“ congruent to T modulo M”.** It is written as

Docsity.com

First Lecture: Nasir M M 18 June 23, 2012

If *S, T* and *M* satisfy above mentioned equation
then two further relations exist:
*1. S, T* and *M* have same remainder after division
by M. For example,
27 = 7 mod 5,

5

2 5

5

27

5

2 1

5

7

*2.* Also, the following relation exists:

* (S – T) = I ×M. *

Where, *I *is an integer.

For instance, 27 = 7 mod 5

It means that *27 – 7 = 4 ×5.*

*MmodTS *

**Random Number Generation
**

Docsity.com

First Lecture: Nasir M M 19 June 23, 2012

4 4

10 10

*S
*

*Find the least positive residue if T and M are given as
*

*10 and 4 respectively. *

*Solution:*

Let us evaluate the *S* when *T* = 10 and *M* = 4.

Thus, *S *= 10 mod 4,

S = 10 – (2)4 (truncated integer division) = 2 In the division, only the integer part of the division has been used and the decimal part of the quotient has been discarded.

**Example
**

*(S – T) = I ×M. *

Where, *I *is an integer.

Docsity.com

First Lecture: Nasir M M 20 June 23, 2012

*M
M
*

*T
TMTS
*

*n
nn
*

*n *

)(mod

The power residue is defined as

25)1(75 5

7 )7( 11

*S*,

45)9(495 5

7 )7(

2 2

2

*S
*

35)68(3433 *S
*

and *S4, S5* are 1 and 0 respectively.

Then what are the residues for T = 7 and M = 5.

SOLUTION:

**Example 2
**

Docsity.com

First Lecture: Nasir M M 21 June 23, 2012

In using power residue concept to obtain random

numbers, a random number *Xn+1* is generated using

following relation:

)2(mod1
*d
*

*nn aXX *

The *Xn *is the previous random number, *d* is number

of significant bits in a word and *a* is a constant

multiplier to obtain the largest possible sequence of

random numbers without repetition.

*d
*

*d
*

*n
n
*

*d
*

*n
*

*aX
aXaX *2

2 )2(mod

Docsity.com