## Tuesday, October 23, 2012

## Thursday, October 18, 2012

### Comparison between OpenBUGS, JAGS, and Stan

I have been playing with the examples of

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

}

*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.

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.

## Tuesday, October 16, 2012

## Monday, October 15, 2012

### 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.

## Sunday, October 14, 2012

### 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.

## Saturday, October 13, 2012

### 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.

## Friday, October 12, 2012

### 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.

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

int

int

real x1[I];

real x2[I];

}

parameters {

real alpha0;

real alpha1;

real alpha12;

real alpha2;

real

real b[I];

}

transformed parameters {

real

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.

## Thursday, October 11, 2012

### 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.

## Wednesday, October 10, 2012

## Tuesday, October 09, 2012

### 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)

Labels:
Bayesian,
missing values,
prediction,
stan

Subscribe to:
Posts (Atom)