如何检查 JAGS 中的收敛性

Question

我已经在 JAGS 中为数百个人的时间序列估计了一个巨大而复杂的时变层次模型。我估计个人的时变参数导致每人 80 个参数。参数是服从多元正态分布的随机效应。当我尝试使用 coda::gelman.diag(samples) 检查各个参数的收敛性时，我的计算机运行出现内存问题，R 在 12 小时后崩溃。是否足以检查多元正态分布的均值和精度的收敛性，或者是否有必要检查个体参数的收敛性（遵循上述分布）？这里的约定是什么？

Answer 1

这个问题很可能是因为您正在尝试为每个参数计算 gelman rubin 诊断，结果是运行 RAM 不足。对此有一些变通办法，即为每个参数一次应用一个函数。我假设您有 RAM 来将结果存储在矩阵中，因为您的计算机可以存储模型的“样本”。

这是 rjags 中的一个非常简单的示例模型，以及您如何一次处理一个参数。

library(rjags)

# load the model
data(LINE)

# recompile it
LINE$recompile()

# get some samples
LINE.out <- coda.samples(LINE, c("alpha","beta","sigma"), n.iter=1000)

# Get the number of parameters
npar <- ncol(LINE.out[[1]])

# and the names of the parameters
par_names <- colnames(LINE.out[[1]])

# A matrix to store all of the results
gdag <- matrix(NA, ncol = 2, nrow = npar)

# give rows informative names
row.names(gdag) <- par_names

# give the columns informative names (based on output of gelman.diag())
colnames(gdag) <- c("PointEst", "UpperC.I.")

# progress bar
pb <- txtProgressBar(0, npar)

# for loop to work through things one at a time
for(par in 1:npar){
  setTxtProgressBar(pb, par)
  tmp <- lapply(LINE.out, function(x) x[,par])
  gdag[par,] <- gelman.diag(tmp)[[1]]
}

如果出于某种原因，您无法将所有内容存储在 gdag 中，您可以使用 write.table() 将结果附加到 csv。像这样的东西应该有用。

library(rjags)

# load the model
data(LINE)

# recompile it
LINE$recompile()

# get some samples
LINE.out <- coda.samples(LINE, c("alpha","beta","sigma"), n.iter=1000)

# Get the number of parameters
npar <- ncol(LINE.out[[1]])

# and the names of the parameters
par_names <- colnames(LINE.out[[1]])

# Create the headers of the table
write.table(
  matrix(c("PointEst", "UpperC.I."), ncol = 2, nrow = 1),
  "my_gelman_rubins.csv",
  sep = ",",
  col.names = FALSE,
  row.names = TRUE
)

# progress bar
pb <- txtProgressBar(0, npar)

for(par in 1:npar){
  setTxtProgressBar(pb, par)
  tmp <- lapply(LINE.out, function(x) x[,par])
  gdag <- gelman.diag(tmp)[[1]]
  row.names(gdag) <- par_names[par]
  write.table(
    gdag,
    "my_gelman_rubins.csv",
    append = TRUE,
    sep = ",",
    col.names = FALSE,
    row.names = TRUE
  )
}

后一种方法会花费更长的时间，但还有一个额外的好处，即可以将结果实际保存到文件中。因此，如果你确实有内存分配问题，你可以在你需要开始的任何参数上恢复你的诊断计算。例如，如果您在第 501 个参数上遇到错误，您只需将 for 循环修改为 for(par in 501:npar)

如何检查 JAGS 中的收敛性

How to check for convergence in JAGS

jags

r2jags

convergence

rjags