_rcpp_sandbox.R

# rcpp sandbox

gibbs_r <- function(n,thin){
  # R実装のギブスサンプラー
  mat <- matrix(0,nrow=n,ncol=2)
  x <- 0
  y <- 0
  for(i in 1:n){
    for(j in 1:thin){
      x <- rgamma(1, 3, 1/(y*y+4) )
      y <- rgamma(1, 1/(x+1), 1/(y*y+4) )
    }
    mat[i, ] <- c(x,y)
  }
  return(mat)
}

library("Rcpp")
library("inline")

# CPP実装のギブスサンプラー
src <-"
Rcpp::NumericMatrix mat(as<int>(n), 2);
double x=0, y=0;
for(int i=0; i < as<int>(n) ; ++i){
  for(int j=0; j < as<int>(thin) ; ++j){
    x = R::rgamma(3.0, 1.0/(y*y+4));
    y = R::rgamma(1.0/(x+1), 1.0/sqrt(2*x+2));
  }
  mat(i, 0) = x;
  mat(i, 1) = y;
}
return(mat);
"

# src <- 'return(Rcpp::wrap(n));'
gibbs_rcpp <- cxxfunction(signature(n = "integer", thin = "integer" ),src,plugin="Rcpp")


library("rbenchmark")
n <- 1000
thin <- 10

benchmark(gibbs_r(n,thin),
          gibbs_rcpp(n,thin),
          columns = c("test","replications","elapsed","relative"),
          order="relative",
          replications=100)

output.md

                 test replications elapsed relative
2 gibbs_rcpp(n, thin)          100    0.47    1.000
1    gibbs_r(n, thin)          100    8.12   17.277

• 17倍の速度（という意味？

参考
http://www.slideshare.net/masakitsuda940/rcpp-36654397