JAGS

Looks like JAGS 4.0 is out.

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. 

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.

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.

JAGS resources

Very good resource list of using JAGS for Bayesian computing

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.

JAGS tutorial

Here is a nice tutorial for JAGS (pdf format).

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.

JAGS tutorial

This is a nice tutorial for JAGS. Here is another one.

Runjags

To my knowledge, there are three interface packages between R and JAGS: rjags, R2jags, and runjags. I just realized that with runjags, you don't have to have a separate JAGS model file but to embed it in your R source file. This certainly has its appeal and might be considered as a unique advantage of runjags.

JAGS tutorial

Here is a JAGS tutorial, very helpful.

A simple JAGS demo

Here is a simple JAGS demo. The "progress.bar="gui"" option for the "coda.samples()" function is very convenient, but was not documented in the manual. 

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.

JAGS' glm module

I did not realize that new GLM module for JAGS 2.0 can be so effective in reducing the running time. Using the example I posed it yesterday (an example from the glmmBUGS package), without loading the glm module, it took 20.223 minutes to draw 10,000 samples from the posterior distribution; with the glm module loaded, it took 11.869 to get the same results:

Without glm module:
   user  system elapsed
 20.170   0.020  20.223 

With glm module:
   user  system elapsed
 11.040   0.000  11.869 

JAGS for Bayesian analysis

One of the many reasons I like Jackman's book "Bayesian Analysis for the Social Sciences" is that it relies on JAGS as the primary computational platform (as opposed to WinBUGS). From what I know, it is the first Bayesian book that does this.

The new JAGS is out

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

JAGS

Looks like a minor bug fix release of JAGS is out. I wonder whether the planned new major release will have some facilities to handle skewed distributions.

An easy way to get additional distribution in BUGS and JAGS

http://mathstat.helsinki.fi/openbugs/data/Docu/Tricks.html

Also, it seems that the JAGS development is moving right along: http://calvin.iarc.fr/wordpress/?p=15

Rjags on 64-bit Linux

Use this: install.packages("rjags",configure.args="--with-jags-lib=/usr/local/lib")