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:


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. 

No comments: