When Gibbs sampling in R is slow, we can use Rcpp to speed it up by integrating C++ code into R. However, Rcpp code is harder to debug because we cannot run Rcpp code line-by-line like R. It’s paramount that we ensure the Rcpp code is doing what we wants it to.
The solution is to start with small Rcpp programs and make sure that they produce the same results as their pure R counterpart. Following this advice, I’m writing a series of increasingly complex Markov chain Monte Carlo (MCMC) code in Rcpp, eventually building towards my final goal of a custom Metropolis-Hastings algorithm for my work on modeling FDI flow.
The first installment will be a Gibbs sampler of a semi-conjugate normal model, arguably the most common introductory model to Bayesian statistics and MCMC. (Both Hoff’s First Course in Bayesian and Gelman’s Bayesian Data Analysis cover it.)
I’ll follow Hoff’s example in Chapter 6 as a ground truth to check whether my code is working correctly. For Rcpp basics, I read Hadley Wickham’s introduction to Rcpp. I’ll note additional pitfalls along the way of converting R to Rcpp.
Normal model with semi-conjugate priors
A normal model with semi-conjugate priors has the following specifications. For more details and intuition behind the model, see Chapter 5 & 6 in Hoff’s book.
There is a quick intuitive interpretation of the full conditional of ,
. Its posterior mean is the weighted sum of the prior mean and the sample mean . This makes sense if we think of the posterior as combining the information from the prior and the sample.
But how should we weigh the information from the prior and the sample? Intuitively, if there’s more information in the prior, we should weigh it more. The amount of information in the prior is represented by the prior precision, . As the formula shows, the precision is the inverse of variance, because the smaller the variance there is in our prior belief about , the more information we have in the prior. In sum, because precision represents the amount of information available, we should weigh the prior mean by the prior precision . That is indeed the weight of the prior mean in the weighted sum.
Similarly, we weigh the sample mean by the sample precision . Note that we multiply the precision by to capture the fact that the bigger the sample size, the more information there is in the sample.
Pure R Gibbs sampler
Below is a direct translation of the full conditionals below into a Gibbs sampler. Inside the big loop, for the iteration we sample:
We confirm that our algorithm successfully reproduces the result in Hoff (p. 95) with exactly the same posterior distribution of theta.
Rcpp Gibbs sampler
We now rewrite the Gibbs sampler in Rcpp. Putting the R and Rcpp codes side by side, we see only a few differences thanks to Rcpp’s syntactic sugar that replicates R-like syntax in C++. Most of these differences are covered in Hadley’s Intro to Rcpp. Beyond that, there are 2 additional pitfalls when transitioning from R to Rcpp.
Remember to use
1while doing division so that we get double division like in R, not integer division. For example, in Rcpp,
1 / 2 = 0while
1.0 / 2 = 0.5.
rgammain R uses two parameters, location and rate. In contrast,
rgammain Rcpp uses location and scale, which is the inverse of rate.
We now check and confirm that our Gibbs sampler in Rcpp produces exactly the same result as its pure R counterpart.
Comparing the density plots on the left with those on the right, we see that the posterior distributions of and are exactly the same for the R and Rcpp implementations.
The Rcpp code is also a lot faster than R as expected. The median running time of
f_gibbsC is about 38 times faster than
And that’s the one of the simplest Gibbs sampler in Rcpp. In the next installment we will implement the Metropolis-Hastings sampler in Rcpp. You can click on the categories Bayesian-sampler-in-Rcpp to find more entries in this series.