Skip to content

Instantly share code, notes, and snippets.

@yuu-ito
Last active August 29, 2015 14:20
Show Gist options
  • Save yuu-ito/24b75e6f900e589b9d06 to your computer and use it in GitHub Desktop.
Save yuu-ito/24b75e6f900e589b9d06 to your computer and use it in GitHub Desktop.
# 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)
                 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倍の速度(という意味?
@yuu-ito
Copy link
Author

yuu-ito commented May 8, 2015

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment