# Gibbs Sampler in C++

**Rcpp Gallery**, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

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 1.161 1.00 2 RCgibbs(n, thn) 10 45.763 39.42 1 Rgibbs(n, thn) 10 57.995 49.95

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.

**leave a comment**for the author, please follow the link and comment on their blog:

**Rcpp Gallery**.

R-bloggers.com offers

**daily e-mail updates**about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.