在 rjags 中使用 coda 样本进行诊断

Question

我对贝叶斯和 JAGS 都是新手，所以请原谅我的无知。

我收到了一个使用 JAGS 代码编写的 R 脚本（由一位同事）。

这段代码的作者定义了一组尾声样本如下：

codaSamples = coda.samples( jagsModel , variable.names=parameters , 
                            n.iter=nPerChain , thin=thinSteps )

我希望获得以下内容，但收效甚微：

Gelman 诊断：我使用了“show(gelman.diag(codaSamples))”，这对于单个模拟来说很好。但是，对于每个感兴趣的模拟，如何按参数将每个格尔曼诊断输出到文件？更有趣的是，是否可以仅记录 Rhat 值 >1.1 的模拟比例？

密度图：我使用了“show(densplot(codaSamples))”。但是，这会在单独的图上生成每个图（我在模型中有 96 个参数）。是否有“autocorr.plot”的等价物，每页放置几个图？

分位数：我使用了“show (summary(codaSamples))”，但是虽然这给出了每个参数的平均值、SD 和特定百分位数（这是我想要的），但它也给出了 MCMC 矩阵。有什么方法可以只为每个参数指定基本的统计属性吗？

后验分布：有没有办法计算每个参数的给定值（比如零）低于/高于的百分位数？然后总结所有模拟？

提前感谢您提供的任何帮助。

Answer 1

• 使用以下内容记录格尔曼输出：

write.csv(gelman.diag(codaSamples)$psrf, file = "gelman.csv")

• 绘制 mcmc 的密度：mcmc 对象具有特定的 S3 类 mcmc.list，因此 class(codaSamples) 应返回 mcmc.list。此类带有参数的对象有一个plot S3 方法：trace 和density。

绘图（codaSamples，迹线 = FALSE，密度 = TRUE）

这应该符合您的期望。

• summary(codaSamples) 给出参数的后验统计。不确定我是否理解这个问题，但 codaSamples 是 list 的长度：链数，列数：您的估计参数，行数：链的每次迭代。因此，您可以使用此对象或子对象（quantiles，...）计算您想要的每个统计信息。 例如第一个参数的分位数：

分位数(codaSamples[[1]][,1])

[[1]] 因为我采用第一个链和[,1] 作为第一个参数。

要查看对象的结构，您需要使用str()。使用此功能，您可以选择所需的所有子对象。