The Metropolis Algorithm and Error-Computational Physics-Lecture Slides, Slides of Computational Physics

Download The Metropolis Algorithm and Error-Computational Physics-Lecture Slides and more Slides Computational Physics in PDF only on Docsity!

Stochastic MethodsTopic:

The Metropolis Algorithm

Dr. Nasir M Mirza

Computational Physics^ Computational Physics

Email: [email protected]

Docsity.com

Lecture: The Metropolis Algorithm

Monte Carlo Methods^ Monte Carlo Methods

Docsity.com

4

Suppose we wish to develop a pdf, g(r) in three dimensions.

For a walker say at some point, r, we propose a random step (

to some new point r’.

If

r^

is

a

three-dimensional

variable

for

M

particles,

then

the

movement is given as

}. 1 , 1 {

}; 1 , 1 {

}; 1 , 1 { − Δ

= ′

− Δ

= ′

− Δ

= ′

i

i

i

i

i

i

z

z

y

y

x

x

The Metropolis Algorithm^ 4.

After proposing this change, the value of the desired pdf

is

computed at new and old positions, giving g(r’) and g(r) andthe ratio is formed:

r

g

r

g

q^

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

Results of Program for a pdf

Example: As a simple example wehave sampled the following pdf: nwalk

= 5000 = number of walkers

max_bins

= 40 = number of bins

Loops = 200 = for repeated iterations

sin(

(^

=^

x

x

x

g

The results are shown above for bothMC simulations and theoreticalvalues.

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

Error Estimations

In Monte Carlo Method

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

Batch statistics error estimation It may require a lot of effort to stop after each history and compute the

x

's i

and

x

1. The “direct error estimation" should be considered first and if it is not

In many cases, the estimation of means and variances using the methodsdescribed in the previous section is not feasible.feasible, then an alternative that is nearly as good is “batch statistics errorestimation". We present a procedure here^ Split the

N

histories into

n

statistical batches of w=

N/n

histories each. Note

that

n

must be “large enough". A standard choice is

n

= 30. The

accumulated quantity for each of these batches is called

sum

as a sum forj

the j’th batch. The sum

isj

∑=

w i

i

j^

x

sum

1

Then compute the mean value of x

∑=

n j

j

sum

N

x

x^

1

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

The Error^ Remarks:^ Thus, by sampling the exponential three times and adding the threevalues one can sample this pdf. In general the sum of n suchvariables of exponential pdf gives a

az

n

n

e

z

n

a

z

w

− −

1

This method of combining errors effectively increases the value of

n, the

number of statistical batches used in the calculation. In view of the fact thatthe calculated statistics are thought to depend weakly on

n, it is preferable

(but only marginally so for the sake of consistency) to combine the

x^ i

's (the

raw data) into the standard number of statistical batches.^ This is easy to do by initializing the data arrays to the results of theprevious run before the start of a new run.

Docsity.com