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 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.

Tweet- Serialize and Deserialize a C++ Object in Rcpp — Wush Wu
- Stochastic SIR Epidemiological Compartment Model — Christian Gunning
- Call matplotlib from R — Matti Pastell