[Rcppdevel] RcppArmadillo  crude benchmark of various mcmc vs rcppbugs
Dirk Eddelbuettel
edd at debian.org
Fri Apr 27 05:27:55 CEST 2012
On 26 April 2012 at 15:58, Whit Armstrong wrote:
 Slightly off topic, but Dirk kindly allowed me to make a quick
 benchmark post on rcppbugs.
"Allowed" ? Come on, benchmarks are always ok, as are posts on Rcpp. This is
after all the Rcpp list...

 Just a short post to the list showing some performance numbers (for
 running a simple linear regression) of rcppbugs (with an RcppArmadillo
 core) vs MCMCpack (Scythe core) vs JAGS (?? core).

 I was fairly optimistic about the speed of rcppbugs, yet this exceeds
 my expectations by so much that I fear there may be a bug in my
 benchmark code.

 I'm gearing up for my presentation at R/Finance, so I would be
 grateful to anyone who wants to kick the tires a bit. The core has
 been tested (we've been using the pure c++ implementation for over a
 year), but the R binding (via rcppbugs) is completely new.

 Thanks in advance. Code below (you'll need to install MCMCpack,
 rjags, and rcppbugs (now on CRAN) to run it).

 The quick summary is a 5x speedup over MCMCpack and a 77x speedup (I'm
 not sure I believe this yet) over JAGS.
A bit more modest here (twoyear old i7 running Ubuntu amd64)
R> print(round(all.times.ratio,2))
time ratio
lm 0.09 0.08
rcppbugs 1.18 1.00
mcmcpack 4.05 3.45
jags 57.06 48.56
But still pretty good compared to MCMCpack. That is possibly legit  modern
C++ versus a C++ design from a decade ago. Jags is a little harder to
believe. Maybe time to write an email to Martyn?
Thanks for sharing this. Looking forward to your talk next month.
Dirk
 Comments welcome.

 Whit


 library(rcppbugs,quietly=TRUE,verbose=FALSE)
 library(rjags,quietly=TRUE,verbose=FALSE)
 library(MCMCpack,quietly=TRUE,verbose=FALSE)

 NR < 1e2L
 NC < 2L

 y < rnorm(NR,1) + 10
 X < matrix(nr=NR,nc=NC)
 X[,1] < 1
 X[,2] < y + rnorm(NR)/2  10


 lm.time < system.time(lm.res < lm.fit(X,y))
 ##print(coef(lm.res))

 b < mcmc.normal(rnorm(NC),mu=0,tau=0.0001)
 tau.y < mcmc.uniform(sd(as.vector(y)),lower=0,upper=100)
 ##y.hat < deterministic(function(X,b) { X %*% b }, X, b)
 y.hat < linear(X,b)
 y.lik < mcmc.normal(y,mu=y.hat,tau=tau.y,observed=TRUE)
 m < create.model(b, tau.y, y.hat, y.lik)


 iterations < 1e5L
 burn < iterations
 adapt < 1e3L
 thin < 10L

 cat("running rcppbugs model...\n")
 rcppbugs.time < system.time(ans < run.model(m,
 iterations=iterations, burn=burn, adapt=adapt, thin=thin))
 rcppbugs.coefs < apply(ans[["b"]],2,mean)
 ##print(rcppbugs.coefs)
 print(rcppbugs.time)

 cat("running MCMCpack model...\n")
 mcmcpack.time < system.time(mcmcpack.out < MCMCregress(y ~
 X.2,data=data.frame(y=y,X=X),burnin = burn, mcmc = iterations,
 thin=thin))
 print(mcmcpack.time)

 cat("running jags model...\n")

 bug.model < '
 model {
 for (i in 1:NR){
 y[i] ~ dnorm(y.hat[i], tau.y)
 y.hat[i] < inprod(b, X[i,])
 }
 for (j in 1:NC){
 b[j] ~ dnorm(0, .0001)
 }
 tau.y ~ dunif(0, 100)
 }
 '
 bug.file < tempfile(fileext=".bug")
 cat(bug.model,file=bug.file)

 jags.time1 < system.time(
 jags < jags.model(bug.file,
 data = list('X' = X,
 'y' = y,
 'NR' = NR,
 'NC' = NC),
 n.chains = 1,
 n.adapt = adapt,
 quiet=TRUE
 ))
 jags.time2 < system.time(update(jags, n.iter=burn, progress.bar="none"))
 jags.time3 < system.time(jags.trace <
 jags.samples(jags,c('b','tau.y'),n.iter=iterations,thin=thin,progress.bar="none"))

 jags.time < jags.time1 + jags.time2 + jags.time3
 print(jags.time)
 jags.coefs < apply(jags.trace$b,1,mean)

 ##print(jags.trace)

 cat("rcppbugs speedup:\n")
 all.times < rbind(lm.time["elapsed"],rcppbugs.time["elapsed"],mcmcpack.time["elapsed"],jags.time["elapsed"])
 all.times.ratio < cbind(all.times,all.times/rcppbugs.time["elapsed"])
 colnames(all.times.ratio) < c("time","ratio")
 rownames(all.times.ratio) < c("lm","rcppbugs","mcmcpack","jags")
 print(round(all.times.ratio,2))

 cat("coef comparison:\n")
 coef.compare <
 rbind(coef(lm.res),rcppbugs.coefs,summary(mcmcpack.out)$statistics[,"Mean"][1:2],jags.coefs)
 rownames(coef.compare) < c("lm","rcppbugs","mcmcpack","jags")
 print(coef.compare)
 _______________________________________________
 Rcppdevel mailing list
 Rcppdevel at lists.rforge.rproject.org
 https://lists.rforge.rproject.org/cgibin/mailman/listinfo/rcppdevel

R/Finance 2012 Conference on May 11 and 12, 2012 at UIC in Chicago, IL
See agenda, registration details and more at http://www.RinFinance.com
More information about the Rcppdevel
mailing list