Rejection Technique-Computational Physics-Lecture Slides, Slides of Computational Physics

Download Rejection Technique-Computational Physics-Lecture Slides and more Slides Computational Physics in PDF only on Docsity!

Stochastic Methods^ Topic: Rejection Technique

Dr. Nasir M Mirza

Computational Physics^ Computational Physics

Email: [email protected]

Docsity.com

Lecture: Rejection Technique

Monte Carlo Methods^ Monte Carlo Methods

Docsity.com

The procedure for rejection method is following: 1.^

First find out x

max^

where probability is maximum.

2.^

Scale the probability distribution function by its maximum value obtaining a new distribution function,

f ( x ) =

p ( x

) /p (

x max

which has a maximum value of 1 which occurs at

x^ =

x max

Rejection Technique •^

Clearly, this method works only if the probability distributionfunction is not infinite anywhere;

-^

and when it is not difficult to determine the location of themaximum value.

-^

When it is not possible to determine the maximum easily, thenoverestimating it will work as well, but less efficiently

Docsity.com

5

Rejection Techniques •^ So, find the maximum value of the pdf, g(x). Then make a sort ofnormalization as below:

)^.

xg g max xh

•^

Then Choose a random number of value x uniformly in the range of [a, b] and compute h(x).

-^

Choose a second uniformly distributed random number y.

-^

If h(x)

≥^

y, then accept the value of x,

•^

otherwise return to step 1. If the point chosen falls with in the twolines below the curve, it is part of thatarea and is accepted. If it is falling abovethen^

it is not part of the area and is rejected.

0.^

0.^

1.^

1.^

2.^

2.^

1 .0 0 .8 0 .6 0 .4 0 .2Normalized g(x)0.

x -v a lu e

Docsity.com

7

Example 1: Rejection Technique Let us consider a pdf:

1 0 ; (^141) )(

2

< <

=^

x

x

xf

π

1

0 ;

. 1

1

1 0

< <

=

ξ

ξ x The algorithm to sample such a pdf

using rejection technique is

following:

0.^ 0.^ 1.^

1.4 1.2 1.0 0.8 f(x) 0.6 0.4 0.2 0.

x-value

It is a monotonically decreasing function on (0, 1) range and itsmaximum value is 4/

π.

2

2 1

x

ξ^1

; 1

1 4

. 3

2

1

x^ o

x If^

≤π

Accept x

, else repeat step 1.o

Docsity.com

Example 1: continued The algorithm can be rephrased as Another much more precise way is:

;

. 1

1 0

ξ= x

2

2

x^ o

If^

Then

repeat step 1.

uniform

and

generate

2 1

ξ ξ

2 1 1 2

ξ^ ξ ξ

+ x o

else

repeat

If

Docsity.com

10

Example 1: Rejection Technique N = 90000 ;max_bins

= 80; for results below:

0.^ 0.^ 1.^

1.4 1.2 1.0 0.8 f(x) 0.6 0.4 0.2 0.

x-value

Docsity.com

Example 2: Rejection Technique Consider a singular probability distribution function:

1 0 ; (^11) 2 )(

2

< <

− =^

x

x

xf

π

This function is unbounded. Straight forward rejection technique

is not

appropriate. Instead, we write:

(^

)^

(^

)^

1 0 ; 1

1

1

2 )(

< <

−

=^

x

x x

xf

π

; 1 2

1 ) (^

x

x g Let^

−

(^

)^

1 0 4 ; (^11) 4 )( )( ,^

< <

≤

=^

x

x

xf xg

Then

π

π

This function is has now a least upper bound that is equal to 4/

π^ and

this can be sampled using rejection technique.

Docsity.com

% Program name; rejection_tech2.m % Sampling from a pd Function % function is 2/(pisqrt(1-xx)), 0l_limit)&(g(kk)<=u_limit)) jj=ii; end end ibin(jj) = ibin(jj) + 1; % one more for the bin end**

Example 2: MATLAB PROGRAM

Docsity.com

N = 50000 ; max_bins

= 1000;

for results.

Results for Example 2

Docsity.com

Direct method:^ The cumulative probability distribution functions in this case are:

Inverting gives:

where the

ri^ are random numbers on the range [

;^ 1].

The code segment that would accomplish this looks like:

rho = rho_0 * sqrt(rand

)

phi = 2e0 * pi * rand x = rho

• cos(phi) y = rho

• sin(phi)

where

rand

is a function that return a random number uniformly on the range

Example 3: Circularly Collimated Parallel Beam [0 ,^ 1].

Docsity.com

17

In this technique, a point is chosen randomly within the square

x

<^ 1;^ - 1 <

y <

If this point lies within a circle with unit radius the point is

accepted

and the

x^ and

y^ values scaled by the collimation radius

ρ.^0

Rejection method^ The code segment that would accomplish this looks like:

x = 2e0 * rng() -

1e

y = 2e0 * rng() -

1e

IF (x2 + y2 .gt. 1e0) goto

1

x = rho_0 * x y = rho_0 * y

Example 3: Circularly Collimated Parallel Beam

Docsity.com

Example 4: Point source collimated to a planar circle The normalized probability distribution in this case is: This is a separable probability distribution of the form: The cumulative probability distribution functions in this case are:

π φ θ θ θθ θ φ π φθ φθ

2 0 0

sincos 1 2 ), (^

≤ ≤ ≤ ≤

− =^

o

d o d dd p Where

θ^ is the polar angle,

θo^

is the collimation angle and

φ^ is the azimuthal

angle.

element angle solidal

differenti

dd =φθ

θsin

Docsity.com

cos_theta

= 1e0 -

rand * (1e0 -

cos_theta_0)

theta = acos(cos_theta) sin_theta

= sin(theta) phi = 2e0 * pi * rand u = sin_theta

• cos(phi)

! u is sin(theta)*cos(phi),

! the x-axis direction cosine v = sin_theta

• sin(phi)

! v is sin(theta)*sin(phi),

! the y-axis direction cosine w = cos_theta

! w is cos(theta),

! the z-axis direction cosine x = z_0 * u/w

! x = z_0 * tan(theta)*cos(phi)

y = z_0 * v/w

! y = z_0 * tan(theta)*sin(phi)

where the

r are random numbers on the i^

range [

;^ 1].

Example 4: Point source collimated to a planar circle^ The code segment that would accomplishthis looks like:

Docsity.com