The BUGS Book

THE BUGS book. No need to say more.

Lubuntu

I played with Lubuntu for a while and really liked its minimalist and yet elegant design. I think it should be my main desktop OS.

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. 

Processing Stan samples with ggmcmc

Very helpful post.

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.

Insightful analysis of the recent quarrel between China and Japan

I've seen lots of economic, political, and social analysis about it, but this cultural analysis is very interesting and insightful.

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 and INLA

Since both Stan and INLA provide their own implementation of the classical BUGS examples, it is very illustrative and educational to compare the two sets of codes.

For example, the INLA code for the "seeds" example (my of my favorites) is like this:


library(INLA)

data(Seeds)
formula = r ~ x1*x2+f(plate,model="iid")
mod.seeds = inla(formula,data=Seeds,family="binomial",Ntrials=n)

## improved estimation of the hyperparameters
h.seeds = inla.hyperpar(mod.seeds)


while the Stan code is like this:


data {
    int I;
    int n[I];
    int N[I];
    real x1[I];
    real x2[I];
}
parameters {
    real alpha0;
    real alpha1;
    real alpha12;
    real alpha2;
    real tau;
    real b[I];
}
transformed parameters {
    real sigma;
    sigma  <- 1.0 / sqrt(tau);
}
model {
   alpha0 ~ normal(0.0,1.0E3);
   alpha1 ~ normal(0.0,1.0E3);
   alpha2 ~ normal(0.0,1.0E3);
   alpha12 ~ normal(0.0,1.0E3);
   tau ~ gamma(1.0E-3,1.0E-3);

   b ~ normal(0.0, sigma);

   for (i in 1:I) {
      n[i] ~ binomial(N[i], inv_logit(alpha0 
                                      + alpha1 * x1[i] 
                                      + alpha2 * x2[i]
                                      + alpha12 * x1[i] * x2[i] 
                                      + b[i]) );
   }
}

INLA has a more R-like syntax and is much more terse, but Stan is much more flexible and can handle very complicated models that INLA may have troubles with. 

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.

A few Stan tutorials

Here is a list of Stan tutorials:

1. Comparing Stan to JAGS for Bayesian Inference 
2. Review of STAN: off-the-shelf Hamiltonian MCMC
3. Stan versus JAGS
4. R, Stan and NUTS: A Hierarchical Linear Model Example Linear Model Example

Prediction, missing data, etc. in Stan

library(rstan)

N <- 1001

N_miss <- ceiling(N / 10)

N_obs <- N - N_miss

mu <- 3

sigma <- 2

y_obs <- rnorm(N_obs, mu, sigma)

missing_data_code <-

'

data {

  int N_obs;

  int N_miss;

  real y_obs[N_obs];

}

parameters {

  real mu;

  real sigma;

  real y_miss[N_miss];

}

model {

  // add prior on mu and sigma here if you want

  y_obs ~ normal(mu,sigma);

  y_miss ~ normal(mu,sigma);

}

generated quantities {

  real y_diff;

  y_diff <- y_miss[101] - y_miss[1];

}

'

results <- stan(model_code = missing_data_code,

                data = list(N_obs = N_obs, N_miss = N_miss, y_obs = y_obs))

y_diff <- apply(extract(results, c("y_miss[1]", "y_miss[101]")), 1:2, diff)