# Fitting a simple JAGS model with RSTAN [closed]

I'm trying to fit a simple exponential model for left censored data using RSTAN to replicate something I did in JAGS.

The JAGS model is:

model{
for(i in 1:N){
isAboveLOD[i] ~ dinterval(x[i], LOD[i])
x[i] ~ dexp(lambda)
}

lambda ~ dgamma(0.001, 0.001)
}


A full explanation of the simulated data can be found at http://jamescurran.co.nz/2014/06/bayesian-modelling-of-left-censored-data-using-jags/, but the essence of it is that the data had 9691 censored values (out of 10,000) and that the data was originally sampled from an Exp(1.05) distribution, with left censoring performed at log(29).

When I run this model with RSTAN, I essentially get no movement in the chain, and hence the estimate of lambda is way off (0.3 in my first run). I'm sure the error is elementary, but I'd like your advice. My RSTAN code (including data simulation) is given below.

expCode = '
data {
int<lower=0> nObs;
int<lower=0> nCens;
real<lower=0> yObs[nObs];
real<lower=0> U;
}
parameters {
real<lower=0, upper=U> yCens[nCens];
real<lower=0> lambda;
}
model {
for (n in 1:nObs)
yObs[n] ~ exponential(lambda);
for (n in 1:nCens)
yCens[n] ~ exponential(lambda);
lambda ~ gamma(0.001, 0.001);
}
'

set.seed(35202)
x = rexp(10000, rate = 1.05)

## set all the censored values to NA's
x[x < log(29)] = NA

yObs = x[!is.na(x)]
yCens = x[is.na(x)]

dataList = list(nObs = length(yObs), nCens = length(yCens), U = log(29), yObs = yObs)
library(rstan)
fit = stan(model_code = expCode, data = dataList, iter = 1000, chains = 1)

• In the future, there are lots and lots of questions like this already asked on the stan-users mailing list. Might be worth checking that one out. Jul 14, 2014 at 15:07

In BUGS/JAGS, the order in which statements are written does not matter. In Stan statements execute in the order in which they are written (see Stan 2.2.0 Reference Manual, pg. 405). Thus your last statement is in the wrong place: lambda is sampled from a gamma distribution, but that happens after the previous statements, so it's sampled 'in the air'.

Furthermore, in Stan

• if no prior is specified for a parameter, it is implicitly given a uniform prior on its support; thus you can define a vague gamma prior for $\lambda$, but that's not needed, real<lower=0> lambda is enough (pg. 8);
• the loop over sampling statements can be vectorized by replacing for (n in 1:N) y[n] ~ ... with the equivalent vectorized form y ~ ... (pg. 9).

If you replace your model block with

  model {
yObs ~ exponential(lambda);
yCens ~ exponential(lambda);
}


you get:

> print(fit, pars="lambda", digits=3)
...
mean se_mean    sd  2.5%   25%   50%   75% 97.5% n_eff  Rhat
lambda 1.033   0.001 0.015 1.005 1.023 1.034 1.044 1.064   382 0.999
...


EDIT: Please, stop upvoting my answer ;-) There is something wrong in what I've written, but I need a true computer to check what happens here (I'm using a small netbook now). I'll edit my answer as soon as possible.

• Thanks so much! I'm sorry I can't upvote this but my reputation isn't high enough! Jun 28, 2014 at 5:27
• I up voted it for you. I have a question though: in some of the Stan examples in the manual (for example, p. 34), a variable is first used before its prior appears (e.g. alpha ~ normal(0,sigma_alpha); appearing before sigma_alpha ~ cauchy(0,5);). Why is this sensible concerning the order of statements matters? Shouldn't the prior for a parameter always appear before the sampling statement for the parameter?
– jona
Jul 12, 2014 at 15:54
• @jona, you're right. I'm edit my answer as soon as possible. Jul 14, 2014 at 14:44