Gibbs Sampler in C++

Dirk Eddelbuettel — written Aug 3, 2013 — source

Markov Chain Monte Carlo (MCMC) is a popular simulation method. As it is somewhat demanding, it is also frequently used to benchmark different implementations or algorithms.

One particular algorithm has been compared a number of times, starting with an article by Darren Wilkinson, and Darren’s follow–up article which in turns responded in part to our article.

This post simply refreshes the implementation using Rcpp attributes.

First we look at the R version, and its byte-compiled variant.

## Here is the actual Gibbs Sampler
## This is Darren Wilkinsons R code (with the corrected variance)
Rgibbs <- function(N,thin) {
    mat <- matrix(0,ncol=2,nrow=N)
    x <- 0
    y <- 0
    for (i in 1:N) {
        for (j in 1:thin) {
            x <- rgamma(1,3,y*y+4)
            y <- rnorm(1,1/(x+1),1/sqrt(2*(x+1)))
        }
        mat[i,] <- c(x,y)
    }
    mat
}

## We can also let the R compiler on this R function
library(compiler)
RCgibbs <- cmpfun(Rgibbs)

Creating a version in C++ is very straightforward thanks to Rcpp and Rcpp Attributes. It transfers the integer arguments n (number of simulations) and thn (number of extra thinning simulations), initializes the R random number generator for us (eg no need to manual nstantiate the RNGScope object), and returns the result matrix

Also, since the initial posts were written, we made the (scalar) RNGs of the R API available directly via the R namespace. This is both little nice to read than the poor-man’s pseudo-namespace in C via the Rf_ prefix, and actually provides proper C++ namespace.

// load Rcpp
#include <Rcpp.h>

using namespace Rcpp;		// shorthand

// [[Rcpp::export]]
NumericMatrix RcppGibbs(int n, int thn) {

    int i,j;
    NumericMatrix mat(n, 2);

    // The rest of the code follows the R version
    double x=0, y=0;

    for (i=0; i<n; i++) {
        for (j=0; j<thn; j++) {
            x = R::rgamma(3.0,1.0/(y*y+4));
            y = R::rnorm(1.0/(x+1),1.0/sqrt(2*x+2));
        }
        mat(i,0) = x;
        mat(i,1) = y;
    }

    return mat;             // Return to R
}

With the functions in place, we can re-run the benchmark.

library(rbenchmark)
n <- 2000
thn <- 200
benchmark(Rgibbs(n, thn),
          RCgibbs(n, thn),
          RcppGibbs(n, thn),
          columns=c("test", "replications", "elapsed", "relative"),
          order="relative",
          replications=10)

               test replications elapsed relative
3 RcppGibbs(n, thn)           10   0.371    1.000
2   RCgibbs(n, thn)           10  13.677   36.865
1    Rgibbs(n, thn)           10  13.843   37.313

As we have seen before, the C++ version is about 50 times faster, and around 40 times faster than the byte-compiled version.

A related article here on the Rcpp Gallery looks into timing different RNG implementation as this study revealed that the generator for Gamma-distributed random number in R is not particularly fast.

tags: rng sugar

Extending R with C++ and Fortran — Dirk Eddelbuettel and JBrandon Duck-Mayr
Benchmarking Rcpp code with RcppClock — Zach DeBruine
Simulation Smoother using RcppArmadillo — Tomasz Woźniak

Gibbs Sampler in C++

Related Articles