[Rcpp-devel] RcppArmadillo -- crude benchmark of various mcmcvs rcppbugs
Smith, Dale
Dale.Smith at fiserv.com
Fri Apr 27 17:26:30 CEST 2012
As the MCMCpack and rcppbugs packages were built on Linux, I can't install them on my Windows box. I'll have to get the sources and try to make them when I have time.
Thanks,
Dale Smith, Ph.D.
Senior Financial Quantitative Analyst
Risk & Compliance
Fiserv.
107 Technology Park
Norcross, GA 30092
Office: 678-375-5315
Mobile: 678-982-6599
Mail: dale.smith at fiserv.com
www.fiserv.com
-----Original Message-----
From: rcpp-devel-bounces at r-forge.wu-wien.ac.at [mailto:rcpp-devel-bounces at r-forge.wu-wien.ac.at] On Behalf Of Dirk Eddelbuettel
Sent: Thursday, April 26, 2012 11:28 PM
To: Whit Armstrong
Cc: rcpp-devel
Subject: Re: [Rcpp-devel] RcppArmadillo -- crude benchmark of various mcmcvs rcppbugs
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 (two-year 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["elaps
| ed"],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[,"M
| ean"][1:2],jags.coefs)
| rownames(coef.compare) <- c("lm","rcppbugs","mcmcpack","jags")
| print(coef.compare)
| _______________________________________________
| Rcpp-devel mailing list
| Rcpp-devel at lists.r-forge.r-project.org
| https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-deve
| l
--
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 _______________________________________________
Rcpp-devel mailing list
Rcpp-devel at lists.r-forge.r-project.org
https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
More information about the Rcpp-devel
mailing list