如何在R2OpenBugs中获得蒙特卡罗错误?

Question

是否有人在运行R2OpenBugs的贝叶斯模型时设法获得一个参数的蒙特卡罗错误?

在OpenBugs的标准输出中提供，但是在R2OpenBugs下运行时，日志文件中没有MC error。是否有办法要求R2OpenBugs计算MC错误?或者有一种人工计算的方法?如果你有什么办法，请告诉我。谢谢你！

以下是R2OpenBugs的标准日志输出:

$stats
mean      sd val2.5pc    median val97.5pc sample
beta0      1.04700 0.13250   0.8130   1.03800   1.30500   1500
beta1     -0.31440 0.18850  -0.6776  -0.31890   0.03473   1500
beta2     -0.05437 0.05369  -0.1648  -0.05408   0.04838   1500
deviance 588.70000 7.87600 575.3000 587.50000 606.90000   1500
$DIC
Dbar  Dhat   DIC    pD
t     588.7 570.9 606.5 17.78
total 588.7 570.9 606.5 17.78

Answer 1

计算蒙特卡罗标准误差(MCSE)的一种简单方法是将链的标准差除以有效样本数的平方根。在您的输出中提供了标准偏差，但是当您打印模型输出时，有效样本大小应该以n.eff(最右边的列)给出—或者至少这是我从以下内容中得到的印象:

https://cran.r-project.org/web/packages/R2OpenBUGS/vignettes/R2OpenBUGS.pdf

我不再使用OpenBugs了，所以不能轻易地为你检查，但是应该有一些东西表明有效的样本大小(这与你采样的迭代次数不同，因为它也考虑到了由于链内相关性而导致的信息损失)。

否则，您可以通过提取原始MCMC链来获得它，然后使用coda包(?coda:: effecvesize)计算有效样本量，或者直接使用LaplacesDemon::MCSE来计算蒙特卡洛标准误差。有关更多信息，请参阅:

https://rdrr.io/cran/LaplacesDemon/man/MCSE.html

指出,一些人(包括我)建议关注有效样本量直接,而不是看着MCSE老"thumb&quot规则;MCSE应小于样本标准差的5%相当于说有效样本量应至少为400(1/0.05^2)。但意见不一:)

Answer 2

MCMC-error被命名为Time-series SE，可以在coda对象摘要的统计部分中找到:

library(R2OpenBUGS)
library(coda)
my_result <- bugs(...., codaPg = TRUE)
my_coda <- read.bugs(my_result)
summary(my_coda$statistics)