Transcript Lesson2

Lesson 2: RNGs and Choosing from
distributions
•
Random number generators (RNG)
•
•
•
•
•
Choosing from distributions
•
•
•
•
•
•
General need for random numbers
Pseudo-random vs. random
Linear congruential generators
Extension of period of generator by combining several generators
Discrete
Continuous: Direct
Continuous: Rejection
Probability mixing
Multidimensional distributions
Coding event problems: FORTRAN code MCBASE.FOR
1
Basic view of MC process
•
We begin with the same black box drawing:
•
Lesson 1 was concerned with the outputs.
Lesson 2 is concerned with the inputs.
•
•
•
The quality of outputs depends on quality of inputs.
The rest of the course is about the black box
2
The quest for randomness
In the early days (1950s) of computer-based Monte
Carlo, there was work done on MECHANICAL or
ELECTRONIC random number generators to give
us “true” binary randomness (0 or 1)=digits of
fraction
•
•
•
•
Many of them “filtered” the answer to correct for
non-50/50 defects
•
•
•
•
Relay flops
Decay rates (odd or even) in
Etc.
Two successive identical bits = 1 OR two successive
different bits = 0
Closer to 50/50 but (at least) twice as slow
Abandoned as too slow, too unwieldy, and—
ultimately—undesirable
3
Pseudorandom generators (RNG)
We cannot have (nor do we really want) truly
random variables.
Any random number generator from a computer is
not truly random (since it is repeatable).
•
•
•
•
•
Referred to as “pseudorandom”
This repeatability is actually very important to us for
validation reasons. Also, this repeatability is the
basis of some Monte Carlo perturbation methods.
The most common form is that of "linear
congruential generators" which generate a series of
bounded integers, with each one using the one
before it with the formula:
in 1  ain  b (mod m)
in 1
 n 1 
m
4
Possible difficulties
•
•
•
•
Periodicity: Each number depends uniquely on the one before it.
So once a number repeats, the rest of the series will be the
same as it was before. How many unique numbers you get is
the PERIOD.
It can be shown that a "full" period of m can be assured by
choosing:
•
a(mod 4)=1
•
m as a power of 2
•
b odd.
The finiteness of the string also means that certain domains of
the problem may be unreachable. (That is, there will be “gaps”
of width 1/m.) [STARS]
We have to be careful not to depend on the later digits of a
uniform deviate: they may be less random than the early digits.
(They may be the “remnants” of repeating digits of fractions.)
5
Example: m=8, a=5, b=3, i0=1
#
0
1
i
1
0
ai+b
8
3
mod m
0
3

0
.375
2
3
4
3
2
5
18
13
28
2
5
4
.25
.625
.5
5
6
7
4
7
6
23
38
33
7
6
1
.875
.75
.125
6
Extending the period
Even with “full period” L.C.G. random number generators, we still
have the related problems:
•
•
•
The period can be no longer than m
m must be a number that can be represented as an integer:
m  2 N 1
where N is the number of bits used to store integers (usu. 16. 32, or 64)
•
Fortunately, we can solve this problem by combining the output of
several uncorrelated L.C.G.’s:
  FRAC1  2     N 
•
•
The resulting period is the PRODUCT of the individual periods
The following slide shows how KENO and MCNP implement this
7
KENO RNG (FORTRAN 77)
DOUBLE PRECISION FUNCTION FLTRN()
INTEGER X, Y, Z, X0, Y0, Z0, XM, YM, ZM
DOUBLE PRECISION V, X1, Y1, Z1
EXTERNAL RANDNUM
COMMON /QRNDOM/ X, Y, Z
SAVE /QRNDOM/
PARAMETER ( MX = 157, MY = 146, MZ = 142 )
PARAMETER ( X0 = 32363, Y0 = 31727, Z0 = 31657 )
PARAMETER ( X1 = X0, Y1 = Y0, Z1 = Z0 )
X
= MOD(MX*X,X0)
Y
= MOD(MY*Y,Y0)
Z
= MOD(MZ*Z,Z0)
V
= X/X1 + Y/Y1 + Z/Z1
FLTRN = V - INT(V)
RETURN
END
8
MCNP5 RNG (FORTRAN 95)
integer, private, parameter :: R8 = selected_real_kind( 15, 307 )
integer, private, parameter :: I8 = selected_int_kind( 18 )
...
integer(I8), private, SAVE ::
&
& RN_MULT
= 5_I8**19,
&
& RN_ADD
= 0_I8,
&
& RN_STRIDE = 152917,
&
& RN_SEED0
= 5_I8**19,
&
& RN_MOD
= 2_I8**48,
&
& RN_MASK
= 2_I8**48 - 1_I8, &
& RN_PERIOD = 2_I8**46
real(R8),
private, SAVE :: &
& RN_NORM
= 1.0_R8/2._R8**48
...
function rang()
implicit none
real(R8) :: rang
RN_SEED = iand( iand( RN_MULT*RN_SEED, RN_MASK) + RN_ADD, RN_MASK)
rang
= RN_SEED * RN_NORM
RN_COUNT = RN_COUNT + 1
return
end function rang
9
“Dimensionality”
Monte Carlo theory uses the word “dimension” in a
way that is a little foreign to normal usage, so be
ready
The dimension of a Monte Carlo algorithm is how
many random numbers have to be drawn to
generate one estimate
Examples:
•
•
•
•
•
•
•
•
Coin toss problem: 1 dimension
Our pi problem: 2 dimension
Homework #1 problem of x1+x2>1.4: 2 dimension
Neutron transport problem simulation: ?
Related to the “unit hypercube” point of view
10
“Discrepancy”
Monte Carlo theory also has the special word “discrepancy”
which is important
The discrepancy of a RNG is a measure of the biggest
(absolute value) difference between how many number
SHOULD have been chosen in any “region of RNG space”
(out of N trials) and how many were actually chosen
•
•
•
•
•
•
We use the simpler “star discrepancy” in which we restrict the “region
of RNG” to be a hypercube with the origin as the lower corner
This makes it the biggest difference between what fraction SHOULD
have been chosen in a space from the origin (0,0,…,0) to any point
(x,y,z,…) and the actual number that was ACTUALLY chosen
[blackboard example]
The fraction that should have fallen in is x*y*z…
You only have to check at actual chosen points (with and without the
point)
The expected size of discrepancy is proportional to
•
•
It is the reason for the slow convergence of Monte Carlo
1
N 11
“Discrepancy”
12
Stratified sampling
•
•
•
•
•
One can use this result to reduce variance by
forcing more order onto the random number stream
used as input
As the simplest example, it can be implemented by
FORCING the random number stream to produce
EXACTLY half the numbers below 0.5 and half
above 0.5
In effect, you are dividing the problem into two
subproblems with ½ the width of the domain each
This reduces the discrepancy (proportional to
the STANDARD DEVIATION) of each to ½ its
previous value, which cuts the variance (per
history) in half.
If you subdivide RNG into M equal divisions, the
13
variance is reduced by a factor of M
Quasi-random numbers
The logical extension of the idea of stratified
sampling is to try to make N=M
•
•
•
•
Result: Quasi-random numbers (also called “Low
discrepancy sequences”)
•
•
•
•
Each of the i=1,2,…,N random numbers is located in its
own subdomain ((i-1)/N,i/N)
This looks (and acts) in 1D like a simple Riemann choice
The trick is to make them jump around “randomly”
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.
14
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 “van der Corput”
sequences, which are found by "reflecting" the
digits of the base counting integers about their
radix point.
•
•
•
•
•
•
You pick a base.
Count in that base.
Reflect the number
Interpret as base 10.
15
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 re-ordering
results in the an (almost) 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.
•
•
•
•
•
So, if your problem has more than one decision (e.g., the
pi problem had 2—x and y coordinates) you have to let
each decision have its own RNG
16
Quasi-random numbers (4)
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—2,3,5,7,11,etc.—for
decisions. (“Halton” sequence, after UNC
professor)
(log 10 N ) d 1
Asymptotic estimate of error is
which
N
means closer to ~ 1/N
•
•
•
•
•
•
but deteriorates very quickly with dimensionality
A Java subroutine is given in the public area.
The resulting standard deviations calculated and
printed by standard MC codes are NOT accurate
17
Overview of pdf and cdf
•
Basic definition of probability distribution
function (p.d.f.):
 ( x)dx  Pr{xi lies in ( x, x  dx)}
•
And its integral, the cumulative distribution
function (c.d.f.):
x
  x     ( x' )dx'  Pr{ xi  x}

18
Overview of pdf and cdf (2)
•
Corollaries of these definitions:
x
  x     ( x' )dx'  Pr{ xi  x}

b
 b    a     ( x' )dx'  Pr{a  xi  b}
a


  ( x' )dx'  1.000

19
Mapping ->x using p(x)
•
Our basic technique is to use a unique
•
•
y->x
y= from (0,1) and x from (a,b)
We are going to use the mapping backwards
20
Mapping (2)
Note that:
•
(a)=0
•
(b)=1
•
Function is non-decreasing over domain (a,b)
Our problem reduces to:
•
Finding (x)
•
Inverting to get x(), a formula for turning pseudorandom numbers into numbers distributed according
to desired (x)
21
Mapping (3)
•
We must have:
Pr{ lies in ( ,   d )}  Pr{ x lies in ( x, x  dx)}
d   ( x)dx
d ( x )
  ( x)
dx
x
 ( x)   (a)    ( x' )dx'
a
x
 ( x)    ( x' )dx'   ( x), the C.D.F of (x)
a
22
Resulting general procedure
•
Form CDF:
x
( x)    ( x' )dx'
a
•
Set equal to pseudo-random number:
   (x)
•
Invert to get formula that translates from  to x:
x   1 ( )
23
Uniform distribution
•
For our first distribution, pick x uniformly in range
(a,b):
 ( x) 
•
~ ( x)  1, x  a, b 
~ ( x)
1

b
b
~ ( x) dx


 1 dx
a
a
1

ba
Step 1: Form CDF.
1
xa
( x)    ( x' ) dx'  
dx' 
ba
ba
a
a
x
x
24
Uniform distribution (2)
•
Step 2: Set pseudo-random number to CDF:
xa

ba
•
Step 3: Invert to get x():
x  a  (b  a)
•
Example: Choose m uniformly in (-1,1):
m  1  (1  (1))  2  1
25
Discrete distribution
•
For a discrete distribution, we have N
choices of state i, each with probability  i ,
so:
 ( x)  1 ( x1 )   2 ( x2 )    N  ( xN )
•
Step 1: Form CDF:
0, x  x1
 , x  x  x
1 1
2

x
   2 , x2  x  x3
 ( x)    ( x' ) dx'   1



N
  i  1, xN  x
 i 1
26
Discrete distribution (2)
•
•
Step 2: Set pseudo-random number to
CDF:    (x)
Step 3: Invert to get x():
 x1 , if 0     1
 x , if       
1
1
2
 2


n-1
n

1
x       xn , if   i  ξ    i

i 1
i 1


N-1

 x N , if   i  ξ  1
i 1

27
Discrete distribution (3)
•
Example: Choose among 3 states with
relative probabilities of 4, 5, and 6.
28
Continuous distribution: Direct
•
•
This fits the “pure” form developed before.
Form CDF:
x
( x)    ( x' )dx'
a
•
Set equal to pseudo-random number:
   (x)
•
Invert to get formula that translates from  to x:
x   1 ( )
29
Continuous: Direct (2)
•
Example: Pick x from:
~( x)  e  x , 0  x  
30
Homework from text
31
Homework from text
32
Homework from text
33
Homework from text
34
Homework from text
35
Homework from text
36