¿R2jags no quita la quemadura en parte a veces?

Question

Me enteré de que la función jags() en la R2jags paquete a veces no quitar la quemadura en la parte incluso con la opción n.burnin=##. Aquí es un ejemplo muy simple en R (un simple modelo lineal):

library(R2jags)

N <- 1000
y <- rnorm(N)
x <- rnorm(N)
data <- list("N", "y", "x")

inits <- function(){list(beta0=rnorm(1), beta1=rnorm(1), tau=1)}
parameters <- c("beta0", "beta1", "tau")

el modelo m.bug es como este:

model{
for (i in 1:N){
y[i] ~ dnorm(mu[i], tau)
mu[i] <- beta0 + beta1*x[i]
}

beta0 ~ dnorm(0, 0.00001)
beta1 ~ dnorm(0, 0.00001)
tau ~ dgamma(0.001, 0.001)
sigma2 <- 1/tau
}

el uso de "puntas()" R2jags paquete como este:

m <- jags(data, inits, parameters, "m.bug", 
  n.chains=3, n.iter=2000, n.burnin=1000, n.thin=1)

Mi pregunta:

En la salida de m, la parte posterior de las estimaciones se basan en el número de la derecha de interations (1000), pero si revisamos las traceplot (usando traceplot(m)), el burnin parte parece que todavía existe(por ejemplo, el primer par de valores de "tau" no son convergentes). Por qué?

y sólo hay una barra de progreso (uaually dos, uno para la quemadura, uno para el resto).

también, si yo cambio n.iter=2000, n.burnin=1000 a n.iter=2001000, n.burnin=2000000

el tiempo transcurrido no cambia, lo que es "demasiado rápido" para tantas iteraciones.

ps. Yo solía R versión 2.15.2 y R2jags versión 0.03-08

Answer 1

Yo tenía exactamente el mismo problema! Después de la lectura de la R2jags::jags fuente llegué a la misma conclusión que @Hao Ye, voy a ampliar en más detalle:

En primer lugar, me di cuenta de que esto sucede sólo para algunos modelos:

• si se quita la dgamma de su ejemplo y de uso dunif en su lugar, es decir, si usted modifica su ejemplo de esta manera:

  ...
  inits <- function(){list(beta0=rnorm(1), beta1=rnorm(1), sigma2=1)}
  parameters <- c("beta0", "beta1")
  ...
  sigma2 ~ dunif(0, 100)
  tau <- 1/sigma2
  ...
  

  entonces no va a observar el error

• Error que se observa es también producida por los modelos con

  • combinación de reglas log(lambda) <- ... y dpois(lambda)
  • combinación de reglas log.lambda <- ... y dpois(exp(log.lambda)), pero no por log.lambda <- ... y dpois(log.lambda), y no de log(lambda) <- ... y dnorm(lambda, tau)

Segundo, hay un error en R2jags::jags (de la versión 0.03-12) - como @Hao Habéis dado cuenta, que por error toma de burn-in parámetro y la utiliza como adaptación, pero para aquellos modelos que no necesitan (vea arriba), la fase de adaptación se omite! Ver el código de R2jags::jags (de la versión 0.03-12):

if (n.burnin > 0) {
    n.adapt <- n.burnin
}
else {
    n.adapt <- 100
}
....
m <- jags.model(model.file, data = data, inits = init.values, 
    n.chains = n.chains, n.adapt = 0)
adapt(m, n.iter = n.adapt, by = refresh, progress.bar = progress.bar, 
    end.adaptation = TRUE)
samples <- coda.samples(model = m, variable.names = parameters.to.save, 
    n.iter = (n.iter - n.burnin), thin = n.thin, by = refresh, 
    progress.bar = progress.bar)
....

La segunda llamada - llamada a adapt - sólo lo va a hacer algo en los modelos que necesita la fase de adaptación (no entiendo cómo este comportamiento se define en la ?adapt ayuda - posiblemente algunos indefinido comportamiento por defecto, o es causada por la end.adaptation parámetro? No sé, la documentación parece insuficiente. De todos modos, es como este).

De todos modos, hay sólo una llamada posterior a coda.samples lo que significa que la grabación real-en la segunda fase es la que falta. El adapt se ejecuta en lugar de burn-in, y sólo para algunos modelos.

Yo empecé a usar el paquete runjags en lugar de R2jags porque de este error.

---

Nota: si alguien se decide a corregir este error debe ser hecho de tal manera que permite a las EXIGENCIAS' de longitud por defecto de la fase de adaptación - ver este problema en runjags paquete: http://stackoverflow.com/questions/22555421/runjags-how-to-use-the-jags-default-for-the-length-of-adaptation-phase

Answer 2

Creo que soy un año es demasiado tarde para este partido... (y siento que esta pregunta puede estar en stackoverflow, en lugar)

El problema parece ser que las puntas() función trata n.burnin como el número de la adaptación de los pasos. No todos los modelos requieren la adaptación, que es la razón por la Baoyue ha visto dos barras de progreso para algunos modelos, pero no para el ejemplo aquí.

No estoy seguro de si entrecortado() tiene opciones para la verdadera burn-in (es decir, descartar la porción inicial de la MCMC de ejecución), lo que podría tener para el cojín n.el iter y el manual de descarte. Como alternativa, utilice rjags (escrito por el autor de PUNTAS) en su lugar.

Answer 3

He experimentado este problema en mis propios modelos. Una solución (es decir, evitar hasta que se corrija este error): interruptor de la previa de la tau a sigma. En lugar de

tau~dgamma(0.001,0.001)
sigma2<-1/tau

trate de

sigma2 ~ dgamma(0.001, 0.001)
tau<-1/sigma2

o

sigma2 ~ dgamma(0.001, 0.001)
tau <- 1/pow(sigma2,2)

Esto parece solucionar el problema, y R2jags NO ignorará la quemadura en el período. Por qué esto funciona, no sé. Que nadie vea algún problema con este enfoque? Obviamente en mi ejemplo, usted necesita proporcionar un valor inicial para sigma2, no tau (es decir, "sigma2 = 1").

Answer 4

Ya tengo un montón de análisis en PUNTAS y prefieren el formato de salida de R2jags me decidí a hackear una solución temporal a este problema.

Básicamente he añadido un argumento a la entrecortado() en función de R2jags, n.adaptar, para especificar el número de 'adaptación' iteraciones y se deshizo de el fragmento de código que establece n.adaptar = n.burnin. A continuación, entre el R2jags llamadas a entrecortado.modelo() [si la adaptación se produce] y coda.muestras (), que genera las muestras de la parte posterior, he añadido una llamada a update() que actualiza el modelo para un número de iteraciones igual a la burnin período especificado en las puntas() la llamada. Entonces, la salida de coda.muestras() no contienen las muestras de la burnin período, independientemente de cómo ENTRECORTADO decide manejar la adaptación.

Mi solución que parece funcionar para el juguete de ejemplo proporcionado en la pregunta (es decir, el burn-in se produce y el traceplots ven bien)

Puede descargar la versión modificada de la R2jags paquete aquí. Lo he probado en Windows y Linux (sólo las puntas() función, aunque)

No soy un experto en la creación de paquetes de R, pero tal vez esto va a ser útiles a los demás.

Answer 5

Seguimiento a Jason Hill respuesta, mientras que esto soluciona el problema, no debe hacerlo. Si $\tau$, la precisión, es la inversa de a $\sigma^2$ (la varianza), y $\tau$ es distribuido gamma, entonces se sigue que $\sigma^2$ no será distribuido gamma, pero a la inversa gamma. No creo que ENTRECORTADO admite que la parametrización, que es la razón por la que tantas personas usando $\tau$ en sus modelos.

Wikipedia tiene un resumen de la distribución, si usted está interesado.