Могу ли я полуавтоматизировать диагностику сходимости MCMC для установки длины выгорания?

Я хотел бы автоматизировать выбор выгорания для цепочки MCMC, например, удалив первые n строк на основе диагностики сходимости.

В какой степени этот шаг можно безопасно автоматизировать? Даже если я все еще дважды проверю автокорреляцию, трассировку mcmc и pdf, было бы неплохо автоматизировать выбор длины записи.

Мой вопрос общий, но было бы замечательно, если бы вы могли предоставить конкретные сведения для работы с R mcmc.object; Я использую пакеты rjags и coda в R.

Вот один из подходов к автоматизации. Обратная связь высоко ценится. Это попытка заменить первоначальный визуальный осмотр вычислением с последующим последующим визуальным осмотром в соответствии со стандартной практикой.

Это решение на самом деле включает в себя два потенциальных решения: сначала вычислить выгорание, чтобы удалить длину цепи до того, как будет достигнут некоторый порог, а затем использовать автокорреляционную матрицу для вычисления интервала утонения.

1. рассчитать вектор максимального медианного диагностического коэффициента усадки Гельмана-Рубина (grsf) для всех переменных в
2. найти минимальное количество выборок, при котором grsf по всем переменным опускается ниже некоторого порога, например, 1,1 в примере, возможно, на практике ниже
3. суб образец цепей от этой точки до конца цепочки
4. тонкая цепь с помощью автокорреляции наиболее автокоррелированной цепи
5. визуально подтвердить сходимость на графиках трассировки, автокорреляции и плотности

Объект mcmc можно скачать здесь: jags.out.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)

--Обновить--

Как реализовано в R, вычисление матрицы автокорреляции происходит медленнее, чем хотелось бы (> 15 минут в некоторых случаях), в меньшей степени, как и вычисление коэффициента сжатия GR. Есть вопрос о том, как ускорить шаг 4 на стеке потока здесь

- обновить часть 2--

дополнительные ответы:

1. Невозможно диагностировать конвергенцию, только диагностировать отсутствие конвергенции (Brooks, Giudici, Philippe, 2003)

2. Функция autorun.jags из пакета runjags автоматизирует расчет длины пробега и диагностику сходимости. Он не начинает мониторинг цепи, пока диагностика рубина Гельмана не станет ниже 1,05; он вычисляет длину цепи с использованием диагностики Raftery и Lewis.

3. Гельман и др. (Gelman 2004 Bayesian Data Analysis, p. 295, Gelman and Shirley, 2010 ) утверждают, что они используют консервативный подход к отбрасыванию 1-й половины цепочки. Хотя это относительно простое решение, на практике этого достаточно, чтобы решить проблему для моего конкретного набора моделей и данных.

#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)