0

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

4

1 回答 1

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(). 像这样的东西应该工作。

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)

于 2021-03-31T13:24:48.320 回答