<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/multiplicative-congruential-generators-in-r/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://www.aaronschlegel.com/multiplicative-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 2 of 2 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>Multiplicative congruential generators, also known as Lehmer random number generators, is a type of linear congruential generator for generating pseudorandom numbers in <img alt="U(0, 1)" title="U(0, 1)" src="https://s0.wp.com/latex.php?latex=U%280%2C+1%29&bg=ffffff&fg=000&s=0" class="latex">. The multiplicative congruential generator, often abbreviated as MLCG or MCG, is defined as a recurrence relation similar to the LCG with <img alt="c = 0" title="c = 0" src="https://s0.wp.com/latex.php?latex=c+%3D+0&bg=ffffff&fg=000&s=0" class="latex">.</p>
<img alt="\large{X_{i+1} = aX_i \space \text{mod} \space m} " title="\large{X_{i+1} = aX_i \space \text{mod} \space m} " src="https://s0.wp.com/latex.php?latex=%5Clarge%7BX_%7Bi%2B1%7D+%3D+aX_i+%5Cspace+%5Ctext%7Bmod%7D+%5Cspace+m%7D+&bg=ffffff&fg=000&s=0" class="latex">
<p>Unlike the LCG, the parameters <img alt="a" title="a" src="https://s0.wp.com/latex.php?latex=a&bg=ffffff&fg=000&s=0" class="latex"> and <img alt="m" title="m" src="https://s0.wp.com/latex.php?latex=m&bg=ffffff&fg=000&s=0" class="latex"> for multiplicative congruential generators are more restricted and the initial seed <img alt="X_0" title="X_0" src="https://s0.wp.com/latex.php?latex=X_0&bg=ffffff&fg=000&s=0" class="latex"> must be relatively prime to the modulus <img alt="m" title="m" src="https://s0.wp.com/latex.php?latex=m&bg=ffffff&fg=000&s=0" class="latex"> (the greatest common divisor between <img alt="X_0" title="X_0" src="https://s0.wp.com/latex.php?latex=X_0&bg=ffffff&fg=000&s=0" class="latex"> and <img alt="m" title="m" src="https://s0.wp.com/latex.php?latex=m&bg=ffffff&fg=000&s=0" class="latex"> is <img alt="0" title="0" src="https://s0.wp.com/latex.php?latex=0&bg=ffffff&fg=000&s=0" class="latex">). The current parameters in common use are <img alt="m = 2^{31} - 1 = 2,147,483,647 \text{and} a = 7^5 = 16,807" title="m = 2^{31} - 1 = 2,147,483,647 \text{and} a = 7^5 = 16,807" src="https://s0.wp.com/latex.php?latex=m+%3D+2%5E%7B31%7D+-+1+%3D+2%2C147%2C483%2C647+%5Ctext%7Band%7D+a+%3D+7%5E5+%3D+16%2C807&bg=ffffff&fg=000&s=0" class="latex">. However, in a correspondence from the <em>Communications of the ACM</em>, Park, Miller and Stockmeyer changed the value of the parameter <img alt="a" title="a" src="https://s0.wp.com/latex.php?latex=a&bg=ffffff&fg=000&s=0" class="latex">, stating:</p>
<blockquote><p>
The minimal standard Lehmer generator we advocated had a modulus of m = 2^31 – 1 and a multiplier of a = 16807. Relative to this particular choice of multiplier, we wrote “… if this paper were to be written again in a few years it is quite possible that we would advocate a different multiplier ….” We are now prepared to do so. That is, we now advocate a = 48271 and, indeed, have done so “officially” since July 1990. This new advocacy is consistent with the discussion on page 1198 of [10]. There is nothing wrong with 16807; we now believe, however, that 48271 is a little better (with q = 44488, r = 3399).
</p>
</blockquote>
<h6>Multiplicative Congruential Generators with Schrage’s Method</h6>
<p>When using a large prime modulus <img alt="m" title="m" src="https://s0.wp.com/latex.php?latex=m&bg=ffffff&fg=000&s=0" class="latex"> such as <img alt="2^{31} - 1" title="2^{31} - 1" src="https://s0.wp.com/latex.php?latex=2%5E%7B31%7D+-+1&bg=ffffff&fg=000&s=0" class="latex">, the multiplicative congruential generator can overflow. Schrage’s method was invented to overcome the possibility of overflow. We can check the parameters in use satisfy this condition:</p>
<pre>
a &lt;- 48271
m &lt;- 2 ** 31 - 1

a * (m %% a) &lt; m
## [1] TRUE
</pre>
<p>Schrage’s method restates the modulus <img alt="m" title="m" src="https://s0.wp.com/latex.php?latex=m&bg=ffffff&fg=000&s=0" class="latex"> as a decomposition <img alt="m = aq + r" title="m = aq + r" src="https://s0.wp.com/latex.php?latex=m+%3D+aq+%2B+r&bg=ffffff&fg=000&s=0" class="latex"> where <img alt="r = m \space \text{mod} \space a" title="r = m \space \text{mod} \space a" src="https://s0.wp.com/latex.php?latex=r+%3D+m+%5Cspace+%5Ctext%7Bmod%7D+%5Cspace+a&bg=ffffff&fg=000&s=0" class="latex"> and <img alt="q = m / a" title="q = m / a" src="https://s0.wp.com/latex.php?latex=q+%3D+m+%2F+a&bg=ffffff&fg=000&s=0" class="latex">.</p>
<img alt="ax \space \text{mod} \space m = \begin{cases} a(x \space \text{mod} \space q) - r\frac{x}{q} &amp; \text{if} \space x \space \text{is} \geq 0 \\ a(x \space \text{mod} \space q) - r\frac{x}{q} + m &amp; \text{if} \space x \space \text{is} \leq 0 \end{cases} " title="ax \space \text{mod} \space m = \begin{cases} a(x \space \text{mod} \space q) - r\frac{x}{q} &amp; \text{if} \space x \space \text{is} \geq 0 \\ a(x \space \text{mod} \space q) - r\frac{x}{q} + m &amp; \text{if} \space x \space \text{is} \leq 0 \end{cases} " src="https://s0.wp.com/latex.php?latex=ax+%5Cspace+%5Ctext%7Bmod%7D+%5Cspace+m+%3D+%5Cbegin%7Bcases%7D+a%28x+%5Cspace+%5Ctext%7Bmod%7D+%5Cspace+q%29+-+r%5Cfrac%7Bx%7D%7Bq%7D+%26+%5Ctext%7Bif%7D+%5Cspace+x+%5Cspace+%5Ctext%7Bis%7D+%5Cgeq+0+%5C%5C+a%28x+%5Cspace+%5Ctext%7Bmod%7D+%5Cspace+q%29+-+r%5Cfrac%7Bx%7D%7Bq%7D+%2B+m+%26+%5Ctext%7Bif%7D+%5Cspace+x+%5Cspace+%5Ctext%7Bis%7D+%5Cleq+0+%5Cend%7Bcases%7D+&bg=ffffff&fg=000&s=0" class="latex">
<h6>Multiplicative Congruential Generator in R</h6>
<p>We can implement a Lehmer random number generator in R using the parameters mentioned earlier.</p>
<pre>
lehmer.rng &lt;- function(n=10) {
  
  rng &lt;- vector(length = n)
  
  m &lt;- 2147483647
  a &lt;- 48271
  q &lt;- 44488
  r &lt;- 3399
  
  # Set the seed using the current system time in microseconds. 
  # The initial seed value must be coprime to the modulus m, 
  # which we are not really concerning ourselves with for this example.
  d &lt;- as.numeric(Sys.time())
  
  for (i in 1:n) {
    h &lt;- d / q
    l &lt;- d %% q
    t &lt;- a * l - r * h
    if (t &lt; 0) {
      d &lt;- t
    }
    else {
      d &lt;- t + m
    }
    
    rng[i] &lt;- d / m
  }
  return(rng)
}
</pre>
<pre>
# Print the first 10 randomly generated numbers
lehmer.rng()
##  [1] 0.68635675 0.12657390 0.84869106 0.16614698 0.08108171 0.89533896
##  [7] 0.90708773 0.03195725 0.60847522 0.70736551
</pre>
<p>Plotting our multiplicative congruential generator in three dimensions allows us to visualize the apparent ‘randomness’ of the generator. As before, we 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 our Lehmer RNG 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;- lehmer.rng(n[i])
    y &lt;- lehmer.rng(n[i])
    z &lt;- lehmer.rng(n[i])
    
    scatter3D(x, y, z, colvar = NULL, pch=20, cex = 0.5, theta=20, main = paste('n = ', n[i]))
  }
}, movie.name = 'lehmer.gif')
</pre>
<p><a rel="nofollow" href="http://www.aaronschlegel.com/?attachment_id=22626" target="_blank"><img alt="Multiplicative Congruential Generators at several values of n" src="https://i1.wp.com/www.aaronschlegel.com/wp-content/uploads/lehmer.gif?resize=480%2C480" class="wp-image-22626" data-recalc-dims="1"></a></p>
<p>The generator appears to be generating suitably random numbers demonstrated by the increasing swarm of points as <img alt="n" title="n" src="https://s0.wp.com/latex.php?latex=n&bg=ffffff&fg=000&s=0" class="latex"> increases.</p>
<h6>References</h6>
<p>Anne Gille-Genest (March 1, 2012). Implementation of the Pseudo-Random Number Generators and the Low Discrepancy Sequences.</p>
<p>Saucier, R. (2000). Computer Generation of Statistical Distributions (1st ed.). Aberdeen, MD. Army Research Lab.</p>
<p>Stephen K. Park; Keith W. Miller; Paul K. Stockmeyer (1988). “Technical Correspondence”. Communications of the ACM. 36 (7): 105–110.</p>
<p>The post <a rel="nofollow" href="http://www.aaronschlegel.com/multiplicative-congruential-generators-r/" target="_blank">Multiplicative Congruential Generators in 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/multiplicative-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/q3dxHyNQsxE" height="1"></a></div>