Direct Sampling in Monte Carlo Methods for Computational Physics, Slides of Computational Physics

Download Direct Sampling in Monte Carlo Methods for Computational Physics and more Slides Computational Physics in PDF only on Docsity!

Stochastic Methods

Topic: Direct Sampling

Dr. Nasir M Mirza

Computational Physics^ Computational Physics

Email: [email protected]

Lecture Eight:

Direct Sampling

Monte Carlo Methods^ Monte Carlo Methods

There are two types of Random Sample Populations:

Sampling with replacement

Sampling without replacement

Sampling Terminology & Concepts^ Another method is called Stratified random sampling: it divides the population

Sampling Random Variables Monte

Carlo

calculation

is

a

numerical

stochastic

process.

It

is

usually required that random variables be drawn from distributionfunctions that derive the process.

For example, the integral is given.

∫^

dx

x

g

x

f^

Its Evaluation requires that values of x must be drawn using

f(x)

and

the

value

of

average

of

g(x)

over

such

set

of

x

values

be

calculated. This is the process of “sampling x from f(x)” and it is an essential techniquefor any Monte Carlo simulation.^ We may also say that the sampling process is an algorithm thatcan produce a sequence of values of x (random variables) suchthat for any W that belongs to a universal set:

. 1

) (

}

{^

≤

= Ω ∈

∫^

dx x f

x P

k

Example 1: Direct Sampling Techniques For univariant

distributions, there is a general inversion technique

that we describe here.

x

x

x

x

x F

Suppose that x is uniform variable and its

cumulative distribution

is

F(x):

Then

for (0, 1) range, the cumulative distribution function for

x is

determined by solving the following equation:

uniform

is

where

x

F

ξ^

(^

( )

uniform is

where

F

x

ξ^

(^1) −

=

Example 2:

Direct Sampling

Let us perform the sampling of function f(x).

Suppose that the given

pdf for y is an

exponential probability distribution.

∞ < <

−

=^

y

y

y f^

0 )

exp(

) (

λ

λ

(^

y F Then the cumulative distribution is The

F(y)

represents the probability that a random variable from

f(y)

lies between 0 and y.

λ^

=^ ∫

exp( 1 )

exp(

) (

0

y

y

y F

y

log(

log(

exp(

y

y

ExponentialDistributionSampling

0

1

2

3

4

5

80 60 40 20 0 140 120 100

N = 10000No. of Bins = 500lambda = 1.0f(y) = lambdaexp(-lambday)

f(y)

Monte Carlo MethodTheory y

Example 2: N = 10000, No. of Bins = 500, Lambda = 1.

Direct Sampling

0

1

2

3

4

5

5 0^0 2 5 0 2 0 0 1 5 0 1 0 0

Example 2: N = 10000, No. of Bins = 500, Lambda = 2.

Direct Sampling

ExponentialDistributionSampling

N = 90000 ;

max_bins = 200;

a = 10.0; b = 2.0; nm = 2; ss = a./b; rand('state', 0) % initialize the generator to zero^ for j=1:max_bins % initialize ibin(j) = 0; xmid(j) = 0; ntheory(j)=0; end for i= 1:N

% Start Monte Carlo loop

yy = rand;

% choose random F

power = -1./(nm - 1.0); ff(i) = ss.(yy.^power - 1.);*

% calculate function

end max_ff = 5; % find max value of func estimated size_bin = max_ff/max_bins; % find size of bins for kk=1: N l_limit=0; u_limit=0; for ii=1: max_bins l_limit = (ii-1)size_bin; u_limit = iisize_bin; if((ff(kk)>l_limit)&(ff(kk)<=u_limit)) jj=ii; end end ibin(jj) = ibin(jj) + 1; % one more for the bin end for k = 1:max_bins^ % compute mid-point^ xmid(k) = (k - 0.5)size_bin;^ % compute the expected pdf^ ddd =(a + bxmid(k)).^nm ;^ fx = (nm - 1.).b.(a.^(nm -1.))./ddd^ ntheory(k) = 2Nfxsize_bin;*

% no of point by theory

end^ plot(xmid, ibin,'r+',xmid, ntheory,'b.')

The Example 3:Direct Samplingprogram inMATLAB

0

1

2

3

4

5

0 1000900800700600500400300200100 f(x)

Monte Carlo MethodTheory x

Parameters Used: a =10,

b = 1,

n = 2

N = 10000, No. of bins

=

200

Range for x

== (0, 5) (^

)

.

1

−

x n

bx a

ba n

x f^

n n

Example 3:

Direct Sampling

Transforming back,

(^

)^

2

2

r

u

where

du

e

drr

e^

u

r^

=^

−

−

)

2 cos(

ln

2

1

πξ

ξ − = x

2

=^

u

r

with

; 1

0 ; 1

0

2

1

<

<

< <

and

θ π^

d

and

2

θ^

sin

cos

r

x

r

x^

The sampling functions are then given as

)

2 sin(

ln

2

1

πξ

ξ − = y

Example 4: Sampling a Gaussian Distribution This method is called Box-Muller method and gives two independentGaussian variables. It is slow method.

0

1

2

3

4

5

0 400 300 200 100

Example 4: Sampling a Gaussian Distribution Using a MATLAB program, we have following:^ Parameters Used:^ N = 10000,^ No. of bins = 200^ Range for x = (0, 5)

% Program name; ex_sampling0.m % Sampling from a pdf: f(x) = 2x % in interval [0, 1] % random numbers in [0, 1]^ N = 10000 ;^ Max_bins = 40; % initialize the generator to zero rand('state', 0)^ for j=1: Max_bins+

% initialize

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

% Start Monte Carlo loop

xy = rand

% choose random F

rr = sqrt(xy);

% calculate function

jj = int8(40.0rr + 1.)*

% calculate bin index

ibin(jj) = ibin(jj) + 1 ;

% one more for the bin

end

for k = 1: Max_bins

%

compute mid-point

xmid(k) = 0.025(k-1) + 0.0125 fx = 2.xmid(k)**

% compute the expected pdf

ntheory(k) = Nfx0.**

% no of point by theory

end plot(xmid, ibin,'r+',xmid, ntheory,'b')

Example 5: Sampling a pdf

0 2500 2000 1500 1000 500

No of histories = 50000,No. of Bins

= 40

pdf =

f(x)

= 2x, [0, 1]

f(x)

x-value

Monte Carlo resultsTheoretical pdf

We developed a MATLABprogram for sampling from apdf

with the following form:

. 1

0

2 ) (^

≤ ≤

=

x

x

x Example 5: Sampling a pdf f^ For about 50 thousandhistories the resultsfrom Monte Carlosimulations startedmatching very closelyto the actual pdf.