3

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

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

这段代码的作者已经定义了一组 coda 样本，如下所示：

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 矩阵。无论如何，是否可以为每个参数指定基本的统计属性？

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

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

4

1 回答 1

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

write.csv（gelman.diag（codaSamples）$psrf，文件=“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(). 使用此功能，您可以选择所需的所有子对象。

于 2014-06-12T14:24:00.137 回答