# Introduction-Monte Carlo Methods-Modeling and Simulation-Lecture, Lecture notes for Mathematical Modeling and Simulation. Pakistan Institute of Engineering and Applied Sciences, Islamabad (PIEAS)

PDF (557 KB)
21 pages
1000+Number of visits
Description
Dr. Nasir M Mirza delivered this lecture at Pakistan Institute of Engineering and Applied Sciences, Islamabad (PIEAS) to cover following points of Modeling and Simulation course: Introduction, Monte, Carlo, Methods, Algo...
20 points
this document
Preview3 pages / 21
Slide 1

First Lecture: Nasir M M 1 June 23, 2012

Room A-114,

Department of Physics & Applied Mathematics,

Pakistan Institute of Engineering & Applied Sciences, PIEAS

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 2rA 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

• 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.

• Convergence problem

• Robust – subtle effect can be lost

The aim of Monte Carlo method

Problem 1: Generating samples {xj} 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

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:

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