机器算法验证 - R2jags 有时不能部分去除烧伤？ - 吾爱随笔录

R2jags 有时不能部分去除烧伤？

机器算法验证 r 锯齿

2022-04-01 06:35:46

我发现即使使用选项，包中的功能jags()有时R2jags也不能部分消除刻录n.burnin=##。这是R（一个简单的线性模型）中的一个非常简单的示例：

library(R2jags)

N <- 1000
y <- rnorm(N)
x <- rnorm(N)
data <- list("N", "y", "x")

inits <- function(){list(beta0=rnorm(1), beta1=rnorm(1), tau=1)}
parameters <- c("beta0", "beta1", "tau")

模型m.bug是这样的：

model{
for (i in 1:N){
y[i] ~ dnorm(mu[i], tau)
mu[i] <- beta0 + beta1*x[i]
}

beta0 ~ dnorm(0, 0.00001)
beta1 ~ dnorm(0, 0.00001)
tau ~ dgamma(0.001, 0.001)
sigma2 <- 1/tau
}

在包中使用“jags()”，R2jags如下所示：

m <- jags(data, inits, parameters, "m.bug", 
  n.chains=3, n.iter=2000, n.burnin=1000, n.thin=1)

我的问题：

在的输出中m，后验估计基于正确的交互次数（1000），但如果我们检查跟踪图（使用traceplot(m)），老化部分似乎仍然存在（例如，“tau”的前几个值没有收敛） . 为什么？

并且只有一个进度条（通常两个，一个用于老化，一个用于其余）。

另外，如果我n.iter=2000, n.burnin=1000改为n.iter=2001000, n.burnin=2000000

经过的时间不会改变，这对于这么多迭代来说“太快了”。

附言。我使用R了 2.15.2 版和R2jags0.03-08 版

4个回答

我有完全相同的问题！在阅读了R2jags::jags源代码后，我得出了与@Hao Ye 相同的结论，我将详细展开：

首先，我注意到这只发生在某些模型上：

如果您dgamma从示例中删除并dunif改用，即如果您以这种方式修改示例：

...
inits <- function(){list(beta0=rnorm(1), beta1=rnorm(1), sigma2=1)}
parameters <- c("beta0", "beta1")
...
sigma2 ~ dunif(0, 100)
tau <- 1/sigma2
...

那么您将不会观察到错误

您观察到的错误也由具有
- 规则log(lambda) <- ...和组合dpois(lambda)
- 规则的组合log.lambda <- ...and dpois(exp(log.lambda))，但不是 for log.lambda <- ...and dpois(log.lambda)，而不是 for log(lambda) <- ...anddnorm(lambda, tau)

其次，（0.03-12版）中存在一个错误- 正如@Hao Ye注意到的那样，R2jags::jags它错误地将老化参数用作适配，但是对于那些不需要它的模型（见上文），适配阶段被跳过！参见R2jags::jags（0.03-12版）的代码：

if (n.burnin > 0) {
    n.adapt <- n.burnin
}
else {
    n.adapt <- 100
}
....
m <- jags.model(model.file, data = data, inits = init.values, 
    n.chains = n.chains, n.adapt = 0)
adapt(m, n.iter = n.adapt, by = refresh, progress.bar = progress.bar, 
    end.adaptation = TRUE)
samples <- coda.samples(model = m, variable.names = parameters.to.save, 
    n.iter = (n.iter - n.burnin), thin = n.thin, by = refresh, 
    progress.bar = progress.bar)
....

第二次调用 - 调用adapt- 只会在需要适应阶段的模型中执行某些操作（我不明白帮助中如何定义此行为?adapt- 可能是一些未定义的默认行为，还是由end.adaptation参数引起的？我不知道'不知道，文档好像不足。反正就是这样）。

无论如何，只有一个后续调用coda.samples意味着缺少真正的老化阶段。adapt运行而不是老化，并且仅适用于某些型号。

runjags由于这个错误，我开始使用包而不是 R2jags。

注意：如果有人决定修复此错误，则应以允许 JAGS 的默认适应阶段长度的方式完成 - 请参阅 runjags 包中的此问题：https ://stackoverflow.com/questions/22555421/ runjags-how-to-use-the-jags-default-for-the-length-of-adaptation-phase

我想我参加这个派对已经晚了一年......（我觉得这个问题可能属于stackoverflow，而不是）

问题似乎是 jags() 函数将 n.burnin 视为适应步骤的数量。并不是所有的模型都需要适配，这就是为什么宝月看到了一些模型的两个进度条，但这里的例子没有。

我不确定 jags() 是否具有真正老化的选项（即丢弃 MCMC 运行的初始部分），因此您可能必须填充 n.iter 并手动丢弃。或者，改用 rjags（由 JAGS 的作者编写）。

由于我在 JAGS 中进行了很多分析，并且更喜欢 R2jags 的输出格式，因此我决定一起破解这个问题的临时解决方案。

本质上，我在 R2jags 中的 jags() 函数中添加了一个参数 n.adapt，以指定“适应”迭代的次数，并去掉了设置 n.adapt = n.burnin 的代码块。然后在 R2jags 调用 jags.model() [适应发生的地方] 和 coda.samples() 从后验生成样本之间，我添加了对 update() 的调用，它更新模型的迭代次数等于jags() 调用中指定的老化时间。然后，无论 JAGS 决定如何处理适应， coda.samples() 的输出都不会包含来自老化期的那些样本。

我的解决方案似乎适用于问题中提供的玩具示例（即发生老化并且跟踪图看起来不错）

您可以在此处下载修改后的 R2jags 包。我已经在 Windows 和 Linux 上测试过它（不过只有 jags() 函数）

我不是构建 R 包的专家，但也许这对其他人有帮助。

我在自己的模型中也遇到过这个问题。一个解决方案（即，在修复此错误之前解决问题）：将先验从 tau 切换到 sigma。代替

tau~dgamma(0.001,0.001)
sigma2<-1/tau

尝试

sigma2 ~ dgamma(0.001, 0.001)
tau<-1/sigma2

或者

sigma2 ~ dgamma(0.001, 0.001)
tau <- 1/pow(sigma2,2)

这似乎解决了问题，并且 R2jags 不会忽略老化期。为什么会这样，我不知道。有人看到这种方法有什么问题吗？显然，在我的示例中，您需要为 sigma2 提供初始值，而不是 tau（即“sigma2 = 1”）。

其它你可能感兴趣的问题

上一篇什么是小型数据集建模的好技术？下一篇在卡方检验中不满足假设时该怎么办？