2. Random number Generators

2.1. Random Number Generation

Earlier we talked about generating a random uniform (0, 1) number. So, how this can be done and why it is needed?

Any simulation requires random numbers to obtain random observations from specific probability distributions. For example, an observation from an exponential distribution that can be viewed as a realisation of a random service time in a queueing system (e.g., the cafeteria Wishbone).

  1. tables - e.g., RAND 1955

http://www.rand.org/publications/classics/randomdigits

11164 36318 75061 37674 26320 75100 10431 20418 19228 91792

21215 91791 76831 58678 87054 31687 93205 43685 19732 08468

10438 44482 66558 37649 08882 90870 12462 41810 01806 02977

36792 26236 33266 66583 60881 97395 20461 36742 02852 50564

73944 04773 12032 51414 82384 38370 00249 80709 72605 67497

49563 12872 14063 93104 78483 72717 68714 18048 25005 04151

64208 48237 41701 73117 33242 42314 83049 21933 92813 04763

51486 72875 38605 29341 80749 80151 33835 52602 79147 08868

99756 26360 64516 17971 48478 09610 04638 17141 09227 10606

71325 55217 13015 72907 00431 45117 33827 92873 02953 85474

65285 97198 12138 53010 94601 15838 16805 61004 43516 17020

17264 57327 38224 29301 31381 38109 34976 65692 98566 29550

95639 99754 31199 92558 68368 04985 51092 37780 40261 14479

61555 76404 86210 11808 12841 45147 97438 60022 12645 62000

78137 98768 04689 87130 79225 08153 84967 64539 79493 74917 \[ \]

  1. computer generated

The first method for random number generation was Von Neuman’s mid-square method:

Square previous number and extract the middle digits. For example, assume 10 digit-numbers and the previous was

\[5772156649 \rightarrow (5772156649)^2 \rightarrow 33317792380594909291\]
\(7923805949\) is the new number

\[ \]

A random number generator is an algorithm to produce sequences of numbers that follow specific probability distribution and possess the appearance of randomness. The method is acceptable if produces numbers that are:

  1. uniformly distributed

  2. statistically independent

  3. reproducible

  4. non-repeating for any desired length


We will use the following terminology:

  1. Random number - means a random observation from some (discrete or continuous) uniform distribution.

  2. Random observation - means an observation from another (not the uniform) probability distribution, e.g., exponential \[ \]

\[ \] Usually an integer random number (IRN) is produced and, if it is needed, converted to uniform random number (URN). Let \(k\) be an integer random number \(\in [0, n]\), then URN \(r \in (0, 1)\) can be produced as follows:

\begin{align*} URN &= \left\lbrace \begin{array}{ccc} & \frac{k}{n} \sim U(0,1)& \\ & \frac {k + \frac{1}{2}}{n + 1} \sim U(0,1) & \text{if $n$ is small} \\ \end{array} \right. \\ \end{align*}

\(n\) is small if \(m < 30\).

In reality, the discrete uniform integers \([0, n]\) are not random. They are called pseudo-random numbers, because they are

  1. predictable and

  2. reproducible.

2.2. Congruential Methods for IRN Generation

NOTE

Given two positive numbers, a (the dividend) and m (the divisor), a modulo m (abbreviated as a mod m) is the remainder of the division of a by m. For example, 7 mod 3 = 1; 12 mod 3 = 0; 23 mod 3 = 2

end note

2.2.1. Mixed congruential method

Generates sequence of random numbers over the range \(0,1,...,m-1\) \begin{align*} x_{n+1} &= \left( ax_{n} + c\right)(mod(m)) \end{align*}

\(a\), \(c\) and \(m\) are positive integers \(a<m\) and \(c<m\).

The initial random number is called seed \(( x_{0} )\) \begin{align*} \Rightarrow& x_{n+1} \; \text{is the remainder of } \; \frac{ax_{n} + c}{m} \\ \Rightarrow& x_{n+1} \in \left[ 0,1,...,m-1 \right] \end{align*}

which is the full cycle length.

Example: Let \(m = 8\), \(a = 5\), \(c = 7\), \(x_{0} = 4\). \begin{align*} &\text{Integer RN} &\text{Uniform RN} \\ x_{1} &=>\frac{(5\times4 + 7)}{8}=>3 & \frac{3+\frac{1}{2}}{8}=0.4375 \\ x_{2} &=> \frac{(5\times3 + 7)}{8}=>6 & \frac{6+\frac{1}{2}}{8}=0.8125 \\ x_{3} &=> \frac{(5\times6 + 7)}{8}=>5 & \frac{5+\frac{1}{2}}{8}=0.6875 \\ x_{4} &=> \frac{(5\times5 + 7)}{8}=>0 & \frac{0+\frac{1}{2}}{8}=0.0625 \\ x_{5} &=> \frac{(5\times0 + 7)}{8}=>7 & \frac{7+\frac{1}{2}}{8}=0.9375 \\ x_{6} &=> \frac{(5\times7 + 7)}{8}=>2 & \frac{2+\frac{1}{2}}{8}=0.3125 \\ x_{7} &=> \frac{(5\times2 + 7)}{8}=>1 & \frac{1+\frac{1}{2}}{8}=0.1875 \\ x_{8} &=> \frac{(5\times1 + 7)}{8}=>4 & \Rightarrow \text{back to } x_{0} \\ \end{align*}

Each of these URN is the middle point of one of the eight equal-sized intervals \((0, 0.125)\),\((0.125, 0.250)\),…\((0.875, 1.0)\), i.e., it is a rough approximation for real \(U(0,1)\).

Example: Let \(m = 8\), \(a = 5\), \(c = 7\), \(x_{0} = 0\). \begin{align*} x_{1} &=> \frac{(5\times0 + 7)}{8}=>7 & x_{6} &=> \frac{(5\times3 + 7)}{8}=>6 \\ x_{2} &=> \frac{(5\times7 + 7)}{8}=>2 & x_{7} &=> \frac{(5\times6 + 7)}{8}=>5 \\ x_{3} &=> \frac{(5\times2 + 7)}{8}=>1 & x_{8} &=> \frac{(5\times5 + 7)}{8}=>0 \\ x_{4} &=> \frac{(5\times1 + 7)}{8}=>4 & x_{9} &=> \frac{(5\times0 + 7)}{8}=>7 \\ x_{5} &=> \frac{(5\times4 + 7)}{8}=>3 & \\ \end{align*}

2.2.2. Multiplicative Method

\[c = 0\]

2.2.3. Additive Method

\[a = 1\] and replaces \(c\) by some random number preceeding \(x_n\) in the sequence, for example \(x_{n-1}.\) \[x_{n+1} = x_{n} + \mbox{some preceding } x_{i}\mbox{'s}.\] It needs more than one seed!

Example Fibonacci’s numbers/sequence

By definition, the first two numbers in the Fibonacci sequence are 1 and 1, or 0 and 1, depending on the chosen starting point of the sequence, and each subsequent number is the sum of the previous two. Consider initial seeds 1 and 1. Then

\[ 1,\;1,\;2,\;3,\;5,\;8,\;13,\;21,\;34,\;55,\;89,\;144,\; \ldots\;\]

is the Fibonacci numbers/sequence.

The R code to produce the Fibonacci sequence with initial seeds 1 and 1 is below

fn  <- rep(0, 15)
fn[1] <- 1
fn[2] <- 1
for (i in 3:15){
     fn[i] <- fn[i-1] + fn[i-2]
     }
print(fn)
##  [1]   1   1   2   3   5   8  13  21  34  55  89 144 233 377 610
  1. For a binary computer with a word size of \(b\) bits
    • the usual choice for \(m\) is \(m = 2^{b}\). This is the total number of non-negative integers that can be expressed within the capacity of the word size.
    • For this \(m\), we can ensure that each possible number occurs exactly once before any number is repeated by selecting \(a\) to be \(a = 1, 5, 9, 13, \ldots ..\)
    • and \(c\) to be \(c = 1, 3, 5, 7, \ldots..\).
  2. For a decimal computer with a word size of \(d\) digits
    • \(m = 10^{d}\)
    • \(a = 1, 21, 41, 61,\ldots\) and
    • \(c = 1, 3, 7, 9, 11, 13, 17,\ldots\) (all odds except multiples of \(5\))

2.2.4. Examples

Write an R code to produce pseudo-random numbers using a mixed congruential generator.

randnumber <- function(a,c,m,x0=0){
            x1 <- (a*x0 +c)%%m
            return(x1)}
x1 <- randnumber(5,7,8,x0=0)
print(x1)
## [1] 7

Next, generate the full cycle of the generator starting from any initial seed (which is a number between \(0\) and \((m-1)\)).

x0 <- 3 # choice of the seed
# xp - x previous
# xn - x next
xp <- x0
all <- c(x0)
for(i in 2:8){
   xn <- randnumber(5,7,8,all[i-1])
   all <- c(all, xn)
}
print(all)
## [1] 3 6 5 0 7 2 1 4