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! 2 Lecture: Rejection Technique Monte Carlo Methods docsity.com 3 • 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 to find. • Then Another approach is used and it is called the rejection method. Rejection method docsity.com 6 The efficiency of the rejection technique is defined as: ∫= b a dxxp xp )( )( 1 max ε Rejection Technique p( x) X docsity.com 7 Example 1: Rejection Technique Let us consider a pdf: 10; 1 14)( 2 <<+ = x x xf π 10;.1 110 <<= ξξx The algorithm to sample such a pdf using rejection technique is following: 0.0 0.5 1.0 1.5 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 f(x ) x-value It is a monotonically decreasing function on (0, 1) range and its maximum value is 4/π. 10;4.2 221 <<= ξξπ x 1ξ;1 14.3 21 ox xIf + ≤ π Accept xo , else repeat step 1. docsity.com 8 Example 1: continued The algorithm can be rephrased as Another much more precise way is: ;.1 10 ξ=x ; 1 14.2 22 ox If + > π ξ Then repeat step 1. uniformandgenerate ;.1 21 ξξ . ;1)1(.2 1 2 12 ξ ξξ = >+ oxelse repeatIf docsity.com 11 Example 2: Rejection Technique Consider a singular probability distribution function: 10; 1 12)( 2 << − = x x xf π This function is unbounded. Straight forward rejection technique is not appropriate. Instead, we write: ( ) ( ) 10; 11 12)( << −+ = x xx xf π ; 12 1)( x xgLet − = ( ) 10;4 1 14 )( )(, <<≤ + = x xxg xfThen ππ This function is has now a least upper bound that is equal to 4/π and this can be sampled using rejection technique. docsity.com 12 Then Choice for function is STEP 2: If STEP 1: First generate two uniform random numbers. Then sample x from g(x): Accept the x; else reject it and repeat step 1. Step 2 can also be written as ; 1 1 /4 )(/)()( x xgxfxh + == π ; 12 1)( x xg − = That is use sampling function, from cdf: ;1 21ξ−=x ; 1 1 2 x+ ≤ξ ( ) ;1122 ≤+ xξ Example 2: Rejection Technique docsity.com 13 % Program name; rejection_tech2.m % Sampling from a pd Function % function is 2/(pi*sqrt(1-x*x)), 0<x<1 % random numbers from rand function % in interval [0, 1] N = 50000 ; max_bins = 1000; 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 g(i) = 0.0; r1 = rand; x = 1.0 - r1^2; r2 = rand; gg =r2*r2*(1.0 + x); if(gg <= 1.0) g(i) = x; end end max_f = max(g); % find max value of func estimated size_bin = max_f/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 = ii*size_bin; if((g(kk)>l_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 16 Direct method: The cumulative probability distribution functions in this case are: Inverting gives: where the ri are random numbers on the range [0; 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 [0, 1]. Example 3: Circularly Collimated Parallel Beam docsity.com 17 In this technique, a point is chosen randomly within the square -1 < x < 1; -1 < y < 1. 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 . The code segment that would accomplish this looks like: Rejection method x = 2e0 * rng() - 1e0 y = 2e0 * rng() - 1e0 IF (x**2 + y**2 .gt. 1e0) goto 1 x = rho_0 * x y = rho_0 * y Example 3: Circularly Collimated Parallel Beam docsity.com 18 • 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 is done does not make use of the inherent • cylindrical symmetry, then the rejection method is about twice as fast as the direct method • It is due to the fact that the trigonometric functions are not employed in the rejection method. Which One is better? Example 3: Circularly Collimated Parallel Beam docsity.com 21 • Imagine that the probability distribution function is too difficult to integrate and invert, • Then it is ruling out the direct approach without a great deal of numerical analysis, and • that it is “spiky", making the rejection method inefficient. • Many probability distributions have this objectionable character. The mixed method is a combination of the previous two methods. It is used when we have following problems: Mixed methods docsity.com 22 Mixed methods The mixed method is a combination of the previous two methods. However, imagine that the probability distribution function can be factored as follows: p(x) = f(x)g(x) where f(x) is an invertible function that contains most of the “spikiness", and g(x) is relatively flat but contains most of the mathematical complexity. The recipe is as follows: 1. Normalize f(x) producing f(x) / fmax such that 1max)/)(( =∫ b a dxfxf docsity.com 23 2. Normalise g(x) producing g(x)/gmax such that this normalized function remains less than or equal to one in the range [a, b]. 3. Using the direct method described previously, choose an x using normalized f(x) as the probability distribution function. 4. Using this x, apply the rejection technique using normalized g(x). That is, choose a random number, r, uniformly in the range [0; 1]. If normalized g(x) is less or equal to r, accept x, otherwise go back to step 3. Remarks: With some effort, any mathematically complex, spiky function can be factored in this manner. The art is the appropriate choice of normalized f(x) that leaves a normalized g(x) that is nearly flat. Mixed methods docsity.com