ISE 525 Generating Random Variables Sources of randomness
ISE 525: Generating Random Variables
Sources of randomness in a computer? • Methods for generating random numbers: – Time of day (Seconds since midnight) – 10438901, 98714982747, 87819374327498, 1237477, 657418, – Gamma ray counters – Rand Tables – CEFT generator • Closed Eye Finger Typing
Pseudo Random Numbers • Pseudo random numbers: – Generated algorithmically. – Can be replicated with original parameters.
LCGs • A linear congruential generator (LCG) represents one of the oldest and best-known pseudorandom number generator algorithms. • The generator is defined by the recurrence relation: – Xn+1 = (a Xn + c) modulo m • where Xn is the sequence of pseudorandom values, and • The period of a general LCG is at most m, and for some choices of a much less than that. The LCG will have a full iff: – (a) andc are relatively prime, – a-1 is divisible by all prime factors of m – a – 1 is a multiple of 4 if m is a multiple of 4.
LCGs continued • While LCGs are capable of producing decent pseudo random numbers, they are extremely sensitive to the choice of the coefficients c, m, and a. • Historically, poor choices had led to ineffective implementations of LCGs. A particularly illustrative example of this is RANDU which was widely used in the early 1970 s and resulted in many results that are currently in doubt because of the use of this poor LCG.
Multiply With Carry Generator • Simple to implement. • In its most common form, a lag-r MWC generator requires a base b, a multiplier a, and a set of r+1 random seed (starting) values, consisting of r residues of b i. e. {x 0, x 1, x 2 , . . . , xr− 1}, and an initial carry cr− 1 < a.
Mersenne Twisters • Mersenne primes: 243, 112, 609 − 1 • Mersenne twisters: – It has a period of 219937 − 1 – It has a very high order of dimensional equidistribution. – It passes numerous tests for statistical randomness. • Diehard tests
Diehard Tests • Tests for random number generators: – Birthday spacings: The spacings between the points should be asymptotically Poisson distributed.
Die. Hard Tests • Overlapping permutations: Analyze sequences of five consecutive random numbers. The 5! possible orderings should occur with statistically equal probability. • Parking lot test: Randomly place unit circles in a 100 x 100 square. If the circle overlaps an existing one, try again. After 12, 000 tries, the number of successfully "parked" circles should follow a certain normal distribution.
Random Parking Problem • In 1958, a Hungarian mathematician, Alfréd Rényi, presented an interesting problem: "We place at random on the interval (0, x ) unit intervals not having points in common. What is the expected number of intervals we can thus place on (0, x )? ". By assuming the unit interval a car on the road, the problem is often called "random car parking problem". Rényi (1958) solved this problem analytically and has obtained the solution • Rényi, A. (1958) On a one-dimensional problem concerning random space-filling: Publications of Mathematical Institute of Hungarian Academy of Sciences, 3, 109 -127.
Die. Hard Tests • Runs test: Generate a long sequence of random floats on [0, 1). Count ascending and descending runs. The counts should follow a certain distribution. • Overlapping sums test: Generate a long sequence of random floats on [0, 1). Add sequences of 100 consecutive floats. The sums should be normally distributed with characteristic mean and sigma.
Die. Hard Tests • Minimum distance test: Randomly place 8, 000 points in a 10, 000 x 10, 000 square, then find the minimum distance between the pairs. The square of this distance should be exponentially distributed with a certain mean. • Random spheres test: Randomly choose 4, 000 points in a cube of edge 1, 000. Center a sphere on each point, whose radius is the minimum distance to another point. The smallest sphere's volume should be exponentially distributed with a certain mean.
Die. Hard Tests – The squeeze test: Multiply 231 by random floats on [0, 1) until you reach 1. Repeat this 100, 000 times. The number of floats needed to reach 1 should follow a certain distribution. – Monkey tests: Treat sequences of some number of bits as "words". Count the overlapping words in a stream. The number of "words" that don't appear should follow a known distribution.
Multiply and Carry RNG
Approaches for RN generation • • • Inverse transform Composition Convolution Acceptance-rejection Special properties
Inverse transform methods • Suppose X is continuous with cumulative distribution function (CDF) F(x) = P(X <= x) for all real numbers x that is strictly increasing over all x. • Step 2 involves solving the equation F(X) = U for X; the solution is written • X = F – 1(U), i. e. , we must invert the CDF F • Inverting F might be easy (exponential), or difficult (normal) in which case numerical methods might be necessary (and worthwhile—can be made “exact” up to machine accuracy) • Algorithm: – 1. Generate U ~ U(0, 1) (random-number generator) – 2. Find X such that F(X) = U and return this value X
Inverse Transform Method • Uniform [a, b]
ITM • Exponential distribution:
ITM (Weibull)
ITM (Weibull)
ITM (Discrete cases) Algorithm: 1. Generate U ~ U(0, 1) (random-number generator) 2. Find the smallest positive integer I such that U <= F(x. I) 3. Return X = x. I
ITM (Discrete case)
ITM (Discrete case)
Simulation of a deck of cards
Composition method for generating random variables • This applies when the distribution function F from which the RN is desired can be expressed as a convex combination of other distribution functions, F 1 , F 2 , F 3 etc. • Specifically, we assume for all x, F(x) can be written as: • This also works if :
Composition • In general, the composition algorithm is: – Generate a positive random integer J such that • P(J=j) = p_j for j = 1, 2, … • Return X with distribution function Fj • Example: the double-exponential (or Laplace) distribution has the density f(x) = 0. 5 e|x| for all x > 0.
Laplace distribution
Composition • We can express the density as: • Where IA denotes the indicator function of the set A, defined by: • So, the algorithm for generating laplacian functions is: – First, generate U 1 and U 2 as IID U(0, 1). If U 1 <. 5, return X = ln(U 2 ), else return X = -ln(U 2 )
Example using Excel
Example 8. 4 • For 0<a<1, the right trapeziodal distribution has the density:
Convolutions • If X = Y 1 + Y 2 + Y 3 + Y 4 + Y 5 … • Use the following algorithm for generating X: – Generate Y 1, Y 2, Y 1, …, Ym IID with distribution function G. – Return X = Y 1 + Y 2 + …+ Ym • Example: – Erlang: the m-Erlang random variable X with mean can be defined as the sum of m IID Random variables with a common mean of
Acceptance-Rejection • Specify a function t that majorizes the density function f, i. e. • Now, • However, is a density function.
Acceptance Rejection • The algorithm for generating random variables using acceptance rejection is: – Generate Y having the density r. – Generate U ~ U(0, 1) independent of Y – If U <= f(Y)/t(Y), return X=Y; otherwise try again. • Example: Consider the beta(4, 3) distribution on the unit interval:
1 4 7 10 13 16 19 22 25 28 31 34 37 40 43 46 49 52 55 58 61 64 67 70 73 76 79 82 85 88 91 94 97 100 2. 5 2 1. 5 1 0. 5 0 Inverse Transform Method ?
Acceptance Rejection So, we first generate Y~U(0, 1). f(Y) = 60*Y^3(1 -Y)^2/2. 0736 If U < f(Y), return Y, else reject Y and regenerate 12 10 8 6 Series 1 4 2 0 1 2 3 4 5 6 7 8 9 10
Acceptance Rejection continued • The acceptance rejection is used mostly for generating bounded random variables – although if the variable is not bounded, it can still be used. • The choice of the majorizing function is critical for the efficiency of the method. • The majorizing function can be composed of segments of several functions – piecewise linear for instance. 2. 5 2 1. 5 1 0. 5 1 4 7 10 13 16 19 22 25 28 31 34 37 40 43 46 49 52 55 58 61 64 67 70 73 76 79 82 85 88 91 94 97 100 0
Acceptance Rejection continued • Adding a minorizing function 2. 5 2 1. 5 1 0. 5 1 4 7 10 13 16 19 22 25 28 31 34 37 40 43 46 49 52 55 58 61 64 67 70 73 76 79 82 85 88 91 94 97 100 0
Pros and cons ?
Special Properties • If Y ~N(0, 1), Y 2 has a distribution. • If Z 1 , Z 2 … Zm ~ then Y=Z 1 + Z 2 + … + Zm ~ • We exploit the special relationships between distributions to generate random variables.
Generating continuous RVs • • Uniform: Exponential: M-Erlang: Gamma
Generating Normal Variates • Normal: N(m, s 2 ) can be obtained from N(0, 1 ). – Box and Muller method: • • Generate U 1 and U 2 X 1 = sqrt(-2 ln U 1) cos(2 PI U 2) X 2 = sqrt(-2 ln U 2) sin( 2 PI U 2) What if an LCG is used for U 1 and U 2 ?
Generating Normal Variates • Generate U 1 and U 2. Vi = 2 Ui -1. • W=V 12 +V 22. • If W > 1, go back to step 1, else let • X 1 = V 1*Y 1, X 2 = V 2*Y 2 • X 1 and X 2 are N(0, 1) • Excel comparison:
Bezier Distributions • Beziers are spline curves that fit arbitrary distributions, given information about percentiles. • They can also be constructed to match moments of a given distribution.
Bezier curves • Given points P 1 and P 2 , a linear Bézier curve is simply a linear bezier curve is the straight line between the two points: • A cubic bezier is
Generating Bezier Variates
Triangular Distributions 1 0 1 2
Triangular Distributions
Special Discrete Distributions • Arbitrary Discrete Distribution: – Generate U ~U(0, 1) – Return the nonnegative integer X = 1, satisfying:
Geometric distribution • Number of failures before a success.
Generating Poisson RVs • algorithm Poisson Random Number (Knuth): init: Let L ← e−λ, k ← 0 and p ← 1. • do: k ← k + 1. Generate uniform random number u in [0, 1] and let p ← p × u. • while p ≥ L. • return k − 1.
Example of Poisson Distribution Generation
Generating Binomial RVs If Yi is geometrically distributed with parameter p, then the smallest integer x for which is binomially distributed with parameters (p, n).
Random Processes • • Markov Chains Semi-Markov Chains Renewal Processes Martingales
Properties of Markov Chains
Generating random transition matrices
Example:
- Slides: 61