Intro

A random number generator (RNG) generates a series of numbers which are pseudo random. There are many efficient algorithms (usually take less than a dozen of arithmetic/logical operations) of RNGs. Statistical tests need to be done to evaluate the randomness of a series of numbers \(\{R_i\}\) generated by an RNG. A simple visual test (similar to The Pickover Test) is to plot \((R_i, R_{i+k})\), for \(k = 1,2,3\). Here there is no apparent patterns in the scatter plot, indicating the RNG is good.

library(scales)
set.seed(10)   # set the RNG seed
n = 3000 # total number of values to be tested
r = runif(n)  # generate a series of n pseduo-random numbers from the R default RNG

## a visual test for RNGs
rng_test <- function(r){
  # r: a series of random numbers
  par(mar=c(5, 4, 4, 8), xpd=TRUE)
  plot(r[1:(n-1)], r[2:n], type = 'p', pch = 16, cex = 0.4, col = alpha(2, 0.5),
       xlab = bquote(R[i]), ylab = bquote(R[i+k]))
  points(r[1:(n-2)], r[3:n], pch = 16, cex = 0.4, col = alpha(3, 0.5))
  points(r[1:(n-3)], r[4:n], pch = 16, cex = 0.4, col = alpha(4, 0.5))
  legend("topright", inset=c(-0.3,0), title = "k",legend = c("1", "2", "3"), 
         col=c(2:4), pch= c(16,16,16))
}

rng_test(r)

Keep in mind that the plot would look different for RNGs of different types of distribution (uniform, exponential, Gaussian).

Uniform distribution of discrete numbers

In R, we can generate discrete random numbers in a uniform distribution by using the sample function. Note that the sampling needs to be done with replacement (i.e., “replace = T”), otherwise the random numbers would be correlated.

n = 10000
r_binary = sample(x = 2, size = n, replace = T)
hist(r_binary, freq = F, xlab = "r", breaks = 2, main = "Histogram of randomly sampled 1 and 2")

Exponential distribution

From an RNG that generates random numbers \(x\) from a uniform distribution of (0, 1), we can easily obtain random numbers \(y\) from an exponential distribution. (See the lecture presentation for the derivation)

\[ y = -ln(x) \tag{1}\]

n = 10000
r = runif(n)
y = -log(r) # random numbers from an exponential distribution
hist(r, freq = F, xlab = "r", main = "Histogram of random numbers from a uniform distribution")

hist(y, freq = F, xlab = "y", main = "Histogram of random numbers from an exponential distribution")

Gaussian distribution

To obtain random numbers from a Gaussian distribution, we use the Box-Muller method. We first sample two uniform random numbers from (-1, 1), \(x\) and \(y\), which satisfies \(0<R^2<1\), where \(R^2 = x^2 + y^2\). If not, another two numbers \(x\) and \(y\) are sampled.

\[\begin{equation} \begin{cases} r_1 = x\sqrt{-2ln(R^2)/R^2} \\ r_2 = y\sqrt{-2ln(R^2)/R^2} \tag{2} \end{cases} \end{equation}\]

# RNG for a Gaussian distribution, of course you don't need this. But it's interesting to know how it works.
ran_gaussian <- function(){
  repeat{
    x = 2*runif(2) -1 # generate two random numbers uniformly from [-1, 1]
    R2 = sum(x*x)
    if(R2 < 1.0 & abs(R2) > 10^-6 ) break
  }
  fac = sqrt(-2*log(R2)/R2)
  return(x*fac)
}

n = 10000
y = unlist(replicate(n,ran_gaussian()))
hist(y[1,], freq = F, xlab = "r_1", main = "Histogram of random numbers (set 1) from a Gaussian distribution")

hist(y[2,], freq = F, xlab = "r_2", main = "Histogram of random numbers (set 2) from a Gaussian distribution")