<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/combined-linear-congruential-generators-with-r/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://www.aaronschlegel.com/combined-linear-congruential-generators-r/"> R – Aaron Schlegel</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<div>Part of 3 in the series <a rel="nofollow" href="http://www.aaronschlegel.com/series/random-number-generation/" title="Random Number Generation" target="_blank">Random Number Generation</a></div>
<p>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:</p>
<img alt=" X_i \equiv \Bigg(\sum^k_{j=1} (-1)^{j-1} Y_{i,j} \Bigg) \space (\text{mod} \space (m_1 - 1)) " title=" X_i \equiv \Bigg(\sum^k_{j=1} (-1)^{j-1} Y_{i,j} \Bigg) \space (\text{mod} \space (m_1 - 1)) " src="https://s0.wp.com/latex.php?latex=+X_i+%5Cequiv+%5CBigg%28%5Csum%5Ek_%7Bj%3D1%7D+%28-1%29%5E%7Bj-1%7D+Y_%7Bi%2Cj%7D+%5CBigg%29+%5Cspace+%28%5Ctext%7Bmod%7D+%5Cspace+%28m_1+-+1%29%29+&bg=ffffff&fg=000&s=0" class="latex">
<p>Where <img alt="m_1" title="m_1" src="https://s0.wp.com/latex.php?latex=m_1&bg=ffffff&fg=000&s=0" class="latex"> is the modulus of the LCG, <img alt="Y_{i,j}" title="Y_{i,j}" src="https://s0.wp.com/latex.php?latex=Y_%7Bi%2Cj%7D&bg=ffffff&fg=000&s=0" class="latex"> is the <img alt="i" title="i" src="https://s0.wp.com/latex.php?latex=i&bg=ffffff&fg=000&s=0" class="latex">th input from the <img alt="j" title="j" src="https://s0.wp.com/latex.php?latex=j&bg=ffffff&fg=000&s=0" class="latex">th LCG and <img alt="X_i" title="X_i" src="https://s0.wp.com/latex.php?latex=X_i&bg=ffffff&fg=000&s=0" class="latex"> is the <img alt="i" title="i" src="https://s0.wp.com/latex.php?latex=i&bg=ffffff&fg=000&s=0" class="latex">th random generated value.</p>
<p>L’Ecuyer describes a combined linear generator that utilizes two LCGs in <em>Efficient and Portable Combined Random Number Generators</em> for 32-bit processors. To be precise, the congruential generators used are actually multiplicative since <img alt="c_1 = c_2 = 0" title="c_1 = c_2 = 0" src="https://s0.wp.com/latex.php?latex=c_1+%3D+c_2+%3D+0&bg=ffffff&fg=000&s=0" class="latex">. The parameters used for the MCGs are:</p>
<img alt=" a_1 = 40014 \qquad m_1 = 2147483563 \qquad a_2 = 40692 \qquad m_2 = 2147483399 " title=" a_1 = 40014 \qquad m_1 = 2147483563 \qquad a_2 = 40692 \qquad m_2 = 2147483399 " src="https://s0.wp.com/latex.php?latex=+a_1+%3D+40014+%5Cqquad+m_1+%3D+2147483563+%5Cqquad+a_2+%3D+40692+%5Cqquad+m_2+%3D+2147483399+&bg=ffffff&fg=000&s=0" class="latex">
<p>The combined linear congruential generator algorithm proposed by L’Ecuyer can be described with the following steps:</p>
<p>The two MCGs, <img alt="Y_{0,1}, \space Y_{0,2}" title="Y_{0,1}, \space Y_{0,2}" src="https://s0.wp.com/latex.php?latex=Y_%7B0%2C1%7D%2C+%5Cspace+Y_%7B0%2C2%7D&bg=ffffff&fg=000&s=0" class="latex">, are seeded. The seed values are recommended to be in the range <img alt="[1, m_1 - 1] \space \text{and} \space [1, m_2 - 1]" title="[1, m_1 - 1] \space \text{and} \space [1, m_2 - 1]" src="https://s0.wp.com/latex.php?latex=%5B1%2C+m_1+-+1%5D+%5Cspace+%5Ctext%7Band%7D+%5Cspace+%5B1%2C+m_2+-+1%5D&bg=ffffff&fg=000&s=0" class="latex">, respectively.</p>
<p>Next, the two MCGs are evaluated using the algorithm above:</p>
<img alt=" Y_{i+1,1} = a_1 \times Y_{i,1} (\text{mod} \space m_1) \qquad Y_{i+1,2} = a_1 \times Y_{i,2} (\text{mod} \space m_2) " title=" Y_{i+1,1} = a_1 \times Y_{i,1} (\text{mod} \space m_1) \qquad Y_{i+1,2} = a_1 \times Y_{i,2} (\text{mod} \space m_2) " src="https://s0.wp.com/latex.php?latex=+Y_%7Bi%2B1%2C1%7D+%3D+a_1+%5Ctimes+Y_%7Bi%2C1%7D+%28%5Ctext%7Bmod%7D+%5Cspace+m_1%29+%5Cqquad+Y_%7Bi%2B1%2C2%7D+%3D+a_1+%5Ctimes+Y_%7Bi%2C2%7D+%28%5Ctext%7Bmod%7D+%5Cspace+m_2%29+&bg=ffffff&fg=000&s=0" class="latex">
<p>With <img alt="Y_{i+1,1} \space \text{and} \space Y_{i+1,2}" title="Y_{i+1,1} \space \text{and} \space Y_{i+1,2}" src="https://s0.wp.com/latex.php?latex=Y_%7Bi%2B1%2C1%7D+%5Cspace+%5Ctext%7Band%7D+%5Cspace+Y_%7Bi%2B1%2C2%7D&bg=ffffff&fg=000&s=0" class="latex"> evaluated, find <img alt="X_{i+1}" title="X_{i+1}" src="https://s0.wp.com/latex.php?latex=X_%7Bi%2B1%7D&bg=ffffff&fg=000&s=0" class="latex"></p>
<img alt=" X_{i+1} = (Y_{i+1,1} - Y_{i+1,2}) \space \text{mod} \space m_1 - 1 " title=" X_{i+1} = (Y_{i+1,1} - Y_{i+1,2}) \space \text{mod} \space m_1 - 1 " src="https://s0.wp.com/latex.php?latex=+X_%7Bi%2B1%7D+%3D+%28Y_%7Bi%2B1%2C1%7D+-+Y_%7Bi%2B1%2C2%7D%29+%5Cspace+%5Ctext%7Bmod%7D+%5Cspace+m_1+-+1+&bg=ffffff&fg=000&s=0" class="latex">
<p>Finally, the random number to be output can be generated:</p>
<p><a rel="nofollow" href="http://www.aaronschlegel.com/?attachment_id=22650" target="_blank"><img alt="" src="https://i1.wp.com/www.aaronschlegel.com/wp-content/uploads/random_variate.png?resize=290%2C104" class="wp-image-22650" data-recalc-dims="1"></a></p>
<h6>Combined Linear Congruential Generators in R</h6>
<p>The following is an implementation of the CLCG algorithm proposed by L’Ecuyer for 32-bit processors. The function utilizes the <code>runif()</code> 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. <a rel="nofollow" href="https://www.johndcook.com/blog/2016/01/29/random-number-generator-seed-mistakes/" target="_blank">John Cook</a> 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.</p>
<pre>
clcg.32bit &lt;- function(n=10) {
  
  rng &lt;- vector(length = n)
  
  a1 &lt;- 40014
  m1 &lt;- 2147483563
  a2 &lt;- 40692
  m2 &lt;- 2147483399
  
  # Seed the two MCGs
  y1 &lt;- runif(1, 1, m1 - 1)
  y2 &lt;- runif(1, 1, m2 - 1)
  
  for (i in 1:n) {
    y1 &lt;- a1 * y1 %% m1
    y2 &lt;- a2 * y2 %% m2
    
    x &lt;- (y1 - y2) %% (m1 - 1)
    
    if (x > 0) {
      rng[i] &lt;- x / m1
    }
    else if (x &lt; 0) {
      rng[i] &lt;- (x / m1) + 1
    }
    else if (x == 0) {
      rng[i] &lt;- (m1 - 1) / m1
    }
  }
  
  return(rng)
}
</pre>
<p>Use our new function to generate ten random values.</p>
<pre>
clcg.32bit()
##  [1] 0.7075798 0.1233885 0.7122606 0.2641577 0.2006503 0.1901152 0.9532196
##  [8] 0.6168034 0.6464510 0.9276738
</pre>
<p>As before with our implementations of <a rel="nofollow" href="http://wp.me/p4aZEo-5SP" target="_blank">linear congruential generators</a> and <a rel="nofollow" href="http://wp.me/p4aZEo-5SV" target="_blank">multiplicative congruential generators</a>, we can visualize the apparent randomness of the generator function by plotting in three dimensions. To do this, generate three random vectors <img alt="x, y, z" title="x, y, z" src="https://s0.wp.com/latex.php?latex=x%2C+y%2C+z&bg=ffffff&fg=000&s=0" class="latex"> with the <code>clcg.32bit()</code> function and plot the points. The <a rel="nofollow" href="https://cran.r-project.org/web/packages/plot3D/index.html" target="_blank">plot3d package</a> is used to create the scatterplot, and the <a rel="nofollow" href="https://cran.r-project.org/web/packages/animation/" target="_blank">animation package</a> is used to animate each scatterplot as the length of the random vectors, <img alt="n" title="n" src="https://s0.wp.com/latex.php?latex=n&bg=ffffff&fg=000&s=0" class="latex">, increases.</p>
<pre>
library(plot3D)
library(animation)
</pre>
<pre>
n &lt;- c(3, 10, 20, 100, 500, 1000, 2000, 5000, 10000, 20000)
saveGIF({
  for (i in 1:length(n)) {
    x &lt;- clcg.32bit(n[i])
    y &lt;- clcg.32bit(n[i])
    z &lt;- 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')
</pre>
<p><a rel="nofollow" href="http://www.aaronschlegel.com/?attachment_id=22635" target="_blank"><img alt="three-dimensional representation of 32-bit combined linear congruential generator" src="https://i2.wp.com/www.aaronschlegel.com/wp-content/uploads/clcg_32bit.gif?resize=480%2C480" class="wp-image-22635" data-recalc-dims="1"></a></p>
<p>As <img alt="n" title="n" src="https://s0.wp.com/latex.php?latex=n&bg=ffffff&fg=000&s=0" class="latex"> increases, it appears the function generates suitably random values for simulation.</p>
<p>L’Ecuyer also gives a combined linear congruential generator for use with 16-bit processors. This CLCG uses three MCGs with the following parameters:</p>
<img alt=" a_1 = 157 \qquad m_1 = 32363 \qquad a_2 = 146 \qquad m_2 = 31727 \qquad a_3 = 142 \qquad m_3 = 31657 " title=" a_1 = 157 \qquad m_1 = 32363 \qquad a_2 = 146 \qquad m_2 = 31727 \qquad a_3 = 142 \qquad m_3 = 31657 " src="https://s0.wp.com/latex.php?latex=+a_1+%3D+157+%5Cqquad+m_1+%3D+32363+%5Cqquad+a_2+%3D+146+%5Cqquad+m_2+%3D+31727+%5Cqquad+a_3+%3D+142+%5Cqquad+m_3+%3D+31657+&bg=ffffff&fg=000&s=0" class="latex">
<p>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.</p>
<pre>
clcg.16bit &lt;- function(n=10) {
  rng &lt;- vector(length = n)
  
  a1 &lt;- 157
  m1 &lt;- 32363
  a2 &lt;- 146
  m2 &lt;- 31727
  a3 &lt;- 142
  m3 &lt;- 31657
  
  y1 &lt;- runif(1, 1, m1 - 1)
  y2 &lt;- runif(1, 1, m2 - 1)
  y3 &lt;- runif(1, 1, m3 - 1)
  
  for (i in 1:n) {
    y1 &lt;- a1 * y1 %% m1
    y2 &lt;- a2 * y2 %% m2
    y3 &lt;- a3 * y3 %% m3
    
    x &lt;- (y1 - y2 - y3) %% (m1 - 1)
    
    if (x > 0) {
      rng[i] &lt;- x / m1
    }
    else if (x &lt; 0) {
      rng[i] &lt;- (x / m1) + 1
    }
    else if (x == 0) {
      rng[i] &lt;- (m1 - 1) / m1
    }
  }
  
  return(rng)
}
</pre>
<p>Generate 10 random values using our newly created <code>clcg.16bit()</code> function.</p>
<pre>
clcg.16bit()
##  [1] 0.82370838 0.10198057 0.69001135 0.76613703 0.58562123 0.71786818
##  [7] 0.13731054 0.88403250 0.08715983 0.35320291
</pre>
<p>As before, let’s visualize the apparent randomness of the function by plotting three randomly generated vectors of increasing sizes in three dimensions.</p>
<pre>
n &lt;- c(3, 10, 20, 100, 500, 1000, 2000, 5000, 10000, 20000)
saveGIF({
  for (i in 1:length(n)) {
    x &lt;- clcg.16bit(n[i])
    y &lt;- clcg.16bit(n[i])
    z &lt;- 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')
</pre>
<p><a rel="nofollow" href="http://www.aaronschlegel.com/?attachment_id=22636" target="_blank"><img alt="three-dimensional representation of 16-bit combined linear congruential generator" src="https://i2.wp.com/www.aaronschlegel.com/wp-content/uploads/clcg_16bit.gif?resize=480%2C480" class="wp-image-22636" data-recalc-dims="1"></a></p>
<h6>References</h6>
<p>Combined Linear Congruential Generator. (2017, July 5). In Wikipedia, The Free Encyclopedia. From <a rel="nofollow" href="https://en.wikipedia.org/w/index.php?title=Combined_Linear_Congruential_Generator&oldid=789099445" target="_blank">https://en.wikipedia.org/w/index.php?title=Combined_Linear_Congruential_Generator&amp;oldid=789099445</a></p>
<p>Pierre L’Ecuyer (1988). Efficient and Portable Combined Random Number Generators. Communications of the ACM. 31: 742–749, 774. <a rel="nofollow" href="doi:10.1145/62959.62969" target="_blank">doi:10.1145/62959.62969</a></p>
<p>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</p>
<p>The post <a rel="nofollow" href="http://www.aaronschlegel.com/combined-linear-congruential-generators-r/" target="_blank">Combined Linear Congruential Generators with R</a> appeared first on <a rel="nofollow" href="http://www.aaronschlegel.com/" target="_blank">Aaron Schlegel</a>.</p>
	<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="http://www.aaronschlegel.com/combined-linear-congruential-generators-r/"> R – Aaron Schlegel</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a rel="nofollow" href="https://www.r-project.org/" title="The R Project for Statistical Computing">R</a> news and <a rel="nofollow" href="https://www.r-bloggers.com/search/tutorial" title="R tutorials">tutorials</a> on topics such as: <a rel="nofollow" href="https://www.r-bloggers.com/search/data%20science" title="Data science">Data science</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data" title="Big Data">Big Data, <a rel="nofollow" href="https://www.r-users.com/" title="R jobs">R jobs</a>, visualization (<a rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2" title="ggplot and ggplot2 tutorials">ggplot2</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/boxplot" title="Boxplots using lattice and ggplot2 tutorials">Boxplots</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/map" title="Maps and gis">maps</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/animation" title="Animation in R">animation</a>), programming (<a rel="nofollow" href="https://www.r-bloggers.com/search/RStudio" title="RStudio IDE for R">RStudio</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/sweave" title="Sweave and literate programming">Sweave</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX" title="LaTeX in R">LaTeX</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/SQL" title="SQL and databases">SQL</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/eclipse" title="Eclipse IDE for R">Eclipse</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/git" title="git and github, Version Control System">git</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/hadoop" title="Large data in R using Hadoop">hadoop</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping" title="Web Scraping of google, facebook, yahoo, twitter and more using R">Web Scraping</a>) statistics (<a rel="nofollow" href="https://www.r-bloggers.com/search/regression" title="Regressions and ANOVA analysis tutorials">regression</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/PCA" title="principal component analysis tutorial">PCA</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/time+series" title="Time series">time series</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/trading" title="finance trading">trading</a>) and more...
<img width="1" alt="" src="http://feeds.feedburner.com/~r/RBloggers/~4/3GT7eMaLojY" height="1"></a></div>