Combined Linear Congruential Generators with R

This article is originally published at http://www.aaronschlegel.com

Combined linear congruential generators, as the name implies, are a type of PRNG (pseudorandom number generator) that combine two or more LCGs (linear congruential generators). The combination of two or more LCGs into one random number generator can result in a marked increase in the period length of the generator which makes them better suited for simulating more complex systems. The combined linear congruential generator algorithm is defined as:

Where m_1 is the modulus of the LCG, Y_{i,j} is the ith input from the jth LCG and X_i is the ith random generated value.

L’Ecuyer describes a combined linear generator that utilizes two LCGs in Efficient and Portable Combined Random Number Generators for 32-bit processors. To be precise, the congruential generators used are actually multiplicative since c_1 = c_2 = 0. The parameters used for the MCGs are:

The combined linear congruential generator algorithm proposed by L’Ecuyer can be described with the following steps:

The two MCGs, Y_{0,1}, \space Y_{0,2}, are seeded. The seed values are recommended to be in the range [1, m_1 - 1] \space \text{and} \space [1, m_2 - 1], respectively.

Next, the two MCGs are evaluated using the algorithm above:

With Y_{i+1,1} \space \text{and} \space Y_{i+1,2} evaluated, find X_{i+1}

Finally, the random number to be output can be generated:

Combined Linear Congruential Generators in R

The following is an implementation of the CLCG algorithm proposed by L’Ecuyer for 32-bit processors. The function utilizes the runif() function from base R to generate the seeds in the recommended ranges. Although there is nothing drastically wrong with seeding a generator using another randomly generated integer, it is not recommended for practical applications for numerous reasons. John Cook gives a good summary of the potential downfalls of seeding a generator with another random integer. However, this is more than fine for the proposes of demonstration.

clcg.32bit < - function(n=10) {
  
  rng <- vector(length = n)
  
  a1 <- 40014
  m1 <- 2147483563
  a2 <- 40692
  m2 <- 2147483399
  
  # Seed the two MCGs
  y1 <- runif(1, 1, m1 - 1)
  y2 <- runif(1, 1, m2 - 1)
  
  for (i in 1:n) {
    y1 <- a1 * y1 %% m1
    y2 <- a2 * y2 %% m2
    
    x <- (y1 - y2) %% (m1 - 1)
    
    if (x > 0) {
      rng[i] < - x / m1
    }
    else if (x < 0) {
      rng[i] <- (x / m1) + 1
    }
    else if (x == 0) {
      rng[i] <- (m1 - 1) / m1
    }
  }
  
  return(rng)
}

Use our new function to generate ten random values.

clcg.32bit()
##  [1] 0.7075798 0.1233885 0.7122606 0.2641577 0.2006503 0.1901152 0.9532196
##  [8] 0.6168034 0.6464510 0.9276738

As before with our implementations of linear congruential generators and multiplicative congruential generators, we can visualize the apparent randomness of the generator function by plotting in three dimensions. To do this, generate three random vectors x, y, z with the clcg.32bit() function and plot the points. The plot3d package is used to create the scatterplot, and the animation package is used to animate each scatterplot as the length of the random vectors, n, increases.

library(plot3D)
library(animation)

n < - c(3, 10, 20, 100, 500, 1000, 2000, 5000, 10000, 20000)
saveGIF({
  for (i in 1:length(n)) {
    x <- clcg.32bit(n[i])
    y <- clcg.32bit(n[i])
    z <- clcg.32bit(n[i])
    
    scatter3D(x, y, z, colvar = NULL, pch=20, cex=0.3, alpha=0.75, theta=20, main = paste('n = ', n[i]))
  }
}, movie.name = 'clcg_32bit.gif')

As n increases, it appears the function generates suitably random values for simulation.

L’Ecuyer also gives a combined linear congruential generator for use with 16-bit processors. This CLCG uses three MCGs with the following parameters:

The algorithm proceeds in the same way as the 32-bit generator. As such, we can implement the algorithm in much the same way as the 32-bit function. Please note that if we were genuinely interested in implementing CLCGs for use with different byte sizes, we would need to take more care than we are giving here for managing the memory constraints and other restrictions.

clcg.16bit < - function(n=10) {
  rng <- vector(length = n)
  
  a1 <- 157
  m1 <- 32363
  a2 <- 146
  m2 <- 31727
  a3 <- 142
  m3 <- 31657
  
  y1 <- runif(1, 1, m1 - 1)
  y2 <- runif(1, 1, m2 - 1)
  y3 <- runif(1, 1, m3 - 1)
  
  for (i in 1:n) {
    y1 <- a1 * y1 %% m1
    y2 <- a2 * y2 %% m2
    y3 <- a3 * y3 %% m3
    
    x <- (y1 - y2 - y3) %% (m1 - 1)
    
    if (x > 0) {
      rng[i] < - x / m1
    }
    else if (x < 0) {
      rng[i] <- (x / m1) + 1
    }
    else if (x == 0) {
      rng[i] <- (m1 - 1) / m1
    }
  }
  
  return(rng)
}

Generate 10 random values using our newly created clcg.16bit() function.

clcg.16bit()
##  [1] 0.82370838 0.10198057 0.69001135 0.76613703 0.58562123 0.71786818
##  [7] 0.13731054 0.88403250 0.08715983 0.35320291

As before, let’s visualize the apparent randomness of the function by plotting three randomly generated vectors of increasing sizes in three dimensions.

n < - c(3, 10, 20, 100, 500, 1000, 2000, 5000, 10000, 20000)
saveGIF({
  for (i in 1:length(n)) {
    x <- clcg.16bit(n[i])
    y <- clcg.16bit(n[i])
    z <- clcg.16bit(n[i])
    
    scatter3D(x, y, z, colvar = NULL, pch=20, cex=0.3, alpha=0.75, theta=20, main = paste('n = ', n[i]))
  }
}, movie.name = 'clcg_16bit.gif')

References

Combined Linear Congruential Generator. (2017, July 5). In Wikipedia, The Free Encyclopedia. From https://en.wikipedia.org/w/index.php?title=Combined_Linear_Congruential_Generator&oldid=789099445

Pierre L’Ecuyer (1988). Efficient and Portable Combined Random Number Generators. Communications of the ACM. 31: 742–749, 774. doi:10.1145/62959.62969

Pierre L’Ecuyer, (1999) Good Parameters and Implementations for Combined Multiple Recursive Random Number Generators. Operations Research 47(1):159-164. doi.org/10.1287/opre.47.1.159

The post Combined Linear Congruential Generators with R appeared first on Aaron Schlegel.

Thanks for visiting r-craft.org
This article is originally published at http://www.aaronschlegel.com
Please visit source website for post related comments.