Random Number Generation II

1
Random Number Generation II
2
Outline

• Recap of Lehmer RNG FAQ
• Multiple streams of random numbers
• Random number variants
• Discrete random variables
• Continuous random variables

3
How do we generate random numbers?

• Lehmers algorithm
• There are many others, but Lehmer is well
  studied, widely used
• Generate a sequence of integers in ?m 1, 2,
  m-1
• x0, x1, x2,
• xi1 a xi mod m
• Modulus m, a large, fixed prime integer e.g.,
  m231-1
• Multiplier a, a fixed integer in ?m
• an initial seed x0 ? ?m if we pick any initial
  seed x0 ? ?m, we are guaranteed this initial seed
  will reappear

4
How do we select m and a?

• Would like m to be large relative to the
  precision being used
• m 231 - 1 a good choice for 32 bit arithmetic
• Select a to be a full period multiplier RNG will
  cycle through all values in ?m (period p m-1)

// brute force algorithm to test if a is full
period mult p 1 x a // assume, initial
seed is 1 while (x ! 1) // cycle through
numbers until repeat p x (a x) m //
careful overflow possible if (p m-1) // a
is a full period multiplier else // a is not a
full period multiplier
5
Can we find all full period multipliers?

• Theorem 2.1.1 If the sequence x0, x1, x2, is
  produced by a Lehmer generator with multiplier a
  and modulus m, then
• xi ai x0 mod m, i0, 1, 2,
• Theorem 2.1.4 faster test for full period
  multiplier If a is any full period multiplier
  relative to the prime modulus m, then each of the
  integers
• ai mod m ? ?m i1, 2, 3, ... m-1
• is also a full period multiplier relative to m if
  and only if the integer i has no prime factors in
  common with the prime factors of m-1 (i.e., i and
  m-1 are relatively prime)

// Given prime modulus m and any full period
multiplier a, // generate all full period
multipliers relative to m i 1 x a //
assume, initial seed is 1 while (x ! 1) //
cycle through numbers until repeat if
(gcd(i,m-1)1) // xaimod m is full period
multiplier i x (a x) m // careful
overflow possible!
6
What about overflow?

• Lehmer RNG with modulus m, multiplier a
• Let m a q r
• q quotient ?m/a? r remainder m mod a
• Let
• ?(x) a (x mod q) - r ?x/q?
• ?(x) ?x/q? - ?ax/m?
• It can be shown
• a x mod m ?(x) m ?(x)
• ?(x) 0 if ?(x) ? ?m , ?(x) 1 if -?(x) ? ?m
• Algorithm maqr is prime, rltq, and x? ?m
• // algorithm to compute ax mod m w/o overflow
• t a (x q) - r (x / q) // ?(x)
• if t gt 0 return t
• else return (t m)

7
Are we done?

• The above only tells us how to create an RNG that
  will have a large period
• It does not guarantee the RNG output will be
  random (e.g., xi1 (xi 1) mod m not
  random!)
• Need to apply statistical tests to validate that
  the RNG gives acceptably random results
• Suggested values for Lehmer (Lemmis/Park)
• m 231 - 1 (2,147,483,647)
• a 48271
• q 44488
• r 3399
• Law/Kelton (p. 412) recommend
• m 231 - 1 (2,147,483,647)
• a 630,360,016
• Note this does not yield rltq, needed for LP
  overflow prevention

8
Outline

• Recap of Lehmer RNG
• Multiple streams of random numbers
• Random number variants
• Discrete random variables
• Continuous random variables

9
Multiple Streams

• It is often desirable to generate multiple
  streams of random numbers, one stream per
  statistical element, e.g., a server in a
  queueing network
• May be desirable for each element to use the same
  random numbers across multiple runs for
  statistical output analysis
• Implementation on parallel computer may preclude
  shared variable
• We could use a family of random number generators
• Simpler approach
• Partition a single random number stream into
  non-overlapping segments each segment becomes a
  separate stream
• First number in each segment becomes seed for
  that segment
• Be careful not to draw numbers beyond your
  segment!

x0, x1, x2, xj-1, xj, x2j-1, x2j, , x3j-1,
x3j,
Stream 1 (seed x0)
Stream 2 (seed xj)
Stream 3 (seed x2j)
Stream 4 (seed x3j)
10
Determining Initial Seeds

• Basic idea define another RNG to generate
  initial seeds for streams by jumping over j
  values in original sequence at each step
• Given Lehmer RNG g(x) ax mod m, value jltm
• Stream x0, x1, x2,
• Jump stream Gj(x) (aj mod m) x mod m
• Use Lehmer RNG with multiplier (aj mod m) to
  create initial seeds
• Jump Stream x0, xj, x2j,
• Example m31, a3, x01, j6
• Jump multiplier 36 mod 31 16
• Sequence 1,3,9,27,19,16,17,20,29,25,13,8,24,10,30
  ,28,22,4,12,
• Jump sequence (initial seeds) 1, 16, 8, 4,
• In general, if you need k streams, subdivide
  sequence into k (approximately) equal sized parts

11
Outline

• Recap of Lehmer RNG
• Multiple streams of random numbers
• Random number variants
• Discrete random variables
• Continuous random variables

12
Definitions and Notation(discrete random
variables)

• Random variable X upper case
• Specific value of random variable x lower case
• Set of all possible values ? calligraphy
• Probability density function f()
• For each value x? ? f(x) P(Xx)
• Cumulative distribution function (cdf) F()
• For each value x? ? F(x) P(Xx) ?t x f(t)
• Random variate an algorithmically generated
  realization of a random variable
• Goal create random variates for random variables
  of a certain distribution f(), using the
  function Rand() that produces a uniformly
  distributed random variate x 0ltxlt1
• Approach Inverse cumulative distribution function

13
Generating Discrete Random Variates
Random variate generation 1. Select u, uniformly
distributed (0,1) 2. Compute F(u) result is
random variate with distribution f()
Example
f(x2)

• Cumulative Distribution Function of X
• F(x) P(Xx)

14
Discrete Random Variate Generation

• Inversion method to generate random variate with
  idf F
• u Random()
• Return F(u)
• Examples
• Bernoulli(p) return 1 w/ probability p, 0 w/
  probability 1-p
• u Random()
• if (u lt 1-p) return 0
• else return 1
• Geometric(p) f(x) px (1-p) ( Bernoulli
  trials until first 0)
• It can be shown F(u) ?ln(1-u) / ln(p)? for 0
  lt u lt 1
• u Random()
• return (long) (log(1.0-u) / log (p))
• Equilikely (a,b) equally likely to select an
  integer in interval a,b
• u Random()
• return a (long) (u (b - a 1))

15
Continuous Random Variates

• Concepts for discrete random variate carry over
  to continuous
• Probability density function f()
• ?ab f(x)dx P(a x b)
• Cumulative distribution function (cdf) F()
• F(x) P(Xx) ?t x f(t)dt
• Inverse distribution function
• Let X be a continuous random variable with cdf
  F() and set of possible values ?. The idf of X
  is the function F-1(0,1)? ? defined for all u
  ?(0,1) as F-1(u)x where x ? ? is the unique
  possible value for which F(x)u.
• Generate random variate as
• u Random()
• return F-1(u)

16
Exponential Random Variates

• Exponential distribution with mean ?
• f(x) (1/?) exp(-x/?) x gt 0
• F(x) 1 - exp(-x/?) x gt 0
• idf F-1(u) -? ln(1-u) 0ltult1
• u Random()
• return -? log(1-u)
• See references for other distributions

17
Concluding Remarks

• RNG selection Use the work of others!
• Look for decent underlying theory
• Look for careful implementation (e.g., overflow)
• Make sure RNG has passed statistical tests
• Multiple Streams
• Be careful to spread sequences out sufficiently
  far to avoid overlaps
• Random number variates
• Inverse method
• Techniques exist to generate distributions youre
  like to need to use in practice