我想自动选择 MCMC 链的老化，例如通过基于收敛诊断删除前 n 行。
这一步可以在多大程度上安全地自动化？即使我仍然仔细检查自相关、mcmc 迹线和 pdf，也可以自动选择老化长度。
我的问题很笼统，但如果您能提供处理 R mcmc.object 的细节，那就太好了；我在 R 中使用 rjags 和 coda 包。
我可以半自动化 MCMC 收敛诊断来设置老化长度吗？
机器算法验证
r
贝叶斯
马尔可夫链蒙特卡罗
2022-03-11 22:57:13
1个回答
这是自动化的一种方法。非常感谢您的反馈。这是一种尝试用计算代替最初的目视检查，然后是随后的目视检查，以符合标准做法。
该解决方案实际上包含了两个可能的解决方案，首先，在达到某个阈值之前计算老化以去除链的长度，然后使用自相关矩阵计算细化间隔。
- 计算所有变量的最大中值 Gelman-Rubin 收敛诊断收缩因子 (grsf) 的向量
- 找到所有变量的 grsf 低于某个阈值的最小样本数，例如示例中的 1.1，在实践中可能更低
- 对从该点到链尾的链进行子采样
- 使用最自相关链的自相关来细化链
- 通过迹线、自相关和密度图直观地确认收敛
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 收缩因子的计算也是如此。这里有一个关于如何加快stackoverflow步骤4的问题
--更新第二部分--
附加答案：
无法诊断收敛性，只能诊断缺乏收敛性（Brooks、Giudici 和 Philippe，2003）
包runjags中的函数 autorun.jags自动计算运行长度和收敛诊断。在 Gelman rubin 诊断低于 1.05 之前，它不会开始监控链；它使用 Raftery 和 Lewis 诊断计算链长。
Gelman 等人（Gelman 2004 贝叶斯数据分析，第 295 页，Gelman 和 Shirley，2010 年）指出，他们使用一种保守的方法来丢弃链的前半部分。虽然是一个相对简单的解决方案，但实际上这足以解决我的特定模型和数据集的问题。
#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)