Nimble 0.4 is out

Quite a number of new features. More details here.

Why WinBUGS and OpenBUGS do not have a future

As this post explains.

The BUGS Book

THE BUGS book. No need to say more.

Comparison between OpenBUGS, JAGS, and Stan

I have been playing with the examples of Introduction to WinBUGS for Ecologists in the past two days.  One thing that concerns me particularly is how well these Bayesian packages handle large data set (i.e., N>10,000), which is the size of the data sets that I work with. So I modified the data generation program provided by the authors and increased the sample sizes:


n1 <- 6000 # Number of females
n2 <- 4000 # Number of males
mu1 <- 105 # Population mean of females
mu2 <- 77.5 # Population mean of males
sigma <- 2.75 # Average population SD of both

n <- n1+n2 # Total sample size
y1 <- rnorm(n1, mu1, sigma) # Data for females
y2 <- rnorm(n2, mu2, sigma) # Date for males
y <- c(y1, y2) # Aggregate both data sets
x <- rep(c(0,1), c(n1, n2)) # Indicator for male

The BUGS program looks like this:


model {
# Priors
 mu1 ~ dnorm(0,0.001)
 mu2 ~ dnorm(0,0.001)
 tau1 <- 1 / ( sigma1 * sigma1)
 sigma1 ~ dunif(0, 1000) # Note: Large var. = Small precision
 tau2 <- 1 / ( sigma2 * sigma2)
 sigma2 ~ dunif(0, 1000)

# Likelihood
 for (i in 1:n1) {
    y1[i] ~ dnorm(mu1, tau1) 
 }

 for (i in 1:n2) {
    y2[i] ~ dnorm(mu2, tau2) 
 }

# Derived quantities
 delta <- mu2 - mu1
}

While the Stan looks like this:

data {

     int n1;

     int n2;

     real y1[n1];

     real y2[n2];

}

parameters {

           real alpha1;

           real alpha2;

           real sigma1;

           real sigma2;

}

model {

    y1 ~ normal(alpha1, sigma1);          

    y2 ~ normal(alpha2, sigma2);

}

generated quantities {

          real dif;

          dif <- sigma1 - sigma2;

}

For 20,000 iterations, OpenBUGS took 7852.813 seconds, JAGS took 2252.273 seconds, Stan took 60-150 seconds (the Stan team is still trying to iron out some bugs, which may explain the wide range of run time). While the run time of both OpenBUGS and JAGS are sensitive to the sample size, the run time of Stan seems much less sensitive to sample size.

However, if the Stan code is written as this:


model {
      for (n in 1:n1)
          y1[n] ~ normal(alpha1, sigma1);
      for (n in 1:n2)
          y2[n] ~ normal(alpha2, sigma2);
}


Then Stan has no apparent performance advantage over JAGS. 

Applied Bayesian Statistics

This book looks interesting.

Stan code for my favorite book

I am learning Stan. So I decided to translate the BUGS code for my favorite applied Bayesian book into  Stan code. I just did the translation for a couple of chapters and I feel I know Stan much better already.

From BUGS with BRugs to JAGS with rjags

John Kruschke provided helpful suggestions about how to convert BUGS programs into JAGS. It would be interesting to see some suggestions about how to convert BUGS or JAGS programs into Stan.

Some thoughts on Stan

Stan is very promising. The glmmBUGS package should be easily extended to produce Stan code in place of or in addition to BUGS/JAGS code, which will makes it even easier for novice uses to get started.

Stan

Stan is a new software package for Bayesian analysis. It also comes with an R interface. According to this benchmark, it may be a viable choice for real world data analysis. I took a similar approach as ADMB by first converting a syntax file into C++ source files and then generating native machine code.

New version of g++, a good reason to upgrade

My Ubuntu 10.04 runs fine on my workstation except for one thing ... its g++4.4.3 does not have the new c++ language feature which is used in software like CppBugs and its R binding called rcppbugs. The idea behind these two packages are very interesting, which gives me a good reason to upgrade to a new version of the OS (Debian, Mint, or Ubuntu).

A practical Bayesian book

Introduction to WinBUGS for Ecologists is a good practical Bayesian book. I find the sample code very useful. The simple example of forming predictions (p. 116) is straightforward and can be easily tailored to be used in different research situations.

Doing Bayesian Data Analysis: A Tutorial with R and BUGS is another useful practical Bayesian book. My copy just arrived today. Can't wait to play with some of the examples there.

CppBugs

A simple way to get CppBugs to work is to copy the "cppbugs" folder to "/usr/local/include". On my Ubuntu 10.04 system, I also need to modify the make file located in the test file directory and replace the "ARMADILLO_LIBS = llapack" with "ARMADILLO_LIBS = -larmadillo".

As a BUGS/JAGS user, I think this project is very promising and extremely interesting. It will be even more attractive if CppBugs can be integrated with R (maybe through Rcpp) in the similar way that R2WinBUGS integrates WinBUGS with R.

The "rube" package

The "rube" package is very promising. It builds upon the "R2WinBUGS" package and adds some very useful features. I hope it further extends to the JAGS engine because JAGS is apparently the future of the BUGS famine.

Use JAGS with ggplot2

Here is a short demo of how to use JAGS with ggplot2.

WinBUGS to ADMB

I found the short article published on the ADMB newsletter about how to translate WinBUGS code into ADMB very educational.

OpenBUGS and rbugs package

I just realized that OpenBUGS 3.10 has a native Linux version and it can be called directly from within R via the package "rbugs". OpenBUGS is probably closer to WinBUGS than JAGS, and it is always nice to alternatives!

It will be fun to do some benchmark test between OpenBUGS and JAGS. 

Making glmmBUGS work with JAGS

To streamline the process of making glmmBUGS to work with JAGS, I made changes to the "writeBugsModel.R" and replaced all "inprod2" with "inprod" and "dflat()" with "dnorm(0.0,1.0E-6)". I will spend some time looking at other sources files and see if anything else should be changed.

A question

I posted a question to the R-Mixed-Model mailing list but have not had any responses yet:
---------------------------------------------------------------------------------------

Dear All,

I wonder if anybody has tried to make glmmBUGS work with JAGS. My
attempt was not successful.

Here is the simple example I copied from the glmmBUGS Vignettes:

------------------------
library(MASS)
data(bacteria)

bacterianew <- bacteria
bacterianew$yInt = as.integer(bacterianew$y == "y")

levels(bacterianew$trt) <- c("placebo", "drug", "drugplus")

library(glmmBUGS)

bacrag <- glmmBUGS(formula = yInt ~ trt + week, data=bacterianew,
effects = "ID",
                   modelFile="model.bug", family="bernoulli")

names(bacrag$ragged)

source("getInits.R")

startingValues = bacrag$startingValues
------------------------

With WinBUGS, it runs well:

------------------------
library(R2WinBUGS)

bacResult1 = bugs(bacrag$ragged, getInits, model.file="model.bug", n.chain=3,
  n.iter=2000, n.burnin=100, parameters=names(getInits()), n.thin=10)
------------------------

But with JAGS, I got error message "Error in FUN(50L[[1L]], ...) :
invalid first argument"

------------------------
library(R2jags)

jags.parms=names(getInits())

bacResult = jags(data=bacrag$ragged, n.chain=3, n.iter=2000,
model.file=model.bug)
------------------------

I hope somebody can help me figuring out how to make this work. Many thanks.

Best,
Shige

---------------------------------------------------------------------------------------

Here is the answer from Jens Åströ
---------------------------------------------------------------------------------------
Hi!

Not an expert by any means, but I got it to run by doing this:

1) change inprod2 to inprod in the model file (JAGS does not have
inprod2 function)

2) change ~dflat() to other uninformative prior, e.g.
~dnorm(0.0,1.0E-6) (Jags does not have dflat distribution)

3) specify/compile the model with e.g.
bac.jags<-jags.model("model.bug",data=bacrag$ragged,n.chains=4)

4) update it, update(bac.jags,1000)

5) Collect coda samples,
bacResult<-coda.samples(temp,names(getInits()),n.iter=10000,thin=10)

This seemed to converge well, (gelman.diag(bacResult))

Good luck
----------------------------------------------------------------------------------------

It works.

The new JAGS is out

Version 2.0 of the cross-platform BUGS language implementation is out.

BUGS language: strength and weakness

Here is an informative introduction to the strength and weakness of the BUGS language in conducting Bayesian data analysis. Since my analysis often utilize large data (tens of thousands, even millions), BUGS may not be a viable option in many cases. I hope the soon-to-come new version of JAGS can overcome some of the weaknesses mentioned in Gelman's post.