如何在 JAGS 中检查收敛

问题描述

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

解决方法

这个问题很可能是因为您正在尝试为每个参数计算 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,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."),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])
  gdag <- gelman.diag(tmp)[[1]]
  row.names(gdag) <- par_names[par]
  write.table(
    gdag,append = TRUE,row.names = TRUE
  )
}

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