<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/better-block-sampling-in-mcmc-with-the-automated-factor-slice-sampler/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://r-nimble.org/better-block-sampling-in-mcmc-with-the-automated-factor-slice-sampler"> R – NIMBLE</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>

<p><br>
</p>
<p>One nice feature of NIMBLE’s MCMC system is that a user can easily write new samplers from R, combine them with NIMBLE’s samplers, and have them automatically compiled to C++ via the NIMBLE compiler. We’ve observed that block sampling using a simple adaptive multivariate random walk Metropolis-Hastings sampler doesn’t always work well in practice, so we decided to implement the Automated Factor Slice sampler (AFSS) of <code><a></a><a target="_blank" rel="nofollow" href="https://r-nimble.org/better-block-sampling-in-mcmc-with-the-automated-factor-slice-sampler#bib-tibbits14">Tibbits, Groendyke, Haran, and Liechty (2014)</a></code> and see how it does on a (somewhat artificial) example with severe posterior correlation problems.</p>
<p>Roughly speaking, the AFSS works by conducting univariate slice sampling in directions determined by the eigenvectors of the marginal posterior covariance matrix for blocks of parameters in a model. So far, we’ve found the AFSS often outperforms random walk block sampling. To compare performance, we look at MCMC efficiency, which we define for each parameter as effective sample size (ESS) divided by computation time. We define overall MCMC efficiency as the minimum MCMC efficiency of all the parameters, because one needs all parameters to be well mixed.</p>
<p>We’ll demonstrate the performance of the AFSS on the correlated state space model described in <code><a></a><a target="_blank" rel="nofollow" href="https://r-nimble.org/better-block-sampling-in-mcmc-with-the-automated-factor-slice-sampler#bib-turek17">Turek, de<br>
Valpine, Paciorek, Anderson-Bergman, and others (2017)</a></code>.</p>
<div>
<h1>Model Creation</h1>
<p>Assume <img title="x_{i}" alt="x_{i}" class="latex" src="http://s0.wp.com/latex.php?latex=x_%7Bi%7D&bg=ffffff&%23038;fg=000000&%23038;s=0"> is the latent state and <img title="y_{i}" alt="y_{i}" class="latex" src="http://s0.wp.com/latex.php?latex=y_%7Bi%7D&bg=ffffff&%23038;fg=000000&%23038;s=0"> is the observation at time <img title="i" alt="i" class="latex" src="http://s0.wp.com/latex.php?latex=i&bg=ffffff&%23038;fg=000000&%23038;s=0"> for <img title="i=1,ldots,100" alt="i=1,ldots,100" class="latex" src="http://s0.wp.com/latex.php?latex=i%3D1%2C%5Cldots%2C100&bg=ffffff&%23038;fg=000000&%23038;s=0">. We define the state space model as</p>
<p><img title=" x_{i} sim N(a cdot x_{i-1} + b, sigma_{PN}) " alt=" x_{i} sim N(a cdot x_{i-1} + b, sigma_{PN}) " class="latex" src="http://s0.wp.com/latex.php?latex=+x_%7Bi%7D+%5Csim+N%28a+%5Ccdot+x_%7Bi-1%7D+%2B+b%2C+%5Csigma_%7BPN%7D%29+&bg=ffffff&%23038;fg=000000&%23038;s=0"><br>
<img title=" y_{i} sim N(x_{i}, sigma_{OE}) " alt=" y_{i} sim N(x_{i}, sigma_{OE}) " class="latex" src="http://s0.wp.com/latex.php?latex=+y_%7Bi%7D+%5Csim+N%28x_%7Bi%7D%2C+%5Csigma_%7BOE%7D%29+&bg=ffffff&%23038;fg=000000&%23038;s=0"></p>
<p>for <img title="i = 2, ldots, 100" alt="i = 2, ldots, 100" class="latex" src="http://s0.wp.com/latex.php?latex=i+%3D+2%2C+%5Cldots%2C+100&bg=ffffff&%23038;fg=000000&%23038;s=0">, with initial states</p>
<p><img title=" x_{1} sim N(frac{b}{1-a}, frac{sigma_{PN}}{sqrt{1-a^2}})" alt=" x_{1} sim N(frac{b}{1-a}, frac{sigma_{PN}}{sqrt{1-a^2}})" class="latex" src="http://s0.wp.com/latex.php?latex=+x_%7B1%7D+%5Csim+N%28%5Cfrac%7Bb%7D%7B1-a%7D%2C+%5Cfrac%7B%5Csigma_%7BPN%7D%7D%7B%5Csqrt%7B1-a%5E2%7D%7D%29&bg=ffffff&%23038;fg=000000&%23038;s=0"><br>
<img title=" y_{1} sim N(x_{1}, sigma_{OE})" alt=" y_{1} sim N(x_{1}, sigma_{OE})" class="latex" src="http://s0.wp.com/latex.php?latex=+y_%7B1%7D+%5Csim+N%28x_%7B1%7D%2C+%5Csigma_%7BOE%7D%29&bg=ffffff&%23038;fg=000000&%23038;s=0"></p>
<p>and prior distributions</p>
<p><img title=" a sim Unif(-0.999, 0.999)" alt=" a sim Unif(-0.999, 0.999)" class="latex" src="http://s0.wp.com/latex.php?latex=+a+%5Csim+Unif%28-0.999%2C+0.999%29&bg=ffffff&%23038;fg=000000&%23038;s=0"><br>
<img title=" b sim N(0, 1000) " alt=" b sim N(0, 1000) " class="latex" src="http://s0.wp.com/latex.php?latex=+b+%5Csim+N%280%2C+1000%29+&bg=ffffff&%23038;fg=000000&%23038;s=0"><br>
<img title=" sigma_{PN} sim Unif(0, 1) " alt=" sigma_{PN} sim Unif(0, 1) " class="latex" src="http://s0.wp.com/latex.php?latex=+%5Csigma_%7BPN%7D+%5Csim+Unif%280%2C+1%29+&bg=ffffff&%23038;fg=000000&%23038;s=0"><br>
<img title=" sigma_{OE} sim Unif(0, 1) " alt=" sigma_{OE} sim Unif(0, 1) " class="latex" src="http://s0.wp.com/latex.php?latex=+%5Csigma_%7BOE%7D+%5Csim+Unif%280%2C+1%29+&bg=ffffff&%23038;fg=000000&%23038;s=0"></p>
<p>where <img title="N(mu, sigma)" alt="N(mu, sigma)" class="latex" src="http://s0.wp.com/latex.php?latex=N%28%5Cmu%2C+%5Csigma%29&bg=ffffff&%23038;fg=000000&%23038;s=0"> denotes a normal distribution with mean <img title="mu" alt="mu" class="latex" src="http://s0.wp.com/latex.php?latex=%5Cmu&bg=ffffff&%23038;fg=000000&%23038;s=0"> and standard deviation <img title="sigma" alt="sigma" class="latex" src="http://s0.wp.com/latex.php?latex=%5Csigma&bg=ffffff&%23038;fg=000000&%23038;s=0">.</p>
<p>A file named <code>model_SSMcorrelated.RData</code> with the BUGS model code, data, constants, and initial values for our model can be downloaded <a target="_blank" rel="nofollow" href="https://github.com/danielturek/automated-blocking-examples/blob/master/data/model_SSMcorrelated.RData">here</a>.</p>
<div>
<div>
<div>
<pre><span>## load the nimble library and set seed</span>
<span>library</span><span>(</span><span>'nimble'</span><span>)</span>
<span>set.seed</span><span>(</span><span>1</span><span>)</span>
<span>load</span><span>(</span><span>'model_SSMcorrelated.RData'</span><span>)</span>
<span>## build and compile the model</span>
<span>stateSpaceModel</span> <span>&lt;-</span> <span>nimbleModel</span><span>(</span><span>code</span> <span>= code,</span>
                              <span>data</span> <span>= data,</span>
                              <span>constants</span> <span>= constants,</span>
                              <span>inits</span> <span>= inits,</span>
                              <span>check</span> <span>=</span> <span>FALSE</span><span>)</span>

<span>C_stateSpaceModel</span> <span>&lt;-</span> <span>compileNimble</span><span>(stateSpaceModel)</span>
</pre>
</div>
</div>
</div>
<div>
<h1>Comparing two MCMC Samplers</h1>
<p>We next compare the performance of two MCMC samplers on the state space model described above. The first sampler we consider is NIMBLE’s <code>RW_block</code> sampler, a Metropolis-Hastings sampler with a multivariate normal proposal distribution. This sampler has an adaptive routine that modifies the proposal covariance to look like the empirical covariance of the posterior samples of the parameters. However, as we shall see below, this proposal covariance adaptation does not lead to efficient sampling for our state space model.</p>
<p>We first build and compile the MCMC algorithm.</p>
<div>
<div>
<div>
<pre><span>RW_mcmcConfig</span> <span>&lt;-</span> <span>configureMCMC</span><span>(stateSpaceModel)</span>
<span>RW_mcmcConfig</span><span>$</span><span>removeSamplers</span><span>(</span><span>c</span><span>(</span><span>'a'</span><span>,</span> <span>'b'</span><span>,</span> <span>'sigOE'</span><span>,</span> <span>'sigPN'</span><span>))</span>
<span>RW_mcmcConfig</span><span>$</span><span>addSampler</span><span>(</span><span>target</span> <span>=</span> <span>c</span><span>(</span><span>'a'</span><span>,</span> <span>'b'</span><span>,</span> <span>'sigOE'</span><span>,</span> <span>'sigPN'</span><span>),</span> <span>type</span> <span>=</span> <span>'RW_block'</span><span>)</span>
<span>RW_mcmc</span> <span>&lt;-</span> <span>buildMCMC</span><span>(RW_mcmcConfig)</span>
<span>C_RW_mcmc</span> <span>&lt;-</span> <span>compileNimble</span><span>(RW_mcmc,</span> <span>project</span> <span>= stateSpaceModel)</span>
</pre>
</div>
</div>
</div>
<p>We next run the compiled MCMC algorithm for 10,000 iterations, recording the overall MCMC efficiency from the posterior output. The overall efficiency here is defined as <img title="min(frac{ESS}{T})" alt="min(frac{ESS}{T})" class="latex" src="http://s0.wp.com/latex.php?latex=min%28%5Cfrac%7BESS%7D%7BT%7D%29&bg=ffffff&%23038;fg=000000&%23038;s=0">, where ESS denotes the effective sample size, and <img title="T" alt="T" class="latex" src="http://s0.wp.com/latex.php?latex=T&bg=ffffff&%23038;fg=000000&%23038;s=0"> the total run-time of the sampling algorithm. The minimum is taken over all parameters that were sampled. We repeat this process 5 times to get a very rough idea of the average minimum efficiency for this combination of model and sampler.</p>
<div>
<div>
<div>
<pre><span>RW_minEfficiency</span> <span>&lt;-</span> <span>numeric</span><span>(</span><span>5</span><span>)</span>
<span>for</span><span>(i</span> <span>in</span> <span>1</span><span>:</span><span>5</span><span>){</span>
  <span>runTime</span> <span>&lt;-</span> <span>system.time</span><span>(C_RW_mcmc</span><span>$</span><span>run</span><span>(</span><span>50000</span><span>,</span> <span>progressBar</span> <span>=</span> <span>FALSE</span><span>))[</span><span>'elapsed'</span><span>]</span>
  <span>RW_mcmcOutput</span> <span>&lt;-</span> <span>as.mcmc</span><span>(</span><span>as.matrix</span><span>(C_RW_mcmc</span><span>$</span><span>mvSamples))</span>
  <span>RW_minEfficiency[i]</span> <span>&lt;-</span> <span>min</span><span>(</span><span>effectiveSize</span><span>(RW_mcmcOutput)</span><span>/</span><span>runTime)</span>
<span>}</span>
<span>summary</span><span>(RW_minEfficiency)</span>
</pre>
</div>
<div>
<pre>##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.3323  0.4800  0.5505  0.7567  0.7341  1.6869
</pre>
</div>
</div>
</div>
<p>Examining a trace plot of the output below, we see that the $a$ and $b$ parameters are mixing especially poorly.</p>
<div>
<div>
<div>
<pre><span>plot</span><span>(RW_mcmcOutput,</span> <span>density</span> <span>=</span> <span>FALSE</span><span>)</span>
</pre>
</div>
</div>
<div><img title="plot of chunk plot-mcmc" alt="plot of chunk plot-mcmc" src="https://i0.wp.com/r-nimble.org/wp-content/uploads/2017/05/plot-mcmc-1.png?w=456&ssl=1" data-recalc-dims="1"></div>
</div>
<p>Plotting the posterior samples of <img title="a" alt="a" class="latex" src="http://s0.wp.com/latex.php?latex=a&bg=ffffff&%23038;fg=000000&%23038;s=0"> against those of <img title="b" alt="b" class="latex" src="http://s0.wp.com/latex.php?latex=b&bg=ffffff&%23038;fg=000000&%23038;s=0"> reveals a strong negative correlation. This presents a problem for the Metropolis-Hastings sampler — we have found that adaptive algorithms used to tune the proposal covariance are often slow to reach a covariance that performs well for blocks of strongly correlated parameters.</p>
<div>
<div>
<div>
<pre><span>plot.default</span><span>(RW_mcmcOutput[,</span><span>'a'</span><span>], RW_mcmcOutput[,</span><span>'b'</span><span>])</span>
</pre>
</div>
</div>
<div><img title="plot of chunk plot-corr" alt="plot of chunk plot-corr" src="https://i1.wp.com/r-nimble.org/wp-content/uploads/2017/05/plot-corr-1.png?w=456&ssl=1" data-recalc-dims="1"></div>
<div>
<div>
<pre><span>cor</span><span>(RW_mcmcOutput[,</span><span>'a'</span><span>], RW_mcmcOutput[,</span><span>'b'</span><span>])</span>
</pre>
</div>
<div>
<pre>## [1] -0.9201277
</pre>
</div>
</div>
</div>
<p>In such situations with strong posterior correlation, we’ve found the AFSS to often run much more efficiently, so we next build and compile an MCMC algorithm using the AFSS sampler. Our hope is that the AFSS sampler will be better able to to produce efficient samples in the face of high posterior correlation.</p>
<div>
<div>
<div>
<pre><span>AFSS_mcmcConfig</span> <span>&lt;-</span> <span>configureMCMC</span><span>(stateSpaceModel)</span>
<span>AFSS_mcmcConfig</span><span>$</span><span>removeSamplers</span><span>(</span><span>c</span><span>(</span><span>'a'</span><span>,</span> <span>'b'</span><span>,</span> <span>'sigOE'</span><span>,</span> <span>'sigPN'</span><span>))</span>
<span>AFSS_mcmcConfig</span><span>$</span><span>addSampler</span><span>(</span><span>target</span> <span>=</span> <span>c</span><span>(</span><span>'a'</span><span>,</span> <span>'b'</span><span>,</span> <span>'sigOE'</span><span>,</span> <span>'sigPN'</span><span>),</span> <span>type</span> <span>=</span> <span>'AF_slice'</span><span>)</span>
<span>AFSS_mcmc</span><span>&lt;-</span> <span>buildMCMC</span><span>(AFSS_mcmcConfig)</span>
<span>C_AFSS_mcmc</span> <span>&lt;-</span> <span>compileNimble</span><span>(AFSS_mcmc,</span> <span>project</span> <span>= stateSpaceModel,</span> <span>resetFunctions</span> <span>=</span> <span>TRUE</span><span>)</span>
</pre>
</div>
</div>
</div>
<p>We again run the AFSS MCMC algorithm 5 times, each with 10,000 MCMC iterations.</p>
<div>
<div>
<div>
<pre><span>AFSS_minEfficiency</span> <span>&lt;-</span> <span>numeric</span><span>(</span><span>5</span><span>)</span>
<span>for</span><span>(i</span> <span>in</span> <span>1</span><span>:</span><span>5</span><span>){</span>
  <span>runTime</span> <span>&lt;-</span> <span>system.time</span><span>(C_AFSS_mcmc</span><span>$</span><span>run</span><span>(</span><span>50000</span><span>,</span> <span>progressBar</span> <span>=</span> <span>FALSE</span><span>))[</span><span>'elapsed'</span><span>]</span>
  <span>AFSS_mcmcOutput</span> <span>&lt;-</span> <span>as.mcmc</span><span>(</span><span>as.matrix</span><span>(C_AFSS_mcmc</span><span>$</span><span>mvSamples))</span>
  <span>AFSS_minEfficiency[i]</span> <span>&lt;-</span> <span>min</span><span>(</span><span>effectiveSize</span><span>(AFSS_mcmcOutput)</span><span>/</span><span>runTime)</span>
<span>}</span>
<span>summary</span><span>(AFSS_minEfficiency)</span>
</pre>
</div>
<div>
<pre>##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   9.467   9.686  10.549  10.889  10.724  14.020
</pre>
</div>
</div>
</div>
<p>Note that the minimum overall efficiency of the AFSS sampler is approximately <code>28</code> times that of the <code>RW_block</code> sampler. Additionally, trace plots from the output of the AFSS sampler show that the <img title="a" alt="a" class="latex" src="http://s0.wp.com/latex.php?latex=a&bg=ffffff&%23038;fg=000000&%23038;s=0"> and <img title="b" alt="b" class="latex" src="http://s0.wp.com/latex.php?latex=b&bg=ffffff&%23038;fg=000000&%23038;s=0"> parameters are mixing much more effectively than they were under the <code>RW_block</code> sampler.</p>
<div>
<div>
<div>
<pre><span>plot</span><span>(AFSS_mcmcOutput,</span> <span>density</span> <span>=</span> <span>FALSE</span><span>)</span>
</pre>
</div>
</div>
<div><img title="plot of chunk plot-mcmc-2" alt="plot of chunk plot-mcmc-2" src="https://i0.wp.com/r-nimble.org/wp-content/uploads/2017/05/plot-mcmc-2-1.png?w=456&ssl=1" data-recalc-dims="1"></div>
</div>
<p><cite><a></a><a target="_blank" rel="nofollow" href="https://r-nimble.org/better-block-sampling-in-mcmc-with-the-automated-factor-slice-sampler#cite-tibbits14">Tibbits, M. M, C. Groendyke, M. Haran, et al.</a><br>
(2014).<br>
“Automated factor slice sampling”.<br>
In: <em>Journal of Computational and Graphical Statistics</em> 23.2, pp. 543–563.</cite></p>
<p><cite><a></a><a target="_blank" rel="nofollow" href="https://r-nimble.org/better-block-sampling-in-mcmc-with-the-automated-factor-slice-sampler#cite-turek17">Turek, D, P. de<br>
Valpine, C. J. Paciorek, et al.</a><br>
(2017).<br>
“Automated parameter blocking for efficient Markov chain Monte Carlo sampling”.<br>
In: <em>Bayesian Analysis</em> 12.2, pp. 465–490.</cite></p>
</div>
<p> </p>
</div>
<p> <img height="1" width="1" src="https://r-nimble.org/?feed-stats-post-id=668"></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="https://r-nimble.org/better-block-sampling-in-mcmc-with-the-automated-factor-slice-sampler"> R – NIMBLE</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/search/tutorial">tutorials</a> on topics such as: <a title="Data science" rel="nofollow" href="https://www.r-bloggers.com/search/data%20science">Data science</a>, <a title="Big Data" rel="nofollow" href="https://www.r-bloggers.com/search/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/search/ggplot2">ggplot2</a>, <a title="Boxplots using lattice and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/boxplot">Boxplots</a>, <a title="Maps and gis" rel="nofollow" href="https://www.r-bloggers.com/search/map">maps</a>, <a title="Animation in R" rel="nofollow" href="https://www.r-bloggers.com/search/animation">animation</a>), programming (<a title="RStudio IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/RStudio">RStudio</a>, <a title="Sweave and literate programming" rel="nofollow" href="https://www.r-bloggers.com/search/sweave">Sweave</a>, <a title="LaTeX in R" rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX">LaTeX</a>, <a title="SQL and databases" rel="nofollow" href="https://www.r-bloggers.com/search/SQL">SQL</a>, <a title="Eclipse IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/eclipse">Eclipse</a>, <a title="git and github, Version Control System" rel="nofollow" href="https://www.r-bloggers.com/search/git">git</a>, <a title="Large data in R using Hadoop" rel="nofollow" href="https://www.r-bloggers.com/search/hadoop">hadoop</a>, <a title="Web Scraping of google, facebook, yahoo, twitter and more using R" rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping">Web Scraping</a>) statistics (<a title="Regressions and ANOVA analysis tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/regression">regression</a>, <a title="principal component analysis tutorial" rel="nofollow" href="https://www.r-bloggers.com/search/PCA">PCA</a>, <a title="Time series" rel="nofollow" href="https://www.r-bloggers.com/search/time+series">time series</a>, <a title="finance trading" rel="nofollow" href="https://www.r-bloggers.com/search/trading">trading</a>) and more...
<img height="1" alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/WSr15GZt-3Q"></a></div>