4.3 Generating Pseudo-Random Numbers

The literature on generating pseudo-random numbers is now extremely vast and it is not our purpose to review it, neither for you to learn how such algorithms work.

4.3.1 Generating Pseudo-Random Numbers in R

R has all the capabilities to generate such numbers. This can be done with the function runif, which takes one input: the number of observations to generate. So for instance:

set.seed(2021)
runif(10)
##  [1] 0.4512674 0.7837798 0.7096822 0.3817443 0.6363238 0.7013460 0.6404389
##  [8] 0.2666797 0.8154215 0.9829869

generates ten random numbers between zero and one. Notice that if we repeat the same code we get the same result since we fixed the so-called seed of the simulation.

set.seed(2021)
runif(10)
##  [1] 0.4512674 0.7837798 0.7096822 0.3817443 0.6363238 0.7013460 0.6404389
##  [8] 0.2666797 0.8154215 0.9829869

Conversely, if we were to simply run the code runif(10) we would get a different result.

runif(10)
##  [1] 0.02726706 0.83749040 0.60324073 0.56745337 0.82005281 0.25157128
##  [7] 0.50549403 0.86753810 0.95818157 0.54569770

4.3.2 Linear Congruential Method

Although we will use the functions already implemented in R, it is useful to at least introduce one of the most classical algorithms to simulate random numbers, called the linear congruential method. This produces a sequence of integers \(x_1,x_2,x_3\) between 0 and \(m-1\) using the recursion: \[ x_{i}=(ax_{i-1}+c)\mod m, \hspace{1cm} \mbox{for } i = 1,2,\dots \] Some comments:

  • \(\mod m\) is the remainder of the integer division by \(m\). For instance \(5 \mod 2\) is one and \(4\mod 2\) is zero.

  • therefore, the algorithm generates integers between 0 and \(m-1\).

  • there are three parameters that need to be chosen \(a, c\) and \(m\).

  • the value \(x_0\) is the seed of the algorithm.

Random numbers between zero and one can be derived by setting \[ u_i= x_i/m. \]

It can be shown that the method works well for specific choices of \(a\), \(c\) and \(m\), which we will not discuss here.

Let’s look at an implementation.

lcm <- function(N, x0, a, c, m){
   x <- rep(0,N)
   x[1] <- x0
   for (i in 2:N) x[i] <- (a*x[i-1]+c)%%m
   u <- x/m
   return(u)
}
lcm(N = 8, x0 = 4, a = 13, c = 0, m = 64)
## [1] 0.0625 0.8125 0.5625 0.3125 0.0625 0.8125 0.5625 0.3125

We can see that this specific choice of parameters is quite bad: it has cycle 4! After 4 numbers the sequence repeats itself and we surely would not like to use this in practice.

In general you should not worry of these issues, R does things properly for you!