Lecture 07 - Random Numbers

Pramod Parajuli
Simulation and Modeling, CS-331
Chapter 7
Uniform Random Generators
Testing of Uniform Random Generators
Methods of Generating Non-Uniform
l
Variables
a
ep
Inversion, Rejection, Composition
itn
1
(C) Pramod Parajuli, 2004
cs
Source: www.csitnepal.com
1
Introduction
Random numbers
R1, R2, R3, . . . .
Must have two essential properties;

- Uniformity
- Independence
Continuous random variable – Random variable X is

continuous if its sample space is a range or
collection of ranges of real values
a l
ep
itn
2
cs
2
Continuous Uniformly Distributed RNs
Each random number Ri is an independent sample

from a continuous uniform distribution between
0 and 1
1 ,0 ≤ x ≤ 1
Therefore, PDF = f ( x) = 
0 , otherwise
1
Expected value, E(R) = ∫ xdx f(x)

0 1
2 1
x
=
2 0 0 1 x
1
=
a l
Variance 2
ep
1
V ( R) = ∫ x 2 dx − [ E ( R)]2 = 1 / 12
itn
3
0 (C) Pramod Parajuli, 2004
cs
3
Generation of RNs – RNGs (RN Generators)
-Should be fast
-Should be portable
-Should have sufficiently long cycle
-Should be replicable
-Should hold ideal statistical properties (uniformity
and independence)
-Linear Congruent Generator
a l
ep
itn
4
cs
4
Linear Congruent Generator (LCG)
-Produces integers between 0 and m-1

xi+1 = (a.xi + c) mod m
i = 0, 1, 2, . .
Initial value i.e. x0 is called seed
a = constant multiplier
c = constant increment
Example;
m = 100, x0 = 27, a = 17, c = 43
For random number between 0 and 1,
a l
= Ri/m
ep
itn
5
cs
5
Multiplicative Congruent Generator (MCG)
-If c = 0
-Example;
A = 13, m = 26 = 64, x0 = 1, 2, 3, and 4
Look for the period/cycle length
Maximum period is 64/4 = 16

Here, 4 is the gap between the elements
a l
ep
itn
6
cs
6
Guidelines for selecting the values
Choose ‘m’ be one more than largest integer that

can be stored in one word of the computer being
used. Problem, how to represent such a
number? Therefore, choose (2maxword-1 – 1) i.e.
in 32 bit number system, choose (231-1)
The seed x0 must be relatively prime to m
The constant multiplier ‘a’ must be selected properly
(prime). Normally, it is selected as;
a = k.8 + 3
or a = k.8 + 5
al
or a = 2max wordsize/2 + 3
ep
or a = 65,539
itn
7
cs
7
Testing Uniform Random Generators
The general idea is to eliminate undesirable

characteristics
Empirical tests – statistical tests, generates uniform
random numbers for tests
Theoretical tests – uses numerical parameters but do
not generate numbers
Frequency test (uniformity test)

Runs test
Autocorrelation test
Gap test
a l
Poker test
ep

itn
8
cs
8
Frequency Tests (uniformity test)
Counts how often numbers in a given range

occur in the sequence to ensure that the
numbers are uniformly distributed
Example;
If there are 5,000 3-digit numbers from 000 to 999, then
we expect that there are about 500 numbers in the
range of 00 to 99. If there are exact 500 numbers in
each 10 segments, then the system is not a random
one
Therefore, how much deviation should we allow from

expected value of 500 and still accept the sequence as
a l
uniformly distributed?
ep
Use chi-square (goodness of fit) test
itn
9
cs
9
Chi-square χ2 test
Statistical test
Checks how well certain observed data fit the
theoretically expected data
First divide the observed data (here, the random
numbers) into k non-overlapping classes, k must
be >= 3
Then count Oi, the number of times the
observed data falls in each class i, i = 1, 2, 3,
4, . .
Then measure how far the observed frequency
a l
deviates from the expected
− 2
ep
k
(O E )
chi 2 = ∑ i i
itn
i =1 Ei 10
cs
10
Chi-square χ2 test
k
(O − E ) 2
chi 2 = ∑ i i
i =1 Ei
Oi is the observed number of sequences of value xi and
Ei is the expected number of occurrences of xi
This is showing the square of the differences the
observed values and the expected values (divided by the
expected values to normalize it)
Now consider two hypotheses:
NULL Hypothesis, H0 : Y matches distribution of X
Hypothesis H1: Y does not match distribution of X
If χ2 is too large, we reject the null hypothesis (i.e. the

distributions do not match)
l
a
If χ2 is small(ish) we do not reject the null hypothesis
ep

(i.e. the distributions may match)
itn
11
cs
11
Chi-square test
If the observed data depart more from the expected
value, the value of chi2 will be larger
In our example, Ei = 500 for i=1, 2, 3, . . . , 10
i Range Oi no. of observed No. of expected (Oi-Ei)2 (Oi-Ei)2/Ei
occurrences occurrences
1 000-099 468 500 1024 2.048
2 100-199 519 500 361 0.722
3 200-299 480 500 400 0.8
4 300-399 495 500 25 0.05
5 400-499 508 500 64 0.128
6 500-599 426 500 5476 10.952
7 600-699 497 500 9 0.018
a l
8 700-799 515 500 225 0.45
ep
9 800-899 463 500 1369 2.738
itn
10 900-999 529 500 841 1.682
12
cs
12
Chi-square test
Chi2 = 19.588
For large values, the hypothesis is rejected if
Χ >Χ2 2
k −1,1−α
where Χ 2 k −1,1−α is upper critical point

3
 2 2 
Χ 2
k −1,1−α ≈ (k − 1)1 − + z1−α 
 9(k − 1) 9(k − 1) 
Where Z1−α is the upper 1−α critical point of the N(0,1)
distribution
a l
Degree of freedom v = k – 1, where k is no. of sets
ep
α = Pr(reject H0 | H0 true), literally – level of significance
itn
13
cs
13
Chi-square test
Let’s consider, probability of neglecting H0 is 0.05
Then,
α = 0.05 α = 0.1
1 − α = 0.95 1 − α = 0.9
k = 10 k = 10
k −1 = 9 k −1 = 9
?? ??
Are we going to reject the hypothesis or not?
a l
ep
itn
14
cs
14
Kolmogorov-Smirnov Test
Another kind of frequency test
Compares continuous cdf, F(x), of the uniform
distribution to the empirical cdf, SN(x), of the sample of
the N observations
F ( x ) = x, 0 ≤ x ≤1
If sample from the random-number generator is R1, R2, .
. . , RN, then the empirical cdf, SN(x) is defined by
SN(x)= number of R1, R2, . . . , RN which are ≤ x

N
As N becomes larger, SN(x) should become a better
approximation to F(x), provided that the null hypothesis
a l
is true
ep
itn
15
cs
15
This test is based on the largest absolute deviation
between F(x) and SN(x) over the range of the random
variable
Steps
1. Rank the data from smallest to largest. E.g.
R(1)<=R(2)<=R(3)<= . . . . <=R(N)
2. Compute
+ i 
D = max  − R( i ) 
1≤i ≤ N N
 
 i − 1
l
−
D = max  R(i ) −
a

ep
1≤i ≤ N
 N 
itn
16
cs
16
3. Compute D = max(D+, D-)
4. Determine the critical value Dα from Kolmogorov-
Smirnov Critical values table
5. If D > Dα, the null hypothesis is rejected, otherwise
accepted
Example 7.6, page 267
a l
ep
itn
17
cs
17
Runs test – runs up and runs down
Direct test of independence (only) assumption
Chi-square test might show some numbers are uniformly
distributed but the ordering matters in case of
independence
Look at example on page 270 of Discrete Event System
Simulation, Banks, Nicol, Nelson
The general idea is to look at the patterns of runs up and
runs down
E.g.
0.08 0.18 0.23 0.36 0.42 0.55 0.63 0.72 0.89 0.91
One run – run up
0.08 0.93 0.15 0.96 0.26 0.84 0.28 0.79 0.36 0.57
a l
Nine run – five up, four down
ep
itn
18
cs
18
The run pattern is likely to be somewhere between two
extremes
If N is the number of numbers in the sequence, the
maximum number of runs is N-1 and minimum number
of runs is one
Now, if a is total number of runs in truly random
sequence, the mean and variance of a are given by
2N −1
µa =
3
16 N − 29
σ a2 =
90
For N>20, the distribution of a is approximated by
a l
normal distribution
ep
N (µa ,σ ) 2
itn
a 19
cs
19
Therefore, standardized normal test parameter
a − µa
Z0 =
σa
a − [(2 N − 1) / 3]
Z0 =
(16 N − 29) / 90
Where, Z0 is approximately equal to N(0,1)
Therefore, the hypothesis can not be rejected if
− Zα / 2 ≤ Z 0 ≤ Zα / 2
Class-demonstration;
Example 7.8, page 272 α /2 α /2
a l
ep
itn
20
− Zα / 2 Zα / 2
cs
20
Runs test – runs above and below the mean
Runs up and down don't indicate any absolute

organization of the data (or lack thereof)
We may also want to test the number of runs
(i.e consecutive values) above and below the
mean
Too many or too few is not a good sign
Calculation for mean and variance is more
complex, but we handle this test in more or less
the same way as the runs up and runs down test
Don't want the our number to be too far from the
mean
a l
Use the standard normal distribution to test
ep
itn
21
cs
21
If n1 and n2 is the number of individual observation above
and below the mean and b is total number of runs,
2n1n2 1
mean µb = +
N 2
2n1n2 (2n1n2 − N )
variance of b σb =
2
N 2 ( N − 1)
If n1 or n2 is greater than 20, b is approximately normally
distributed
a l
ep
itn
22
cs
22
Let’s define standardized variable Z0 for testing
b − (2n1n2 / N ) − 1 / 2
Z0 = 1/ 2
 2n1n2 (2n1n2 − N ) 
 N 2
( N − 1) 
 
Failure to reject the hypothesis of independence
occurs when − Zα / 2 ≤ Z 0 ≤ Zα / 2
where α is the level of significance
a l
Class-demonstration
ep
itn
23
cs
23
Runs test – length of runs
Given a sequence of random numbers, we can expect a
certain number of runs of each length 1, 2, 3, …
If the run lengths in our empirical tests differ too much from
the expected distribution, we can reject the generator
Idea is to determine the empirical distribution of runs of
various lengths, then compare it with the expected
distribution given random data
The distributions are compared using the Chi-Square test
Sample from book, page 274

In the example, we saw that run with length two is
occurring. In reality, it must be random
a l
ep
itn
24
cs
24
As in chi-square test for the samples, the runs

themselves should also poses uniformity
Determined by; 2 L [Oi − E (Yi )]
2
χ0 = ∑
i =1 E (Yi )
where,
Yi = number of runs of length in a sequence of
N numbers
L = N – 1 for runs up and down
L = N for runs above and below the
mean
a l
ep
Derivations – on blackboard
itn
25
cs
25
If null hypothesis of independence is true, then

χ02 is approximately chi-square distributed with
L-1 degrees of freedom
Example 7.10 – page 275
a l
ep
itn
26
cs
26
Tests relationship of items m (lag) locations

apart
Autocorrelation is determined for
Ri, Ri+m, Ri+2m, . . . ., Ri+(M+1)m
The value of M is the largest integer that
satisfies;
i+(M+1)m <= N
where N is total number of values in the sequence
A nonzero autocorrelation implies a lack of
independence. Let
l
H 0 : ρ im = 0
a
ep
H1 : ρ im ≠ 0
itn
27
cs
27
For large value of M, the distribution of the

estimator of ρ im ( ρ̂ im ) is approximately normal if
the values Ri, Ri+m, Ri+2m, . . . ., Ri+(M+1)m are
uncorrelated
Therefore;
ρˆ im
Z0 =
σρˆ im
where
1 M 
ρ im =
ˆ  ∑
M + 1  k =0
Ri + km Ri + ( k +1) m  − 0.25

a l
13M + 7
σρˆ im =
ep
12( M + 1)
itn
28
cs
28
Do not reject the null hypothesis of

independence if − Zα / 2 ≤ Z 0 ≤ Zα / 2
If ρ im > 0 , the subsequence is said to exhibit

positive autocorrelation. High random numbers
in the subsequence followed by high, and low
followed by low)
If ρ im < 0 , the subsequence is said to exhibit
negative autocorrelation. Low random numbers
tend to be followed by high ones, and vice versa
For desired property of independence, which
l

a
ep
implies zero autocorrelation
itn
29
cs
29
Gap test
The gap test is used to determine the

significance of the interval between the
recurrences of the same digit
The gap of length x occurs between the
recurrences of some specified digit
E.g.
4, 1, 3, 5, 1, 7, 2, 8, 2, 0, 7, 9, 1, 3, 5, 2, 7, 9,
4, 1, 6, 3, 3, 9, 6, 3, 4, 8, 2, 3, 1, ..
For digit 3, the first gap is 10 and second gap is 7.
a l
ep
P(gap of 10) = P(no 3), . . . (P no 3) P(3)
itn
30
cs
30
Gap test
The probability of the digits not being 3 is 0.9

Therefore, P(gap of 10)= (0.9)10(0.1)
The theoretical frequency distribution for
randomly ordered digits is given by
x
P(gap <= x) = F(x) = 0.1 ∑ ( 0
n =0
. 9) n
= 1 − 0. 9 x +1
Steps
1. Specify the cdf for the theoretical frequency
distribution given by P(gap<=x) based on the
selected class interval width
l
2. Arrange the observed sample of gaps in a cumulative
a
ep
distribution with these same classes
itn
31
cs
31
Gap test
Steps
3. Find D, the maximum deviation between F(x) and
SN(x) as in D =| F ( x) − S N ( x) |
4. Determine the critical value, Dα , from Kolmogorov-
Smirnov Critical Values table for the specified value of
α and the sample size N
5. If the calculated value of D is greater than the
tabulated value of Dα, the null hypothesis of
independence is rejected
a l
ep
itn
32
cs
32
Poker test
Test for interdependence based on the frequency with which
certain digits are repeated in a series of numbers
e.g.
0.255, 0.577, 0.331, 0.414, 0.828, 0.909, 0.303, 0.001..
Here, there are three possibilities
1. The individual numbers can all be different
2. The individual numbers can all be the same
3. There can be one pair of like digits
The probability associated with each of these possibilities is
given as;
P(three different digits) = P(second different from the first) *
P(third different from the first and second)
a l
ep
= (0.9)*(0.8)
itn
= 0.72 33
cs
33
Poker test
P(three like digits) = P(second digit same as the first) *
P(third digit same as the first)
= (0.1) * (0.1) = 0.01
P(exactly one pair) = 1 - 0.72 – 0.01 = 0.27
Alternatively, the last result can be obtained as follows;
 3
P(exactly one pair) =  (0.1)(0.9) = 0.27
 2
Normally used with chi-square test.
a l
ep
itn
34
cs
34
Code snippet
/* Returns values between 0 and 1*/
/* Depends on 32 bit representation for integers */
double randu( seed )
long int *seed;
{
long int a = 16807,
mod = 2147483647; /* 2^31 - 1 */
double dmod = 2147483647.0; /* 2^31 - 1 */
*seed = *seed * a;
if ( *seed < 0 )
*seed += mod;
return( ( double ) *seed / dmod );
a l
}
ep
itn
35
cs
35
Discrete distributions
In real life, the probability of occurring one kind of event/data
has different probability than other kind of event/data
Normally, there is a need for discrete but non-uniform
distribution
Therefore, while generating a uniform random number ‘U’,
the value is compared with the values of y. If the value
falls in a interval y i < U ≤ y i + 1 (I = 0, 1, . . . 4) the
corresponding value of xi+1 is taken as the desired output
No. of items (x) Probability p(x) Cumulative prob. (y)
0 0 0
1 0.10 0.10
2 0.51 0.61
a l
3 0.19 0.80
ep
4 0.15 0.95
itn
36
5 0.05 1.00
cs
36
Discrete distributions
When the number U is being tested, only 10% of the
numbers are validated by using the first entry where as
only 61% of numbers are validated by using second
entry
For such cases, the table can be rearranged as the ordering
doesn’t matter for testing
Cumulative prob. (y) No. of items (x)
0.51 2
0.70 3
0.85 4
0.95 1
1.00 5
a l
Now, the first entry can satisfy 51% of data
ep
itn
37
cs
37
Non-uniform continuously distributed RNs
Random numbers drawn from a distribution that

is continuous & non-uniform
These methods are based on use of sequences
of uniformly distributed random numbers
The algorithm must satisfy
¾ Exactness – the random variates must follow the
desired distribution
¾ Efficient – in terms of time (setup time, marginal
execution time) and space
¾ Less Complexity
Inverse transformation method – most widely
l
a
used
ep
Rejection, Composition, Convolution etc.
itn
38
cs
38
Inverse transformation method

If Ui (i = 1, 2, 3, 4, …) are independent variables,
uniformly distributed over the interval 0 to 1, and
F-1(x) is the inverse of the cumulative distribution
function for random variable x, then the random
variables defined by xi=F-1(Ui) are the random sample
of x.
If desired PDF is
a (1 − x) (0 ≤ x ≤ 1)
f ( x) = 
 0 otherwise
Cumulative;
a l
x2
ep
F ( x) = a ( x − )
2
itn
39
cs
39
Let U = F(x)
x = 1 – (1 – U)1/2
U X
0.1009 69.23
0.3754 88.48
0.0842 66.65
0.9901 142.33
a l
ep
itn
40
cs
40
For exponential distribution

CDF  0 x<0
F ( x) =  − λx
1 − e x≥0
now, find U(0,1)
then,
F(x) = R
or
1 − e − λx = R
1
e − λx = 1 − R β= mean value
λ
− λx = ln(1 − R)
l
R is cumulative
a
 1 
ep
x =  − ln( R) 
 λ 
itn
41
cs
41
In case of discrete data, if the result falls between

two values, then upper value is taken
But in case of continuous system, any value might
occur
In such case, interpolation
is applied
The plot of uniform random
x
numbers and the input
data (x) is;
a l
ep
itn
42
0 U 1.0
cs
42
The general idea is, if a value falls in between two

values, the result is calculated by using lower
value plus portion of the output that is divided
by the input
Graphically,
output Interpolated
} interval
a l
ep
input
itn
43
cs
43
The standard curve and its points do not change
Let’s define the slope of the segment as;
xi +1 − xi (i = 0, 1, 2, 3, . . . N-1)
a = i
yi +1 − yi
now,
if Yi < U < Yi+1 (i = 0, 1, 2, 3, . . . N-1)
then,
x = xi + ai(U-Yi)
U X
0.1009 69.29
a l
0.3754 68.48
ep
0.6842 66.65 . .
itn
44
cs
44
Disadvantage of inverse-transform
1. For normal and gamma distribution, there is no
mathematical expression of F-1
2. For given distribution, it may not be the fastest way
Advantages
1. Even if the distribution is truncated, can generate the
output
2. Facilitates variance reduction that rely on correlation
a l
ep
itn
45
cs
45
Composition
Applies when the distribution function F from which we
wish to generate can be expressed as a convex
combination of other distributions F1, F2, . . .
For all x,
∞
F ( x) = ∑ p j F j ( x)
j =1
pj ≥ 0
∞
∑p
j =1
j =1
a l
Fj is distribution function
ep
itn
46
cs
46
If, pk > 0, pj = 0, j>k
then, it will be finite
∞
pdf = f(x) = ∑ p .f
j =1
j j ( x)
Algorithm
1. Generate positive random number J such that
P( J − j ) = pj for j = 1, 2, 3, 4, . .
2. Return X with distribution function FJ
a l
ep
itn
47
cs
47
Convolution
Probability distribution of a sum of two or more
independent random variables is a convolution of the
distributions of the original variables
Algorithm
Let F be the distribution function of X and G be the
distribution function of Yj
1. Generate Y1, Y2, Y3, . . . Ym, independent and
identically distributed random variables each with
distribution function G
2. Return, X = Y1 + Y2 + Y3 + . . . + Ym
a l
ep
itn
48
cs
48
Rejection method
Applied when the PDF f(x) has a lower and upper limit to
its range a and b and upper bound c.
However, the upper bound ‘c’ is not very much relevant
Steps
1. Compute the values of two independent uniformly
distributed variables U1 and U2
2. Compute, X0 = a + U(b – a)
3. Compute, Y0 = c.U2
4. If Y0 <= f(X0) accept X0 otherwise repeat with two
new uniform variates
a l
ep
itn
49
cs
49
Recall the Monte Carlo method

The probability density function is enclosed in a
rectangle i.e. if a point is accepted, ‘n’ is
incremented by 1
The curve f(x) represents probability density
function. Therefore, the area under the curve
must be 1. But,
c.(a-b) = 1
In such case, the graph must be chosen so that
c.(a-b) and f(x) both are 1
a l
ep
Therefore, is this a valid method?
itn
50
cs
50
Given, X0,
X0 is accepted if,
Y0 <= f(X0)
f(X0)
∆x
y
a l
a X0 b
x
ep
itn
51
cs
51
Let ∆x be small area. If Y0 still is inside the small rectangle
bounded by ∆x , under the curve, now output range will
be,
X0 – ∆x to X0
In such case,
P(X <= X0) = P (Y0 bounded by left of X0)
X0
∫ f ( x)dx( x0 − a )
F(X0) = a
.
c( x0 − a ) (b − a )
l
Therefore,
a
X0
ep
F(X0) = ∫ f ( x)dx
itn
a 52
cs
52
Disadvantage
Two uniform variates must be calculated
The situation will be even more complex if more
iterations need to be performed
Therefore, it is usable only if,
The probability density function is given by a
mathematical function
The PDF must be limited i.e. F(x) = 0 for < a and > b. If
not, then use inverse transformation method
a l
ep
itn
53
cs
53
Sample function
/* C implementation Park & Miller’s random function */
double random ( seed )
long int *seed;
{
long int a = 16807, /* 7^5 */
m = 2147483647, /* 2^31 - 1 */
q = 127773, /* m / a (int divide) */
r = 2836, /* m % a */
lo, hi, test;
double dm = 2147483647;
hi = *seed / q;
lo = *seed % q;
test = a * lo - r * hi;
l
*seed = test > 0 ? test : test + m;
a
ep
return( (double ) *seed / dm );
}
itn
54
cs
54

Lecture 07 - Random Numbers

Uploaded by

Copyright:

Available Formats

Lecture 07 - Random Numbers

Uploaded by

Document Information

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Lecture 07 - Random Numbers

Uploaded by

Copyright:

Available Formats

Pramod Parajuli

Simulation and Modeling, CS-331

Must have two essential properties;

Continuous random variable – Random variable X is

Each random number Ri is an independent sample

Expected value, E(R) = ∫ xdx f(x)

-Linear Congruent Generator

-Produces integers between 0 and m-1

Look for the period/cycle length

Maximum period is 64/4 = 16

Choose ‘m’ be one more than largest integer that

The general idea is to eliminate undesirable

 Frequency test (uniformity test)

 Counts how often numbers in a given range

Therefore, how much deviation should we allow from

 If χ2 is too large, we reject the null hypothesis (i.e. the

where Χ 2 k −1,1−α is upper critical point

Are we going to reject the hypothesis or not?

SN(x)= number of R1, R2, . . . , RN which are ≤ x

Example 7.6, page 267

 Runs up and down don't indicate any absolute

where α is the level of significance

Sample from book, page 274

As in chi-square test for the samples, the runs

If null hypothesis of independence is true, then

Example 7.10 – page 275

 Tests relationship of items m (lag) locations

 For large value of M, the distribution of the

 Do not reject the null hypothesis of

 If ρ im > 0 , the subsequence is said to exhibit

 The gap test is used to determine the

For digit 3, the first gap is 10 and second gap is 7.

 The probability of the digits not being 3 is 0.9

Example 7.13, page 282

Alternatively, the last result can be obtained as follows;

Normally used with chi-square test.

 Random numbers drawn from a distribution that

Inverse transformation method

For exponential distribution

In case of discrete data, if the result falls between

The general idea is, if a value falls in between two

Recall the Monte Carlo method

You might also like

Frequency test (uniformity test)

Counts how often numbers in a given range

If χ2 is too large, we reject the null hypothesis (i.e. the

Runs up and down don't indicate any absolute

Tests relationship of items m (lag) locations

For large value of M, the distribution of the

Do not reject the null hypothesis of

If ρ im > 0 , the subsequence is said to exhibit

The gap test is used to determine the

The probability of the digits not being 3 is 0.9

Random numbers drawn from a distribution that