<p><div><div><div><div data-width="100" data-show-faces="false" data-send="true" data-layout="button_count" data-height="21" data-href="https://www.r-bloggers.com/simulating-from-the-bivariate-normal-distribution-in-r/"></div></div><div><a data-url="https://www.r-bloggers.com/simulating-from-the-bivariate-normal-distribution-in-r/" data-counturl="https://www.r-bloggers.com/simulating-from-the-bivariate-normal-distribution-in-r/" data-count="horizontal" data-text="Simulating from the Bivariate Normal Distribution in R" class="twitter-share-button" data-via="rbloggers" href="https://twitter.com/share"></a></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://blog.revolutionanalytics.com/2016/08/simulating-form-the-bivariate-normal-distribution-in-r-1.html"> Revolutions</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<div>
<p>by Joseph Rickert</p>
<p>My guess is that a good many statistics students first encounter the <a target="_blank" href="http://athenasc.com/Bivariate-Normal.pdf">bivariate Normal distribution</a> as one or two hastily covered pages in an introductory text book, and then don't think much about it again until someone asks them to generate two random variables with a given correlation structure. Fortunately for R users, a little searching on the internet will turn up  several nice tutorials with R code explaining various aspects of the bivariate Normal. For this post, I have gathered together a few examples and tweaked the code a little to make comparisons easier.</p>
<p>Here are five different ways to simulate random samples bivariate Normal distribution with a given mean and covariance matrix.</p>
<p>To set up for the simulations this first block of code defines N, the number of random samples to simulate, the means of the random variables, and and the covariance matrix. It also provides a small function for drawing confidence ellipses on the simulated data.</p>
<p><span>library(mixtools)  #for ellipse</span></p>
<p><span>N &lt;- 200 # Number of random samples</span><br><span>set.seed(123)</span><br><span># Target parameters for univariate normal distributions</span><br><span>rho &lt;- -0.6</span><br><span>mu1 &lt;- 1; s1 &lt;- 2</span><br><span>mu2 &lt;- 1; s2 &lt;- 8</span></p>
<p><span># Parameters for bivariate normal distribution</span><br><span>mu &lt;- c(mu1,mu2) # Mean </span><br><span>sigma &lt;- matrix(c(s1^2, s1*s2*rho, s1*s2*rho, s2^2),<br>           2) # Covariance matrix</span></p>
<p><span># Function to draw ellipse for bivariate normal data</span><br><span>ellipse_bvn &lt;- function(bvn, alpha){</span><br><span>  Xbar &lt;- apply(bvn,2,mean)</span><br><span>  S &lt;- cov(bvn)</span><br><span>  ellipse(Xbar, S, alpha = alpha, col=&quot;red&quot;)</span><br><span>}</span></p>
<p>The first method, the way to go if you just want to get on with it, is to use the <span>mvrnorm()</span> function from the MASS package. </p>
<p><span>library(MASS)</span><br><span>bvn1 &lt;- mvrnorm(N, mu = mu, Sigma = sigma ) # from MASS package<br>colnames(bvn1) &lt;- c(&quot;bvn1_X1&quot;,&quot;bvn1_X2&quot;)</span></p>
<p>It takes so little code to do the simulation it might be possible to tweet in a homework assignment.</p>
<p>A look at the source code for <span>mvrnorm()</span> shows that it uses eignevectors to generate the random samples. The documentation for the function states that this method was selected because it is stabler than the alternative of using a <a target="_blank" href="https://en.wikipedia.org/wiki/Cholesky_decomposition">Cholesky decomposition</a> which might be faster. </p>
<p>For the second method, let's go ahead and directly generate generate bivariate Normal random variates with the Cholesky decomposition. Remember that the Cholesky decomposition of sigma (a positive definite matrix) yields a matrix M such that M times its transpose gives sigma back again. Multiplying M by a matrix of standard random Normal variates and adding the desired mean gives a matrix of the desired random samples. A <a target="_blank" href="https://www2.stat.duke.edu/courses/Spring12/sta104.1/Lectures/Lec22.pdf">lecture</a> from Colin Rundel covers some of the theory.</p>
<p><span>M &lt;- t(chol(sigma))</span><br><span># M %*% t(M)</span><br><span>Z &lt;- matrix(rnorm(2*N),2,N) # 2 rows, N/2 columns</span><br><span>bvn2 &lt;- t(M %*% Z) + matrix(rep(mu,N), byrow=TRUE,ncol=2)<br>colnames(bvn2) &lt;- c(&quot;bvn2_X1&quot;,&quot;bvn2_X2&quot;)<br></span></p>
<p>For the third method we make use of a special property of the bivariate normal that is discussed in almost all of those elementary textbooks.  If X<sub>1</sub> and X<sub>2</sub> are two jointly distributed random variables, then the conditional distribution of X<sub>2</sub> given X<sub>1</sub> is itself normal with: mean = <span>m</span><sub>2</sub> + <span>r</span>(<span>s</span><sub>2</sub>/<span>s</span><sub>1</sub>)(X<sub>1</sub> – m<sub>1</sub>) and variance = (1 – <span>r</span><sup><sub>2</sub></sup>)<span>s</span><sup>2</sup>X<sub>2</sub>.</p>
<p>Hence, a sample from a bivariate Normal distribution can be simulated by first simulating a point from the marginal distribution of one of the random variables and then simulating from the second random variable conditioned on the first. A brief proof of the underlying theorem is available <a target="_blank" href="http://athenasc.com/Bivariate-Normal.pdf">here</a>.</p>
<p><span>rbvn&lt;-function (n, m1, s1, m2, s2, rho) </span><br><span>  {</span><br><span>     X1 &lt;- rnorm(n, mu1, s1)</span><br><span>     X2 &lt;- rnorm(n, mu2 + (s2/s1) * rho * <br>           (X1 – mu1), sqrt((1 – rho^2)*s2^2))</span><br><span>     cbind(X1, X2)</span><br><span>  }</span></p>
<p><span>bvn3 &lt;- rbvn(N,mu1,s1,mu2,s2,rho)<br>colnames(bvn3) &lt;- c(&quot;bvn3_X1&quot;,&quot;bvn3_X2&quot;)<br></span></p>
<p>The fourth method, my favorite, comes from Professor Darren Wiliinson's <a target="_blank" href="http://www.mas.ncl.ac.uk/~ndjw1/teaching/sim/gibbs/gibbs.html">Gibbs Sampler tutorial</a>. This is a very nice idea; using the familiar bivariate Normal distribution to illustrate the basics of the <a target="_blank" href="https://www.psych.umn.edu/faculty/waller/classes/Bayes/Readings/Casella%20and%20George%201992.pdf">Gibbs Sampling Algorithm</a>. Note that this looks very much like the previous method, except that now we are alternately sampling from the full conditional distributions.</p>
<p><span>gibbs&lt;-function (n, mu1, s1, mu2, s2, rho) </span><br><span>{</span><br><span>  mat &lt;- matrix(ncol = 2, nrow = n)</span><br><span>  x &lt;- 0</span><br><span>  y &lt;- 0</span><br><span>  mat[1, ] &lt;- c(x, y)</span><br><span>  for (i in 2:n) {</span><br><span>    x &lt;- rnorm(1, mu1 + <br>          (s1/s2) * rho * (y – mu2), sqrt((1 – rho^2)*s1^2))</span><br><span>    y &lt;- rnorm(1, mu2 + <br>          (s2/s1) * rho * (x – mu1), sqrt((1 – rho^2)*s2^2))</span><br><span>    mat[i, ] &lt;- c(x, y)</span><br><span>  }</span><br><span>  mat</span><br><span>}</span><br><span>bvn4 &lt;- gibbs(N,mu1,s1,mu2,s2,rho)<br>colnames(bvn4) &lt;- c(&quot;bvn4_X1&quot;,&quot;bvn4_X2&quot;)<br></span></p>
<p>The fifth and final way uses the <span>rmvnorm()</span> function form the <a target="_blank" href="https://mran.revolutionanalytics.com/package/mvtnorm/">mvtnorm package</a> with the s<a target="_blank" href="https://en.wikipedia.org/wiki/Singular_value_decomposition">ingular value decomposition</a> method selected. The functions in this package are overkill for what we are doing here, but mvtnorm is probably the package you would want to use if you are calculating probabilities from high dimensional multivariate distributions. It implements numerical methods for carefully calculating the high dimensional integrals involved that are based on some papers by Professor <a target="_blank" href="http://www.sci.wsu.edu/math/faculty/genz/homepage">Alan Genz</a> dating from the early '90s. These methods are briefly explained in the package <a target="_blank" href="https://mran.revolutionanalytics.com/web/packages/mvtnorm/vignettes/MVT_Rnews.pdf">vignette</a>. </p>
<p><span>library (mvtnorm)<br>bvn5 &lt;- mvtnorm::rmvnorm(N,mu,sigma, method=&quot;svd&quot;)</span><span><br>colnames(bvn5) &lt;- c(&quot;bvn5_X1&quot;,&quot;bvn5_X2&quot;)<br></span></p>
<p>Note that I have used the :: operator here to make sure that R uses the <span>rmvnorm()</span> function from the mvtnorm package. There is also a <span>rmvnorm()</span> function in the mixtools package that I used to get the ellipse function. Loading the packages in the wrong order could lead to the rookie mistake of having the function you want inadvertently overwritten.</p>
<p>Next, we plot the results of drawing just 100 random samples for each method. This allows us to see how the algorithms spread data over the sample space as they are just getting started.</p>
<p><span>bvn &lt;- list(bvn1,bvn2,bvn3,bvn4,bvn5)</span></p>
<p><span>par(mfrow=c(3,2))</span><br><span>plot(bvn1, xlab=&quot;X1&quot;,ylab=&quot;X2&quot;,main= &quot;All Samples&quot;)</span><br><span>for(i in 2:5){</span><br><span>  points(bvn[[i]],col=i)</span><br><span>}</span><br><span>for(i in 1:5){</span><br><span>  item &lt;- paste(&quot;bvn&quot;,i,sep=&quot;&quot;)</span><br><span>  plot(bvn[[i]],xlab=&quot;X1&quot;,ylab=&quot;X2&quot;,main=item, col=i)</span><br><span>  ellipse_bvn(bvn[[i]],.5)</span><br><span>  ellipse_bvn(bvn[[i]],.05)</span><br><span>}</span><br><span>par(mfrow=c(1,1))</span></p>
<p><a target="_blank" href="http://revolution-computing.typepad.com/.a/6a010534b1db25970b01bb09243521970d-pi"><img title="Bvn_final" alt="Bvn_final" src="http://revolution-computing.typepad.com/.a/6a010534b1db25970b01bb09243521970d-800wi" border="0"></a></p>
<p>The first plot shows all 500 random samples color coded by the method with which they were generated. The remaining plots show the samples generated by each method. In each of these plots the ellipses mark the 0.5 and 0.95 probability regions, i.e. the area within the ellipses should contain 50% and 95% of the points respectively. Note that bvn4 which uses the Gibbs sampling algorithm looks like all of the rest. In most use cases for the Gibbs it takes the algorithm some time to converge to the target distribution. In our case, we start out with a pretty good guess.</p>
<p>Finally, a word about accuracy: nice coverage of the sample space is not sufficient to produce accurate results. A little experimentation will show that, for all of the methods outlined above, regularly achieving a sample covariance matrix that is close to the target, sigma, requires something on the order of 10,000 samples as is Illustrated below.</p>
<p> <span>> sigma</span><br><span>      [,1] [,2]</span><br><span>[1,] 4.0 -9.6</span><br><span>[2,] -9.6 64.0</span></p>
<p><span>for(i in 1:5){</span><br><span>  print(round(cov(bvn[[i]]),1))</span><br><span>}</span></p>
<p><span>            bvn1_X1 bvn1_X2</span><br><span>bvn1_X1     4.0    -9.5</span><br><span>bvn1_X2    -9.5    63.8</span></p>
<p><span>           bvn2_X1 bvn2_X2</span><br><span>bvn2_X1     3.9    -9.5</span><br><span>bvn2_X2    -9.5    64.5</span></p>
<p><span>          bvn3_X1 bvn3_X2</span><br><span>bvn3_X1     4.1    -9.8</span><br><span>bvn3_X2    -9.8    63.7</span></p>
<p><span>          bvn4_X1 bvn4_X2</span><br><span>bvn4_X1     4.0    -9.7</span><br><span>bvn4_X2    -9.7    64.6</span></p>
<p><span>         bvn5_X1 bvn5_X2</span><br><span>bvn5_X1     4.0    -9.6</span><br><span>bvn5_X2    -9.6    65.3</span></p>
<p>Many people coming to R for the first time find it disconcerting to realize that there are several ways to do some fundamental calculation in R. My take is that rather than being a point of frustration, having multiple options indicates that richness of the R language. A close look at the package documentation will often show that yet another method to do something is a response to some subtle need that was not previously addressed. Enjoy the diversity!</p>
</div>
<p><div><div><div><div data-width="55" data-show-faces="false" data-send="true" data-layout="box_count" data-height="62" data-href="https://www.r-bloggers.com/simulating-from-the-bivariate-normal-distribution-in-r/"></div></div><div><a data-url="https://www.r-bloggers.com/simulating-from-the-bivariate-normal-distribution-in-r/" data-counturl="https://www.r-bloggers.com/simulating-from-the-bivariate-normal-distribution-in-r/" data-count="vertical" data-text="Simulating from the Bivariate Normal Distribution in R" class="twitter-share-button" data-via="rbloggers" href="https://twitter.com/share"></a></div><div></div></div><div></div></div>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="http://blog.revolutionanalytics.com/2016/08/simulating-form-the-bivariate-normal-distribution-in-r-1.html"> Revolutions</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 title="The R Project for Statistical Computing" rel="nofollow" href="https://www.r-project.org/">R</a> news and <a title="R tutorials" rel="nofollow" href="https://www.r-bloggers.com/?s=tutorial">tutorials</a> on topics such as: <a title="Data science" rel="nofollow" href="https://www.r-bloggers.com/?s=data%20science">Data science</a>, <a title="Big Data" rel="nofollow" href="https://www.r-bloggers.com/?s=Big%20Data">Big Data, <a title="R jobs" rel="nofollow" href="https://www.r-users.com/">R jobs</a>, visualization (<a title="ggplot and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/?s=ggplot2">ggplot2</a>, <a title="Boxplots using lattice and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/?s=boxplot">Boxplots</a>, <a title="Maps and gis" rel="nofollow" href="https://www.r-bloggers.com/?s=map">maps</a>, <a title="Animation in R" rel="nofollow" href="https://www.r-bloggers.com/?s=animation">animation</a>), programming (<a title="RStudio IDE for R" rel="nofollow" href="https://www.r-bloggers.com/?s=RStudio">RStudio</a>, <a title="Sweave and literate programming" rel="nofollow" href="https://www.r-bloggers.com/?s=sweave">Sweave</a>, <a title="LaTeX in R" rel="nofollow" href="https://www.r-bloggers.com/?s=LaTeX">LaTeX</a>, <a title="SQL and databases" rel="nofollow" href="https://www.r-bloggers.com/?s=SQL">SQL</a>, <a title="Eclipse IDE for R" rel="nofollow" href="https://www.r-bloggers.com/?s=eclipse">Eclipse</a>, <a title="git and github, Version Control System" rel="nofollow" href="https://www.r-bloggers.com/?s=git">git</a>, <a title="Large data in R using Hadoop" rel="nofollow" href="https://www.r-bloggers.com/?s=hadoop">hadoop</a>, <a title="Web Scraping of google, facebook, yahoo, twitter and more using R" rel="nofollow" href="https://www.r-bloggers.com/?s=Web+Scraping">Web Scraping</a>) statistics (<a title="Regressions and ANOVA analysis tutorials" rel="nofollow" href="https://www.r-bloggers.com/?s=regression">regression</a>, <a title="principal component analysis tutorial" rel="nofollow" href="https://www.r-bloggers.com/?s=PCA">PCA</a>, <a title="Time series" rel="nofollow" href="https://www.r-bloggers.com/?s=time+series">time series</a>, <a title="finance trading" rel="nofollow" href="https://www.r-bloggers.com/?s=trading">trading</a>) and more...
<img height="1" alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/GrdRRnmYzCA"></a></div>