I want to speed up my R implementation of a Metropolis Hasting procedure by replacing the slow parts with functions written in Rcpp.

There are already some examples online using Rcpp to speed up Gibbs sampling and Metropolis Hastings.

However, these examples implement the entire sampler in Rcpp. In contrast, I hope to rewrite only the slower part of my R implementation so that I can keep the niceties of R as much as possible. (I'm trying to avoid premature optimization here.)

However, profiling my code does not reveal any particularly candidate for optimization. Most of the code already uses R's vectorized operation & built-in random number generator.

What parts of the Metropolis Hastings algorithm can be rewritten in Rcpp to achieve a speed up? Or do I indeed have to rewrite the entire loop in Rcpp itself?

On why coding the Metropolis Hastings by hand: One of my parameter is a binary matrix of 0s and 1s, and my understanding is that Stan / JAGS don't accommodate such kind of variable.

On why this is a statistics question: I already know how to use Rcpp or to do code profiling, so I'm not limited by programming ability. The question requires understanding MCMC procedures, especially which part can be ripped out and put in Rcpp, and which part cannot. That understanding is a lot more about statistics than programming.

P/S: I can't answer this question myself partially because I don't understand why R is slower than C++ / Rcpp, and thus can't reason through which part of Metropolis-Hastings would be faster in Rcpp. If the answer touches on the reason why a certain part is better written in Rcpp, it'd be very instructive.

closed as off-topic by Dougal, kjetil b halvorsen, mdewey, Michael Chernick, Stephan Kolassa Feb 23 at 7:03

This question appears to be off-topic. The users who voted to close gave this specific reason:

  • "This question appears to be off-topic because EITHER it is not about statistics, machine learning, data analysis, data mining, or data visualization, OR it focuses on programming, debugging, or performing routine operations within a statistical computing platform. If the latter, you could try the support links we maintain." – mdewey, Michael Chernick, Stephan Kolassa
If this question can be reworded to fit the rules in the help center, please edit the question.

  • 1
    I think you might get better answers on a programming site. Note that what you want to do is to profile your code to see where it is spending all its time as there is no point optimising the other bits (yet). I believe there are many tools in R for doing profiling but I have never used them myself so cannot offer any more. – mdewey Feb 22 at 14:41
  • I debated whether to post on SO, and thought that here is better. I added my reasoning for asking this on Cross Validated. – Heisenberg Feb 22 at 16:08
  • 1
    As it stands, this question is quite general right now, because it really depends on the kind of target distribution you have, and how complicated the M-H ratio is. One thing that I understand of Rcpp is that it loops much faster than R. So since M-H cannot be vectorized and has to be looped, I would imagine coding the whole thing in Rcpp will probably be the best way to speed it up. – Greenparker Feb 22 at 16:18
  • Why do you want to code it by hand rather then using software like JAGS or Stan? MH is a pretty simple algorithm so since lops in C++ are faster then in R then re-wiring it in C++ will speed it up. If you have some special needs, then maybe there are problem specific optimizations, but you didn't give us much details... – Tim Feb 22 at 16:51
  • @Tim I have to code it by hand because one of my parameter is a binary matrix of 0s and 1s, and my understanding is that Stan can't handle that kind of discrete variable. I've added this to the question. Thank you for the suggestion. – Heisenberg Feb 22 at 17:07

To elaborate on my comment about loops, here is a benchmark of looping in R and looping in Cpp. I construct a function in both that sums up numbers upto p. The Rcpp file:

// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]
double summingC(double p)
{
  double add_sum = 0;

  for(double i = 1; i <= p; ++i)
  {
    add_sum = add_sum + i;
  }
 return(add_sum);
}

The R code and benchmarking:

set.seed(100)

library(Rcpp)
library(rbenchmark)

# Loading Rcpp file
sourceCpp("summing.cpp")

add_sum <- 0
summingR <- function(p)
{
for(i in 1:p)
{
    add_sum <- add_sum + i
}
return(add_sum)
}

## benchmarking between Cpp and R
benchmark(summingC(1e7), summingR(1e7))

And the output below indicates that even for doing a simple addition, Cpp is 15 times faster.

> benchmark(summingC(1e7), summingR(1e7))
             test replications elapsed relative user.self sys.self user.child
1 summingC(1e+07)          100   1.505     1.00     1.471    0.004          0
2 summingR(1e+07)          100  22.740    15.11    22.654    0.022          0

So it seems best to put all loops in Rcpp.

Not the answer you're looking for? Browse other questions tagged or ask your own question.