The Metropolis Algorithm-Stochastic Process-Lecture Slides, Slides of Stochastic Processes

Download The Metropolis Algorithm-Stochastic Process-Lecture Slides and more Slides Stochastic Processes in PDF only on Docsity!

Lecture: The Metropolis Algorithm

Monte Carlo Methods^ Monte Carlo Methods

docsity.com

Most

of

the

sampling

technique

can

sample

multi-dimensional

The Metropolis Algorithmfunctions with great difficulty.^ The metropolis algorithm is designed to do it efficiently.^ It is based on the concept of a random walk. In a two-dimensionalclassical random walk, a point is started at the origin and it isallowed to moved one unit in any direction with equal probability. However, a walker is allowed to spend more of its time in thoseregions where the function is sampled with large probability. Thenpdf can be realized. Let us find rules of such walks and develop algorithms for realsituations.

docsity.com

The step is accepted with the probability q. that is

accepted

is

it

then

q

if

Else if q < 1, then a random number

is generated. If the q >

then the step is accepted.

The Metropolis Algorithm^ It is equivalent to finding an equilibrium most quickly and efficiently.

docsity.com

Consider P(r, r’)

as the probability of going from r to r’ and g(r) as

probability distribution function at r.

For equilibrium condition,

(^

r g r r P r g r r P

It says that the rate of flow towards high density over the rate

of flow

away from high density is inverse of the ratios of the densities. Then the P(r’, r) can be written into a product:

) , ( ) , ( ) ,

(^

r r T r r A r r P ′

′

=

′

The A(r’, r) is probability of a move from r to r’ is accepted and T(r,r’) is probability of proposing such a move.

) ( ) , ( ) ( ) ,

(^

r g r r A r g r r A ′

′

=

′

]

min[

(^

q

r

r

A

and

r g r g r r

q^

When q<1, it satisfies the condition of acceptance.

The Metropolis Algorithm

docsity.com

% Program name; metropolis1.m % Sampling from a pd Function % using metropolis algorithm % random numbers from rand function % in interval [0, 1] nwalk = 5000 ;

% number of walkers

max_bins = 40; % number of bins xb = max_bins; delta = 0.

% maximum step size

rand('state', 0)

% initialize the generator to zero

for j=1:max_bins+

% initialize

ibin(j) = 0; ntheory(j)=0; end for i= 1:nwalk

% Start Monte Carlo loop

x(i) = rand;

% choose random F

g(i) = sin(pix(i)); % normalization does not matter end for i=1:nwalk*

j = int8(xbx(i) + 1.); ibin(j) = ibin(j) + 1; end for j=1: max_bins^ xmin = (j-1)/xb;^ xmax = j/xb;^ ntheory(j) = nwalk(cos(pixmin)-cos(pixmax))0.5 end*

Program in

MATLAB

docsity.com

kloop = 0; for k=1:

% loop it for 500 times

k for i=1:nwalk

% for each walker Metropolis algo

xt = x(i) + delta(2.rand - 1.0);**

% propose a new value

if( (xt <= 1.0) & (xt > 0.0) )

% it must be in the range

gt = sin(pixt) ;*

% compute proposed pdf

q = gt/g(i) ; qch = rand;

if(q >= qch)

% accept or reject it

x(i) = xt; g(i) = gt; end end end

if(mod(k,10)==0)

% bin every 10 steps

for j=1:max_bins+

% start binning

xx(j) =j; ibin(j)=0; end for i=1:nwalk^ j = int8(xbx(i) + 1.);*^ ibin(j) = ibin(j) + 1; end kloop = kloop + 1; plot(xx, ibin,'r+',xx, ntheory,'b.') m(kloop)= getframe; end end

Program inMATLAB-continued

docsity.com

Direct error estimation

Assume that the calculation calls for the simulation of

N

particle

histories. ‰

Assign and accumulate the value

xi for the score associated with

the

i'th history, where

1 < i< N.

Assign as well the square of the

score

2 x

for i’th history.

Calculate the mean value of

x:

Assume that

x

is a quantity we calculate during the course of a Monte

Carlo simulation,

i.e.

a scoring variable or simply a “score" or a \tally".

The output of a Monte Carlo calculation is usually useless unless wecan ascribe a probable error to it. The conventional approach to calculating the estimated error is as follows:

∑

=

N i

i x

N

x^

1

1

docsity.com

∑

∑

=

=

− − = − − =

N i

i

N i

i

x^

x x N x x N

s^

1

2

2

1

2

2

) ( 1 1 ) ( 1 1

It is the error in <

x>,

we are seeking, not the “spread" of the distribution

of the

x

. Report the final result as i

x

x> ± s()

Estimate the variance associated with the distribution of the

x

: i

The estimated variance of <

x>

is the standard variance of the

mean:

s^ N

x

s^

2 x

2

)

(^

= > <

Remarks: The true mean and variances are not available to us in general. We mustestimate them. The estimated mean <

x>

calculated in above equation is an

estimate for the true mean and the estimated variance

s

2

calculated in

above equation is an estimate for the true variance

σ

2

Direct error estimation

docsity.com

Discrete Distribution Sampling Estimate the variance associated with the distribution of the

x

: j

∑

∑

=

=

− − = − − =

n j

j

n j

j

x^

x x n x x n

s^

1

2

2

1

2

2

) ( 1 1 ) ( 1 1

The estimated variance of

x

is the standard variance of the mean:

Report the final result as

x

x> ± s()

s^ n

x

s^

2 x

2

)

(^

= > <

Remarks: We use above equations with

n

fixed (at say, 30) because it gives a

reasonable estimate of the error in

x

Any large number will do, as long as we are within the range of applicability ofthe Central limit theorem. There is some evidence that the calculated statisticdepends weakly on the choice of

n

. Therefore, it is important to report how

your statistics were done when you publish your Monte Carlo results.

docsity.com

15

Combining Errors of Independent Runs For

m

independent Monte Carlo runs, it is easy to derive the following

relation:

k

m k

k^

x

N N

x^

∑=

⎞⎟ ⎠

⎛⎜ ⎝

>=

<^

1

Where,

is

the average value of x for the k’th run and the N

k^

is

the number of histories in the k’th run. The total number of historiesare :

k x^ ∑

=

m k

k N

N^

1

Then assuming the first order propagation of independent errors,

it

is also easy to derive that

2

1

2

2

x^ k

m k

k

x^

s

N N

s^

∑

=

⎞⎟ ⎠

⎛⎜ ⎝

=

when m = 2

2 2 2 2 2 1

2

1

2

2

1 1

2

1

x

x

x^

s

N N

s

N N

s

N

N

N

x

N N

x

N N

x

⎞ ⎟ ⎠

⎛⎜ ⎝

⎞⎟ ⎠

⎛⎜ ⎝

=

⎞⎟ ⎠

⎛⎜ ⎝

⎞⎟ ⎠

⎛⎜ ⎝

docsity.com