Introduction to Statistical Computation

In this section we are interested in looking into the process of random number generation in computer systems. Though this may not seem too obvious, computers do not priduce trul random numbers. This is because, a computer by virtue of it’s very design is quite deterministic. I mean to say that the process by which a computer/process generates random numbers will be transigent on the underlying algorithm. So, if one can beat the algorith, one can always realise the generated numbers.

I remember once watching an online lecture from MITCourseware, with the professor giving a lecture on “Stochastic Processes”. He was narrating a story of his time in an investment firm in the 80’s when computers were only beggining to be utilized en masse in firms. A fellow colleague was playing around on one of them and was plotting random process as graphs. The professor narrating the story was telling the class that when he attempted the same (plotting random processes) on another system, seeing the same graph on both systems, he remarked that he thought he’d broken through the fifth dimension.

Of course, I mention this story only to emphasise that the number generation done by systems is deterministic and dependent on the underlying algorithm. As a result this is called a ’Pseudo Random Number Generator“. This will be important to keep in mind as we go through the rest of the textbook.

1. Linear Congruential Generator

We’ll first look into a very basic random number generator called the Linear Congruential Generator. This is now a primitive algortihm, but it serves as a good starting point for purely educational purposes. Understaing the basic notion of pseudo number generation, will at a later point enable us to encode complex(not necessarily difficult) methodologies. The basic algorithm is given below:

1.1 Algorithm

input:

m > 1(modulus)

a \(\in\) {1,2…m-1} (multiplier)

c \(\in\) {0,1…m-1} (increment)

X₀ \(\in\) {0,1…m-1} (Random seed)

output:

a sequence of random numbers X₁, X₂, X₃…:

1.for n = 1, 2, 3…do

2.X_n = (a**X_(n-1)+c) mod m*

3.output X_n

4.end for

LCG <- function(n,m,a,c,X0){
  X <- numeric(length=n)
  Xn <- X0
  for (i in 1:n){
    Xn <- (a*Xn + c) %% m
    X[i] <- Xn
  }
  return(X)
} 

Output

LCG(10,8,5,1,0) \(\newline\) [1] 1 6 7 4 5 2 3 0 1 6

As mentioned in the algorithm, a numeric vector is fed through a linear equation and the result is the random number we get as the output. This is awfully simplistic. Further reasons as to why this is ineffective and roundly dismissed for any practical purposes will be discusse din the next post. @fenner2012a says

References