Ajuste de un modelo JAGS sencillo con RSTAN

Question

Estoy tratando de ajustar un modelo exponencial simple para datos censurados a la izquierda usando RSTAN para replicar algo que hice en JAGS.

El modelo JAGS es:

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

  lambda ~ dgamma(0.001, 0.001)
}

Encontrará una explicación completa de los datos simulados en http://jamescurran.co.nz/2014/06/bayesian-modelling-of-left-censored-data-using-jags/ pero la esencia es que los datos tenían 9691 valores censurados (de 10.000) y que los datos se muestrearon originalmente a partir de una distribución Exp(1,05), con censura a la izquierda realizada en log(29).

Cuando ejecuto este modelo con RSTAN, esencialmente no obtengo ningún movimiento en la cadena y, por tanto, la estimación de lambda es muy errónea (0,3 en mi primera ejecución). Estoy seguro de que el error es elemental, pero me gustaría que me aconsejarais. Mi código RSTAN (incluyendo la simulación de datos) se da a continuación.

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)

Answer 1

En BUGS/JAGS, el orden en que se escriben las declaraciones no importa. En Stan las sentencias se ejecutan en el orden en que están escritas (véase el Manual de referencia de Stan 2.2.0, pág. 405). Por lo tanto, su última declaración está en el lugar equivocado: lambda se muestrea a partir de una distribución gamma, pero eso sucede después de las declaraciones anteriores, por lo que se muestrea "en el aire".

Además, en Stan

• si no se especifica una prioridad para un parámetro, se le asigna implícitamente una prioridad uniforme sobre su soporte. puede definir una prioridad gamma imprecisa para $\lambda$ pero no es necesario, real<lower=0> lambda es suficiente (pág. 8);
• el bucle sobre las sentencias de muestreo puede vectorizarse sustituyendo for (n in 1:N) y[n] ~ ... con la forma vectorizada equivalente y ~ ... (pág. 9).

Si sustituye el bloque modelo por

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

lo consigues:

> 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: Por favor, deja de votar mi respuesta ;-) Hay algo mal en lo que he escrito, pero necesito un ordenador de verdad para comprobar lo que pasa aquí (ahora estoy usando un pequeño netbook). Editaré mi respuesta lo antes posible.

Answer 2

No estoy seguro al 100% de que esta sea la respuesta correcta, así que tómatelo con humor. En Stan las órdenes de declaración importa - tratar de poner el prior en lambda antes de usarlo. Es posible que lambda sea simplemente un muestreo de gamma(0.001, 0.001). O no definir una prioridad en absoluto, entonces usted consigue la distribución uniforme (o la versión impropia de eso) en el apoyo como una prioridad.