To Documents

Random Number Generators

Uniform Random Variable Generators

• Realistic simulations depend on the availability of a realistic pseudo-random number generator.
   
• Random number generators for arbitrary densities (e.g., normal, exponential) depend on having a good uniform random number generator in some interval [0, M].
   
• The earliest attempt at a random number generator is the mid square method, mentioned in a 1250 manuscript by a Franciscan Friar known as brother Edvin.
   
• Here is an R implementation of the the mid square algorithm.
   
• Use this algorithm with the seed 675,248 to generate a sequence of 100 pseudo-random numbers.
   
• Although this is a good first attempt, a histogram of these pseudo-random numbers shows that they are not uniformly distributed.
   
• The most popular uniform random number in use today is the Mersenne Twister, developed in 1997 by Makoto Matsumoto and Takuji Nishimura.
   
• The Mersenne Twister is the default uniform random generator for these languages:
   
  
  R   SAS    Ruby   Python   PHP   Matlab   C++
  
  
• Pseudocode for the Mersenne Twister algorithm
   
• R also has these algorithms available for uniform random number generation:
  

  Wichmann-Hill  Marsaglia-Multicarry  Super-Duper  Knuth-TAOCP
  
  Knuth-TAOCP-2002  L'Ecuyer-CMRG  user-supplied
  

• The seed for a random number generator is used to obtain the next value for a uniform random number generator.
   
• Use this R command to set the seed and algorithm for the uniform random number generator:
  

  set.seed(seed=3584938, kind="Mersenne-Twister")
  

• Mersenne-Twister is the default algorithm, so the following command has the same effect:
  

  set.seed(seed=3584938)
  

• Set the seed explicity with the set.seed function to obtain reproducable output from a script.
   
• Use
  

  call streaminit(seed_value);
  

  to set the seed in a SAS script for the random number generator function rand.

Normal Random Variable Generators

• The R default method for obtaining normal pseudo-random variables is the inverse CDF method. If x is a standard normal random variable and Φ is the standard normal cumulative probability function Φ(a) = P(x ≤ a), then Φ^-1(x) is a standard normal random variable.
   
• We can phrase this in terms of the R functions runif and qnorm:
  

  u <- runif(100, 0, 1)
  n <- qnorm(u, 0, 1)
  qqnorm(n)
  

• R has these algorithms available for generating normal pseudo-random numbers.
  

  Kinderman-Ramage  Buggy Kinderman-Ramage (not for set.seed)  
  
  Ahrens-Dieter  Box-Muller  Inversion (default)   user-supplied
  

Exponential Random Number Generators

• The exponential distribution models various kinds of random arrivals.  For example:
   
  
  1. Customers arriving at a grocery checkout counter.
      
  2. Automobiles arriving a tool booth.
      
  3. Phone calls arriving at a help desk.
      
  4. Print jobs arriving at a networked printer.
      
  5. HTTP GET requests received by an internet server.
      
  6. Radioactive events detected by a Geiger counter.
      
  7. Number of months until a mortgage default.
  
  
• The exponential distribution is also called the memoryless distribution because the length of time you have already waited for the arrival of the next event does not affect the amount of time you will still have to wait for that event.
   
• The probability density of the exponential distribution with arrival rate λ:
   
  
  p(t) = λe^-λt,    t ≥ 0.
  
  The height of this density can be obtained with the R function call dexp(t, lambda).
   
• The cumulative probability function of the exponential distribution with arrival rate λ:
   
  
  F(t) = 1 - e^-λt,    t ≥ 0.
  
  The value F(t) can be computed with the R function call pexp(t, lambda).
   
• Problem: Create a quantile-quantile plot to decide if times between births of babies in a hospital approximate an exponential distribution.