<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/how-to-sample-from-multidimensional-distributions-using-gibbs-sampling/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://appsilondatascience.com/blog/rstats/2017/10/09/gibbs-sampling.html"> Appsilon Data Science Blog</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p><strong>We will show how to perform multivariate random sampling using one of the Markov Chain Monte Carlo (MCMC) algorithms, called the Gibbs sampler.</strong></p>
<hr>
<p><small>Random sampling with rabbit on the bed plane</small></p>
<p><iframe width="450" allowfullscreen src="https://giphy.com/embed/49euH7z1UKBy0"></iframe></p>
<p><a rel="nofollow" href="https://giphy.com/gifs/rabbit-49euH7z1UKBy0" target="_blank">via GIPHY</a></p>
<hr>
<p><strong>To start, what are MCMC algorithms and what are they based on?</strong></p>
<p>Suppose we are interested in generating a random variable  with a distribution of , over .<br>
If we are not able to do this directly, we will be satisfied with generating a sequence of random variables , which in a sense tending to a distribution of . The MCMC approach explains how to do so:</p>
<p>Build a Markov chain , for , whose stationary distribution is .<br>
If the structure is correct, we should expect random variables to converge to .</p>
<p>In addition, we can expect that<br>
for function , occurs:<br>
</p>
<p>with probability equals to 1.</p>
<p>In essence, we want the above equality to occur for any arbitrary random variable .</p>
<hr>
<p><strong>One of the MCMC algorithms that guarantee the above properties is the so-called Gibbs sampler.</strong></p>
<p><strong>Gibbs Sampler – description of the algorithm.</strong></p>
<p>Assumptions:</p>
<ul>
<li> is defined on the product space </li>
<li>We are able to draw from the conditional distributions<br>
,<br>
where<br>
</li>
</ul>
<p>Algorithm steps:</p>
<ol>
<li>Select the initial values </li>
<li>For  repeat:<br>
<blockquote>
<p>For  sample  from distribution </p>
</blockquote>
</li>
<li>Repeat step 2 until the distribution of vector  stabilizes.</li>
<li>The subsequent iterations of point 2 will provide a randomization of .</li>
</ol>
<p>How do we understand point 3?<br>
We are most satisfied when average coordinates stabilize to some accuracy.</p>
<hr>
<p><strong>Intuition in two-dimensional case:</strong><br>
<img alt="Gibbs sampling on the plane." align="center" src="https://i2.wp.com/appsilondatascience.com/blog/assets/article_images/2017-10-02-gibbs-sampling/gibbs.png?w=100%25" data-recalc-dims="1"><br>
<small>Source: <a rel="nofollow" href="http://vcla.stat.ucla.edu/old/MCMC/MCMC_tutorial.htm" target="_blank">[3]</a></small></p>
<hr>
<p><strong>Gibbs sampling for randomization with a two-dimensional normal distribution.</strong></p>
<p>We will sample from the distribution of , where<br>
 and .</p>
<p>Knowing joint density of , it’s easy to show, that:<br>
<br>
</p>
<p>
R implementation:</p>
<figure>
<pre><code data-lang="r"><span>gibbs_normal_sampling</span><span> </span><span>&lt;-</span><span> </span><span>function</span><span>(</span><span>n_iteration</span><span>,</span><span> </span><span>init_point</span><span>,</span><span> </span><span>ro</span><span>)</span><span> </span><span>{</span><span>
</span><span># init point is some numeric vector of length equals to 2
</span><span> </span><span>theta_1</span><span> </span><span>&lt;-</span><span> </span><span>c</span><span>(</span><span>init_point</span><span>[</span><span>1</span><span>],</span><span> </span><span>numeric</span><span>(</span><span>n_iteration</span><span>))</span><span>
</span><span>theta_2</span><span> </span><span>&lt;-</span><span> </span><span>c</span><span>(</span><span>init_point</span><span>[</span><span>2</span><span>],</span><span> </span><span>numeric</span><span>(</span><span>n_iteration</span><span>))</span><span>
</span><span>for</span><span> </span><span>(</span><span>i</span><span> </span><span>in</span><span> </span><span>2</span><span>:</span><span>(</span><span>n_iteration</span><span>+1</span><span>))</span><span> </span><span>{</span><span>
</span><span>theta_1</span><span>[</span><span>i</span><span>]</span><span> </span><span>&lt;-</span><span> </span><span>rnorm</span><span>(</span><span>1</span><span>,</span><span> </span><span>ro</span><span> </span><span>*</span><span> </span><span>theta_2</span><span>[</span><span>i</span><span>-1</span><span>],</span><span> </span><span>sqrt</span><span>(</span><span>1</span><span> </span><span>-</span><span> </span><span>ro</span><span>^</span><span>2</span><span>))</span><span>
</span><span>theta_2</span><span>[</span><span>i</span><span>]</span><span> </span><span>&lt;-</span><span> </span><span>rnorm</span><span>(</span><span>1</span><span>,</span><span> </span><span>ro</span><span> </span><span>*</span><span> </span><span>theta_1</span><span>[</span><span>i</span><span>],</span><span> </span><span>sqrt</span><span>(</span><span>1</span><span> </span><span>-</span><span> </span><span>ro</span><span>^</span><span>2</span><span>))</span><span>
</span><span>}</span><span>
</span><span>list</span><span>(</span><span>theta_1</span><span>,</span><span> </span><span>theta_2</span><span>)</span><span>
</span><span>}</span></code></pre>
</figure>
<p>
Using the above function, let’s see how the 10 points were scored at :<br>
<img alt="Sampling ten points via gibbs sampler - animation." src="https://i0.wp.com/appsilondatascience.com/blog/assets/article_images/2017-10-02-gibbs-sampling/animation.gif?w=456" data-recalc-dims="1"></p>
<p>And for 10000 iterations:<br>
<img alt="Sampling from bivariate normal distribution." src="https://i1.wp.com/appsilondatascience.com/blog/assets/article_images/2017-10-02-gibbs-sampling/sampler.png?w=100%25" data-recalc-dims="1"></p>
<p>We leave you a comparison of how the stability of the parameters changes depending on the selected  parameter.</p>
<hr>
<p><strong>Let’s move on to use the Gibbs sampler to estimate the density parameters.</strong></p>
<p>We will show the use of the Gibbs sampler and bayesian statistics to estimate the mean parameters in the mix of normal distributions.</p>
<p>Assumptions (simplified case):</p>
<ul>
<li>iid. sample  comes from a mixture of normal distributions , where ,  i  are known.</li>
<li>For i=1,2  (a priori distributions) and  with  are independent.</li>
<li> – the classification vector (unobserved) from which Y is derived (when ,  is drawn from , when  then drawn from ).</li>
<li>

</li>
</ul>
<p>With above assumptions it can be shown that:</p>
<ul>
<li>

</li>
<li>

</li>
<li>

</li>
<li>

</li>
</ul>
<p>The form of the algorithm:</p>
<ol>
<li>Choose starting point for the mean </li>
<li>In the -th iteration do:<br>
<blockquote>
<ul>
<li>With the probability:<br> draw  for </li>
<li>Calculate:<br><br></li>
<li>Generate:<br><br></li>
</ul>
</blockquote>
</li>
<li>Perform step 2. until the distribution of vector  stabilizes.</li>
</ol>
<hr>
<p>How to do this in R?</p>
<p>At start let’s generate random sample from mixture of normals with parameters .</p>
<figure>
<pre><code data-lang="r"><span>set.seed</span><span>(</span><span>12345</span><span>)</span><span>
</span><span>mu_1</span><span> </span><span>&lt;-</span><span> </span><span>10</span><span>
</span><span>mu_2</span><span> </span><span>&lt;-</span><span> </span><span>2</span><span>
</span><span>sigma_1</span><span> </span><span>&lt;-</span><span> </span><span>1</span><span>
</span><span>sigma_2</span><span> </span><span>&lt;-</span><span> </span><span>2</span><span>
</span><span>pi_known</span><span> </span><span>&lt;-</span><span> </span><span>0.7</span><span>
</span><span>n</span><span> </span><span>&lt;-</span><span> </span><span>2000</span><span>
</span><span>pi_sampled</span><span> </span><span>&lt;-</span><span> </span><span>rbinom</span><span>(</span><span>n</span><span>,</span><span> </span><span>1</span><span>,</span><span> </span><span>pi_known</span><span>)</span><span>
</span><span>y_1</span><span> </span><span>&lt;-</span><span> </span><span>rnorm</span><span>(</span><span>n</span><span>,</span><span> </span><span>mu_1</span><span>,</span><span> </span><span>sigma_1</span><span>)</span><span>
</span><span>y_2</span><span> </span><span>&lt;-</span><span> </span><span>rnorm</span><span>(</span><span>n</span><span>,</span><span> </span><span>mu_2</span><span>,</span><span> </span><span>sigma_2</span><span>)</span><span>
</span><span>y</span><span> </span><span>&lt;-</span><span> </span><span>(</span><span>1</span><span> </span><span>-</span><span> </span><span>pi_sampled</span><span>)</span><span> </span><span>*</span><span> </span><span>y_1</span><span> </span><span>+</span><span> </span><span>pi_sampled</span><span> </span><span>*</span><span> </span><span>y_2</span></code></pre>
</figure>
<p>Take a quick look at a histogram of our data:</p>
<p><img alt="" src="https://i0.wp.com/appsilondatascience.com/blog/figs/2017-10-02-gibbs-sampling/unnamed-chunk-4-1.png?w=456" data-recalc-dims="1"></p>
<p>The following task is to estimate  and  from above model.</p>
<figure>
<pre><code data-lang="r"><span>mu_init</span><span> </span><span>&lt;-</span><span> </span><span>c</span><span>(</span><span>12</span><span>,</span><span> </span><span>0</span><span>)</span><span>
</span><span>n_iterations</span><span> </span><span>&lt;-</span><span> </span><span>300</span><span>
</span><span>delta</span><span> </span><span>&lt;-</span><span> </span><span>vector</span><span>(</span><span>&quot;list&quot;</span><span>,</span><span> </span><span>n_iterations</span><span>)</span><span>
</span><span>mu_1_vec</span><span> </span><span>&lt;-</span><span> </span><span>c</span><span>(</span><span>mu_init</span><span>[</span><span>1</span><span>],</span><span> </span><span>numeric</span><span>(</span><span>n_iterations</span><span>))</span><span>
</span><span>mu_2_vec</span><span> </span><span>&lt;-</span><span> </span><span>c</span><span>(</span><span>mu_init</span><span>[</span><span>2</span><span>],</span><span> </span><span>numeric</span><span>(</span><span>n_iterations</span><span>))</span><span>
</span><span>delta_probability</span><span> </span><span>&lt;-</span><span> </span><span>function</span><span>(</span><span>i</span><span>,</span><span> </span><span>y</span><span>,</span><span> </span><span>mu_1_vec</span><span>,</span><span> </span><span>mu_2_vec</span><span>,</span><span> </span><span>sigma_1</span><span>,</span><span> </span><span>sigma_2</span><span>,</span><span> </span><span>pi_known</span><span>)</span><span> </span><span>{</span><span>
</span><span>pi_known</span><span> </span><span>*</span><span> </span><span>dnorm</span><span>(</span><span>y</span><span>,</span><span> </span><span>mu_2_vec</span><span>[</span><span>i</span><span> </span><span>-</span><span> </span><span>1</span><span>],</span><span> </span><span>sigma_2</span><span>)</span><span> </span><span>/</span><span>
</span><span>((</span><span>1</span><span>-</span><span> </span><span>pi_known</span><span>)</span><span> </span><span>*</span><span> </span><span>dnorm</span><span>(</span><span>y</span><span>,</span><span> </span><span>mu_1_vec</span><span>[</span><span>i</span><span> </span><span>-</span><span> </span><span>1</span><span>],</span><span> </span><span>sigma_1</span><span>)</span><span> </span><span>+</span><span> </span><span>pi_known</span><span> </span><span>*</span><span> </span><span>dnorm</span><span>(</span><span>y</span><span>,</span><span> </span><span>mu_2_vec</span><span>[</span><span>i</span><span> </span><span>-</span><span> </span><span>1</span><span>],</span><span> </span><span>sigma_2</span><span>))</span><span>
</span><span>}</span><span>
</span><span>mu_1_mean</span><span> </span><span>&lt;-</span><span> </span><span>function</span><span>(</span><span>delta</span><span>,</span><span> </span><span>i</span><span>,</span><span> </span><span>y</span><span>)</span><span> </span><span>{</span><span>
</span><span>sum</span><span>((</span><span>1</span><span> </span><span>-</span><span> </span><span>delta</span><span>[[</span><span>i</span><span> </span><span>-</span><span> </span><span>1</span><span>]])</span><span> </span><span>*</span><span> </span><span>y</span><span>)</span><span> </span><span>/</span><span> </span><span>(</span><span>1</span><span> </span><span>+</span><span> </span><span>sum</span><span>(</span><span>1</span><span> </span><span>-</span><span> </span><span>delta</span><span>[[</span><span>i</span><span> </span><span>-</span><span> </span><span>1</span><span>]]))</span><span>
</span><span>}</span><span>
</span><span>mu_2_mean</span><span> </span><span>&lt;-</span><span> </span><span>function</span><span>(</span><span>delta</span><span>,</span><span> </span><span>i</span><span>,</span><span> </span><span>y</span><span>)</span><span> </span><span>{</span><span>
</span><span>sum</span><span>(</span><span>delta</span><span>[[</span><span>i</span><span> </span><span>-</span><span> </span><span>1</span><span>]]</span><span> </span><span>*</span><span> </span><span>y</span><span>)</span><span> </span><span>/</span><span> </span><span>(</span><span>1</span><span> </span><span>+</span><span> </span><span>sum</span><span>(</span><span>delta</span><span>[[</span><span>i</span><span> </span><span>-</span><span> </span><span>1</span><span>]]))</span><span>
</span><span>}</span><span>
</span><span>for</span><span> </span><span>(</span><span>j</span><span> </span><span>in</span><span> </span><span>2</span><span>:</span><span>(</span><span>n_iterations</span><span> </span><span>+</span><span> </span><span>1</span><span>))</span><span> </span><span>{</span><span>
</span><span>delta</span><span>[[</span><span>j</span><span> </span><span>-</span><span> </span><span>1</span><span>]]</span><span> </span><span>&lt;-</span><span> </span><span>y</span><span> </span><span>%>%</span><span> </span><span>map_int</span><span>(</span><span>
</span><span>~</span><span> </span><span>rbinom</span><span>(</span><span>1</span><span>,</span><span> </span><span>1</span><span>,</span><span> </span><span>delta_probability</span><span>(</span><span>j</span><span>,</span><span> </span><span>.</span><span>,</span><span> </span><span>mu_1_vec</span><span>,</span><span> </span><span>mu_2_vec</span><span>,</span><span> </span><span>sigma_1</span><span>,</span><span> </span><span>sigma_2</span><span>,</span><span> </span><span>pi_known</span><span>)))</span><span>
</span><span>mu_1_vec</span><span>[</span><span>j</span><span>]</span><span> </span><span>&lt;-</span><span> </span><span>rnorm</span><span>(</span><span>1</span><span>,</span><span> </span><span>mu_1_mean</span><span>(</span><span>delta</span><span>,</span><span> </span><span>j</span><span>,</span><span> </span><span>y</span><span>),</span><span> </span><span>sqrt</span><span>(</span><span>1</span><span> </span><span>/</span><span> </span><span>(</span><span>1</span><span> </span><span>+</span><span> </span><span>sum</span><span>(</span><span>1</span><span> </span><span>-</span><span> </span><span>delta</span><span>[[</span><span>j</span><span> </span><span>-</span><span> </span><span>1</span><span>]]))))</span><span>
</span><span>mu_2_vec</span><span>[</span><span>j</span><span>]</span><span> </span><span>&lt;-</span><span> </span><span>rnorm</span><span>(</span><span>1</span><span>,</span><span> </span><span>mu_2_mean</span><span>(</span><span>delta</span><span>,</span><span> </span><span>j</span><span>,</span><span> </span><span>y</span><span>),</span><span> </span><span>sqrt</span><span>(</span><span>1</span><span> </span><span>/</span><span> </span><span>(</span><span>1</span><span> </span><span>+</span><span> </span><span>sum</span><span>(</span><span>delta</span><span>[[</span><span>j</span><span> </span><span>-</span><span> </span><span>1</span><span>]]))))</span><span>
</span><span>}</span></code></pre>
</figure>
<p>
Let’s see the relation of sampled data to known one:</p>
<p>The following plot presents the mean of the  vector at each iteration:</p>
<p></p>
<p><img alt="" src="https://i1.wp.com/appsilondatascience.com/blog/figs/2017-10-02-gibbs-sampling/unnamed-chunk-6-1.png?w=456" data-recalc-dims="1"></p>
<p></p>
<p>Let’s check how mean of parameters  and  stabilize at used algorithm:</p>
<p><img alt="" src="https://i0.wp.com/appsilondatascience.com/blog/figs/2017-10-02-gibbs-sampling/unnamed-chunk-7-1.png?w=456" data-recalc-dims="1"></p>
<p>Note how little iteration was enough to stabilize the parameters.</p>
<p></p>
<p>Finally let’s see estimated density with our initial sample:</p>
<p><img alt="" src="https://i1.wp.com/appsilondatascience.com/blog/figs/2017-10-02-gibbs-sampling/unnamed-chunk-8-1.png?w=456" data-recalc-dims="1"></p>
<hr>
<p>To those concerned about the topic, refer to <a rel="nofollow" href="http://gauss.stat.su.se/rr/RR2006_1.pdf" target="_blank">[1]</a> where you can find a generalization of normal distribution mixture by extending a priori distributions to other parameters.</p>
<p>It is also worth to compare the above algorithm with its deterministic counterpart, Expectation Maximization (EM) algorithm see <a rel="nofollow" href="http://web.stanford.edu/~hastie/ElemStatLearn/" target="_blank">[2]</a>.</p>
<hr>
<p>Write your question and comments below. We’d love to hear what you think.</p>
<hr>
<p>Resources:</p>
<p>[1] http://gauss.stat.su.se/rr/RR2006_1.pdf</p>
<p>[2] http://web.stanford.edu/~hastie/ElemStatLearn/</p>
<p>[3] http://vcla.stat.ucla.edu/old/MCMC/MCMC_tutorial.htm</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://appsilondatascience.com/blog/rstats/2017/10/09/gibbs-sampling.html"> Appsilon Data Science Blog</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/vFKd1L45dB0" height="1"></a></div>