Rejection Technique-Stochastic Process-Lecture Slides, Slides of Stochastic Processes

Download Rejection Technique-Stochastic Process-Lecture Slides and more Slides Stochastic Processes in PDF only on Docsity!

Lecture: Rejection Technique

Monte Carlo Methods^ Monte Carlo Methods

docsity.com

•^

When the invertible cumulative probability distribution function

is

there the direct method is always possible,

-^

There are cases when it is not practical to

calculate

the inverse

because it may contain mathematical structure that is difficult tofind.

-^

Then Another approach is used and it is called the

rejection

method.

Rejection method

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

The efficiency of the rejection technique is defined as:

∫

b a

dxx

p

xp

1 (max

Rejection Technique^ ε

p(x)

X

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

Example1: MATLAB Program

% Program name; rejection_tech1.m % Sampling from a given pd Function using rejection % Technique ; function is 4/(pi(1+xx)), 0l_limit)&(g(kk)<=u_limit)) jj=ii; end end ibin(jj) = ibin(jj) + 1; % one more for the bin end % continued

% the last part of program –

continued

for k = 1:max_bins^ % compute mid-point^ xmid(k) = (k -

0.5)*size_bin;

% compute the expected pdf^ fx

= 4.0/(pi(1 + xmid(k)^2)); ntheory(k) = Nfx*size_bin; end plot(xmid, ibin,'r+',xmid, ntheory,'b.')

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

Then Choice for function is STEP 1: First generate two uniform random numbers. Then^ STEP 2: If

sample

x from g(x):

Accept the x; else reject it and repeat step 1. Step 2 can also be written as

; (^11)

/ 4

) ( /) ( ) (^

x x g x f x

h^

=

=

π

; 1 2

1

) (^

x

x g^

−

=

That is use sampling function, from cdf:

(^2) ξ 1

x

; (^11) 2

x +

≤ ξ

(^

)^

; 1

212

≤ +^ x ξ

Example 2: Rejection Technique

docsity.com

N = 50000 ; max_bins

= 1000;

for results.

Results for Example 2

docsity.com

Example 3 Circularly Collimated Parallel Beam The normalized probability distribution in this case is

π φ ρ ρ φρ ρ πρ φ ρ φρ

2 0

0

1

), (^

2

≤ ≤ ≤ ≤

=^

o

o

dd

d d p Where,

ρ^ is the cylindrical radius,

ρo^

is the collimation beam radius and

φ^ is

the azimuthal

angle.

element

surfaceal

differenti

dd

=φ ρρ This is a separable probability distribution of the form:

(^

2 1

ρ ρ φ ρ φρ

dp

dp

d

d

p^

o

o

d

d p

ρ ρ ρρ ρ ρ ρ^

≤ ≤

=^

0

2 )(

2

1

π φ φ π φ φ^

2 0 (^12) )( 2

≤ ≤

=^

d

d p

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

-^

Actually, both methods are equivalent mathematically.

-^

However, one or the other may have advantages in execution speed depending on other factors in the application.

-^

If the geometry is not cylindrically symmetric or all the scoring that isdone does not make use of the inherent

-^

cylindrical symmetry, then the rejection method is about twice as fast asthe direct method

-^

It is due to the fact that the trigonometric functions are not employed inthe rejection method. Example 3: Circularly Collimated Parallel Beam^ Which One is better?

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

•^

Imagine that the probability distribution function is too difficult tointegrate and invert,

-^

Then it is ruling out the direct approach without a great deal ofnumerical analysis, and

-^

that it is “spiky", making the rejection method inefficient.

-^

Many probability distributions have this objectionable character. Mixed methods The mixed method is a combination of the previous two methods. It is used when we have following problems:

docsity.com