[Rcpp-devel] RcppArmadillo -- crude benchmark of various mcmcvs rcppbugs

Smith, Dale Dale.Smith at Fiserv.com
Mon Apr 30 14:46:13 CEST 2012


Hi Whit,

I installed rcppbugs, MCMCpack, and rjags from CRAN and ran your R script. I was puzzled by a few errors, but sorted that out, I think, as the script was not formatted properly in my email. Here's what I'm using, just so I have the right code

***
...

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)

...

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)

...
***

I looked at the output and found no errors. Here are my results

> print(round(all.times.ratio,2))
          time ratio
lm        0.00  0.00
rcppbugs  0.27  1.00
mcmcpack  1.48  5.48
jags     15.41 57.07

I am running

Windows 7 Professional, SP1, 64 bit
Intel Core i5 M520 @ 2.40GHz

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: Whit Armstrong [mailto:armstrong.whit at gmail.com] 
Sent: Friday, April 27, 2012 11:34 AM
To: Smith, Dale
Cc: Dirk Eddelbuettel; rcpp-devel
Subject: Re: [Rcpp-devel] RcppArmadillo -- crude benchmark of various mcmcvs rcppbugs

The windows build was just accepted this morning by CRAN.

It should be appearing on the mirrors shortly.

It appears the mac builds are being held back because gcc 4.2 has some template issues for modern c++ code.  I'll try to find some #ifdef's to work around this issue.

-Whit


On Fri, Apr 27, 2012 at 11:26 AM, Smith, Dale <Dale.Smith at fiserv.com> wrote:
> 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["ela
> | ps
> | 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-de
> | ve
> | 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-deve
> l


More information about the Rcpp-devel mailing list