Transcript Generating Continuous

some
Generating Continuous
Random Variables
Quasi-random numbers
So far, we learned about pseudo-random
sequences and a common method for generating
them. These will be the inputs to our Monte
Carlo processes.
Alternate source: Quasi-random numbers


–
–

They are not really random, so they have to be used
very carefully.
They are not used very much for Monte Carlo
methods in radiation transport.
But, because they are of such importance to the
general field of Monte Carlo (i.e., beyond
transport methods), you should have a taste of
what they are, how to get and use them.
Quasi-random numbers (2)


"Quasi"-random numbers are not random at all;
they are a deterministic, equal-division,
sequence that are "ordered" in such a way that
they can be used by Monte Carlo methods.
The easiest to generate are Halton sequences,
which are found by "reflecting" the digits of
prime base counting integers about their radix
point.
–
–
–
–
You pick a prime number base.
Count in that base.
Reflect the number
Interpret as base 10.
Quasi-random numbers (3)




After M=BaseN-1 numbers in the sequence, the
sequence consists of the equally spaced fractions
(i/(M+1), i=1,2,3,...M).
Therefore, the sequence is actually deterministic,
not stochastic. (Like a numerical method.)
The beauty of the sequence is that the reordering results in the an even coverage of the
domain (0,1) even if you stop in the middle.
Because you need ALL of the numbers in the
sequence to cover the (0,1) range (which is not
true of the pseudo-random sequence), it is
important that all of the numbers of the sequence
be used on the SAME decision.




For a Monte Carlo process that has more than
one decision, a different sequence must be used
for each decision. (This is why we don’t use
them.)
The most common way this situation is handled
is to use prime numbers in order—
d 1
(log
N
)
2,5,7,11,etc.—for decisions. (“Halton”
10
sequence, after UNC professor) N
Asymptotic estimate of error is
which means closer to ~ 1/N
The resulting standard deviations calculated and
printed by standard MC codes are NOT
accurate:
–
Get estimates by “oversetting”
some
Generating Continuous
Random Variables
Generating Continuous Random Variables
Suppose that X is a continuous random
variable with cdf F(x).
Furthermore, suppose that you can invert F.
Note that if U~unif(0,1), the random variable
F-1(U)
has the same distribution as X!
Generating Continuous Random Variables
Example: To generate exponential rate  rv’s:
- x
f(x)   e ,
pdf:
x0

x
cdf: F(x)    e -u du  1 - e -x ,
x0
0

-1
F (x)  -
1

ln(1 - x)
Now plug in
a uniform!
Generating Continuous Random Variables
The standard normal distribution:
1
f(x) 
e
2
1
- x2
2
,
-  x  
There is no closed form expression for the cdf!
normal
Generating Continuous Random Variables
The Box-Muller Transformation:
Let U1 and U2 be independent unif(0,1) rv’s.
Then
X1  - 2 ln U1 cos(2 U2 )
is normally distributed with mean 0 and variance 1.
So is
X2  - 2 ln U1 sin(2 U2 )
and X1 and X2 are independent!
normal
Generating Continuous Random Variables
Acceptance/Rejection Method:
• want to simulate a rv X with pdf f
• need to find another function g so that
g(x)  f(x)
x
1
• normalize g to a pdf h(x)  g(x)
c
where
c   g(x) dx
• want to simulate a rv X with pdf f
• need to find another function g so that
g(x)  f(x)
x
• normalize g to a pdf h(x) 
where c   g(x) dx
1
g(x)
c
• generate Y with pdf h
• generate U~unif(0,1)
(indep of Y)
• if U<f(Y)/g(Y), “accept” Y, set X=Y
• otherwise, “reject” Y, return to
Proof: (discrete case)


P(X  x)   P X  x and first accept happens on n th try

n 1

f(y1 )
f(yn-1 )
f(x) 


  P Y1  y 1 , U1 
, , Yn-1  y n-1 , Un-1 
, Yn  x, Un 
g(y1 )
g(yn-1 )
g(x) 
n 1 


f(x) 

  P~ no x' s ~  P Yn  x, Un 
g(x) 
n 1



f(x) 

  P~ no x' s ~  P Y  x, U 
g(x) 
n 1


Proof: (continued)

f(x)  
   P(~ no x' s ~)
 P Y  x, U 
g(x)  n1


f(x) 
   no x' s 
 P Y  x, U 
g(x) 

f(x)
 g(x) 
  no x' s 
g(x)
 f(x)   no x' s 
Proof: (continued)
So, using this algorithm,
P(X  x)  f(x)   no x' s 
Sum both sides over x:
 P(X  x)   no x' s    f(x)
x
x
1   no x' s   1

 no x' s   1
What we
wanted!
 P(X  x)  f(x)
Acceptance/Rejection Method Example:
Example:
Target Density:
f(x)  9x e-3x ,
x0
-2x-1 -2x
f2(x)
 2e9e
g(x)
e
Try
g(x)  ce-2x
c  9e
works
-1
Acceptance/Rejection Method Example:
RESULTS: 100,000 draws
Target Density:
f(x)  9x e-3x ,
x0