Puedo semi-automatizar MCMC convergencia de diagnóstico para establecer el burn-in de longitud?

Question

Me gustaría automatizar la elección de burn-in por un MCMC de la cadena, por ejemplo, mediante la eliminación de las n primeras filas en base a una convergencia de diagnóstico.

¿Hasta qué punto puede este paso de forma segura automatizado? Incluso si todavía me doble verificación de la autocorrelación, mcmc de seguimiento, y de los pdfs, sería bueno tener la opción de grabar en longitud automatizado.

Mi pregunta es general, pero sería estupendo si pudiera proporcionarle detalles específicos para tratar con un R mcmc.objeto; estoy usando el rjags y coda paquetes en R.

Answer 1

Aquí es un enfoque en la automatización. Comentarios muy apreciado. Este es un intento de reemplazar inicial de la inspección visual con el cálculo, seguido por la posterior inspección visual, de acuerdo con la práctica estándar.

Esta solución incorpora dos posibles soluciones, en primer lugar, calcular burn-in para eliminar la longitud de la cadena antes de que algunos se alcanza el umbral, y, a continuación, utilizando la autocorrelación de la matriz para calcular el intervalo de adelgazamiento.

1. calcular un vector de la máxima mediana Gelman-Rubin convergencia de diagnóstico de reducir el factor (grsf) para todas las variables en el
2. encontrar el número mínimo de muestras en el que el grsf a través de todas las variables que van por debajo de un cierto umbral, por ejemplo, 1.1, en el ejemplo, tal vez menor en la práctica
3. sub-muestra las cadenas a partir de este punto hasta el final de la cadena
4. delgado de la cadena de uso de la autocorrelación de la mayoría de los autocorrelated de la cadena de
5. confirmar visualmente la convergencia con el seguimiento, autocorrelación y densidad de las parcelas

---

La mcmc objeto puede ser descargado aquí: entrecortado..Rdata

# jags.out is the mcmc.object with m variables
library(coda)    
load('jags.out.Rdata')
# 1. calculate max.gd.vec, 
# max.gd.vec is a vector of the maximum shrink factor
max.gd.vec     <- apply(gelman.plot(jags.out)$shrink[, ,'median'], 1, max)
# 2. will use window() to subsample the jags.out mcmc.object
# 3. start window at min(where max.gd.vec < 1.1, 100) 
window.start   <- max(100, min(as.numeric(names(which(max.gd.vec - 1.1 < 0)))))
jags.out.trunc <- window(jags.out, start = window.start)
# 4. calculate thinning interval
# thin.int is the chain thin interval
# step is very slow 
# 4.1 find n most autocorrelated variables
n = min(3, ncol(acm))
acm             <- autocorr.diag(jags.out.trunc)
acm.subset      <- colnames(acm)[rank(-colSums(acm))][1:n]
jags.out.subset <- jags.out.trunc[,acm.subset]
# 4.2 calculate the thinning interval
# ac.int is the time step interval for autocorrelation matrix
ac.int          <- 500 #set high to reduce computation time
thin.int        <- max(apply(acm2 < 0, 2, function(x) match(T,x)) * ac.int, 50)
# 4.3 thin the chain 
jags.out.thin   <- window(jags.out.trunc, thin = thin.int)
# 5. plots for visual diagnostics
plot(jags.out.thin)
autocorr.plot(jags.win.out.thin)

--update--

Como el implementado en R el cálculo de la autocorrelación de la matriz es más lento de lo que sería deseable (>15 min en algunos casos), en menor medida, por lo que es el cálculo de la GR reducir factor. Hay una pregunta acerca de cómo acelerar el paso 4 en stackoverflow aquí

--update--parte 2

otras respuestas:

1. No es posible diagnosticar la convergencia, sólo para diagnosticar la falta de convergencia (Brooks, Giudici, y Philippe, 2003)

2. La función de ejecución automática.exigencias del paquete runjags automatiza el cálculo de la longitud de ejecución y la convergencia de los diagnósticos. No iniciar la supervisión de la cadena hasta la Gelman rubin de diagnóstico está por debajo de 1.05; calcula la longitud de la cadena mediante el Raftery y Lewis diagnóstico.

3. Gelman et al (Gelman 2004 Bayesiano de Análisis de Datos, pág. 295, Gelman y Shirley, 2010) estado de que se utilice un enfoque conservador de descartar la 1ª mitad de la cadena. Aunque una solución relativamente simple, en la práctica, esto es suficiente para resolver el problema para mi en particular conjunto de modelos y de datos.

---

#code for answer 3
chain.length <- summary(jags.out)$end
jags.out.trunc <- window(jags.out, start = chain.length / 2)
# thin based on autocorrelation if < 50, otherwise ignore
acm <- autocorr.diag(jags.out.trunc, lags = c(1, 5, 10, 15, 25))
# require visual inspection, check acceptance rate
if (acm == 50) stop('check acceptance rate, inspect diagnostic figures') 
thin.int <- min(apply(acm2 < 0, 2, function(x) match(TRUE, x)), 50)
jags.out.thin <- window(jags.out.trunc, thin = thin.int)