Lecture Notes for Simulation

28 February 2005 - Random Numbers

Outline

• Random numbers.

• Generating random numbers.

  • Linear congruential algorithms.

• Validating randomness.

  • The uniformity test, the run test.

• Analyzing empirical results.

• Generating non-uniform random numbers.

  • Inversion, rejection

Why Random Numbers?

• Simulation is a search, and the broader the search the better.

• Even simple simulations have enormous search spaces.

• Good random sampling is an effective replacement for exhaustive search.

  • The happy accident view of searching.

• More practically, randomness is a realistic model of various phenomena.

  • And randomness is a pretty good clay for forming other, non-random behavior.

Random or Pseudo-random?

• Randomness is directly proportional to outcome uncertainty.

  • The more uncertainty, the more random the outcomes are.

• Most computer programs are deterministic, and have (or should have) minimal uncertainty.

• Writing code to generate random numbers seems a fool's errand.

• Unless

  • the code incorporates uncertainty, or

  • almost (or pseudo-) randomness is acceptable.

• True uncertainty is expensive and hard.

• Opt for pseudo-randomness and deal with it.

Requirements

• Better random numbers are

  • drawn from larger possible sets of numbers.

  • part of longer cycles.

  • have smaller autocorrelation effects.

  • more supportive of predictive statistics.

• Better random-number generators are

  • Portable; computer arithmetic makes this tricky.

  • Fast; simulations call RNGs all the time.

  • Repeatable; it should be easy to recreate a random-number sequence.

  • Clonable; simulations rely on separate random-number streams.

Generating Random Numbers

• Rely on existing RNGs after you've validated them.

• Create your own RNG based on existing RNG algorithms (and validate it).

• Never implement your own RNG algorithm.

Linear Congruential Algorithms

• Linear congruential generators are based on one of the earliest and time-tested RNG algorithms.

• X_i = (aX_{i - 1} + c) mod m, where

  m		is the modulus,	m > 0
  a		is the multiplier,	0 <= a < m
  c		is the increment,	0 <= c < m
  X0		is the seed,	0 <= X0 < m

• The choice of m, a, and c is crucial. For a maximum period, the following should hold:

  1.		gcd(a, m) = 1 (a and m are relatively prime)
  2.		(a mod p) = 1 for prime factors p > 2 of m
  3.		If (m mod 4) = 0, then (a mod 4) = 1
  4.		gcd(c, m) = 1

Validating Randomness

• Always validate unfamiliar RNGs.

• Determining a RNG's quality is tricky.

  • It's part of what makes developing RNG algorithms hard.

• A randomness test can help determine an RNG's quality.

  • "Help" because there aren't good measures of randomness quality.

• A randomness test can be either theoretical or empirical.

  • A theoretical randomness test derives results from the workings of the RNG algorithm.

  • An empirical randomness test derives results from the behavior of the generated random numbers.

• Randomness test results are suggestive, rather than definitive.

The Uniformity Test

• The histogram of outcomes from a uniform random variable should be flat.

• Flat histograms are the basis for the empirical uniformity test:

  1. Divide the unit interval [0, 1] into b buckets.

  2. Sort n random numbers into the buckets.

  3. Each bucket count should be around b/n.

  A non-uniform distribution that passes the uniformity test.

• Repeating the uniformity test with different bucket counts can help avoid these problems.

The Run Test

• A run is a maximal sequence of increasing numbers.
  25, 64, 70, 51, 44, 13, 79, 5

• Given a sequence of random numbers, the probability of a run of length k is
  1/k! - 1/(k + 1)!

• This leads to the empirical run test:

  • Given a sequence of n random numbers,

  • count the number of runs of length k,

  • skipping the random numbers that end runs.

  • The expected number of runs of length k is n(1/k! - 1/(k +1)!).

Analyzing Empirical Results

• An empirical test results in an expected result set and an actual result set.

  • A run of the uniformity test with 5 buckets and 1000 samples gives

    results	195	209	196	195	205
    expected	200	200	200	200	200

    Is this good or bad?

• A significance test (a test statistic) provides

  • a judgment about whether the results are good or bad and

  • a probability that the judgment is correct.

The Chi² Significance Test

• Given the actual and expected results from an empirical test run on a outcomes sequence, the chi² significance test (or just chi² test) produces

  1. a difference outcome drawn from a chi² distribution characterizing the difference between the actual and expected results, and

  2. the probability, under the assumption that the original outcome sequence was random, that the difference outcome came from a truly random chi² distribution.

• Notice that there are two chi² distributions involved:

  1. The one derived from the empirical test results that produces the difference outcome.

  2. The one derived from a random sequence to which the difference outcome is compared.

The Chi² Random Variable

• A chi² random variable V characterizes the difference between the actual and expected results of an empirical test.

• The empirical test data are sorted into k categories.

  • k is a parameter to the chi² distribution.

    • The distribution has k - 1 degrees of freedom.

  • Let A_i be the number of actual outcomes in category i.

  • Let E_i be the number of expected outcomes in category i.

• The value of the chi² random variable V_o is
  sum(i = 1 to k, (A_i - E_i)²/E_i)

• chi² random variables are monotone increasing.

  • If V₁ < V₂, then V₁ is closer to expected behavior than V₂.

Evaluating Chi² Random Variables

• Given V = V_o, what can you say about V_o?

• If P(V <= V_o) = 0.95, then 95% of all outcomes from random sequences are at most V_o.

  • And the probability that V_o is from a random sequence is 5%.

• Interestingly, if P(V <= V_o) = 0.05, then 5% of all random-sequence outcomes are at most V_o.

  • V_o is too perfect; it's like flipping H, T, H, T, H, ...

• The values for P(V <= V_o) come from a chi²-distribution table.

Interpreting Chi² Random Variables

• The idea is to throw away the outliers when it's safe to do so.

  • We can't say anything about the non-outliers

• Let S be a sequence of outcomes generating V_o and consider p = P(V <= V_o).

  • If p = 0.01 or 0.99, then declaring S to be non-random carries with it a 1 in 100 chance of being wrong.

  • If p = 0.05 or 0.95, then declaring S to be non-random carries with it a 1 in 20 chance of being wrong.

  • If p = 0.1 or 0.9, then declaring S to be non-random carries with it a 1 in 10 chance of being wrong.

• Repeated runs with different outcome sequences may clarify the interpretation.

Example Chi² Test

• An empirical test produced these results:

  actual	195	209	196	195	205
  expected	200	200	200	200	200

• With five categories, the chi² random-variable value V_o is

  Vo	=	sum(i = 1 to k, (Ai - Ei)2/Ei)
  	=	((195-200)2+(209-200)2+(196-200)2+(195-200)2+(205-200)2)/200
  	=	(25 + 81 + 16 + 25 + 25)/200 = 197/200 = 0.985

• With four degrees of freedom, P(V <= 0.985) = 0.09.

  • Declaring the data non-random carries a 9 in 100 chance of being wrong.

Some Theoretical Randomness Tests

• The serial correlation test for liner congruential generators.

  • Computes the correlation between consecutive LCG values.

  • Defined as 1/a(1-6c/m + 6(c/m)² + e), where e = (a + 6)/m.

• The spectral test groups sequence numbers into points in an n-dimensional space.

  • Randomness quality translates into hypergeometrical properties.

  • The spectral test not only passes good random numbers, but fails bad ones.

Generating Non-Uniform Random Numbers

• Simulations do not live by uniform distributions alone.

  • Many characteristics are more accurately captured by other distributions.

• There are many techniques by which a uniform distribution can be used to generate random numbers from non-uniform distributions.

  1. Inversion

  2. Rejection

  3. Composition

  4. Convolution

  5. Transformation

Inversion

• The cdf F(X) maps samples to values in [0, 1].

• The inverse cdf F^-1(U) maps values in [0, 1] to samples.

• For example, the exponential distribution has cdf 1 - e^-x/m with parameter m.

  • The inverse cdf for the exponential is -m ln(1 - u)

• The procedure is

  • Generate the standard-uniform random number u from [0, 1].

  • Invert u via F^-1(u) and return the result.

• Unfortunately, most distributions don't have invertible cumulative distribution functions.

Rejection

• If a distribution's inverse cdf is too complicated, use a nicer one.

• Given a pdf f(x), find a function g(x) such that

  • f(x) < g(x) for all x (the pseudo-pdf),

  • the integral G(x) of g(x) exists (the pseudo-cdf),

  • and so does the inverse G^-1(x) with area A under G^-1.

Implementing Rejection

• Perform inversion on g(x) and compensate for the difference between g(x) and f(x).

  1. Generate a uniform random number a in [0, A].

  2. Compute x₀ = G^-1(a).

  3. Pick a random point r on the line (x₀, 0), (x₀, g(x₀)).

  4. If r_y <= f(x₀), accept x₀; otherwise reject it and go to step 1.

Convolution

• The values of a distribution's random variable X may be representable as the sums of another distribution's random variables Y_i:
  X = Y₁ + ... + Y_n

• X's distribution is known as the n-fold convolution (or just convolution) of Y's distribution.

• For example, a random variable from the binomial distribution of n flips with success probability p can be represented as the sum of n random variables from the Bernoulli distribution with success probability p.

  • And don't forget that a chi-square distribution with k degrees of freedom is the k-fold convolution of the standard normal distribution squared.

Points to Remember

• Treat your random-number generators well, and they'll treat your simulation well.

• Random-number generators are guilty until proven innocent.

  • Or at least until you understand how guilty they are.

• As much as possible, use other people's work (after validation).

  • Dealing with randomness is a tricky subject.

• Have a few, quick randomness tests in your programmer's toolkit.

  • A test harness to run them is useful too.

This page last modified on 2 April 2005.