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

Question

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

我收到了一个使用 JAGS 代码编写的 R 脚本(由同事发送)。

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

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

我希望获得以下，并取得了有限的成功:

Gelman 诊断:我使用了“show(gelman.diag(codaSamples))”，它适用于单个模拟。但是，对于每个感兴趣的模拟，如何将每个参数的每个 gelman 诊断输出到一个文件中？更有趣的是，是否可以只记录 Rhat 值 >1.1 的模拟比例？

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

分位数:我使用了“show (summary(codaSamples))”，但尽管这给出了每个参数的平均值、SD 和特定百分位数(这正是我想要的)，但它也给出了 MCMC 矩阵。无论如何，是否可以只指定每个参数的基本统计属性？

后验分布:有没有办法计算每个参数的给定值(比如零)低于/高于的百分位数？然后对所有模拟进行总结？

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

Answer 1

• 用类似的东西记录 gelman 输出:

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

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

plot(codaSamples, trace = FALSE, density = TRUE)

这应该符合您的期望。

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

quantile(codaSamples[[1]][,1])

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

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