了解 mcmc 包中使用的误差估计

Question

我在看the MCMC Example tutorial from the package mcmc in R。那里给出了下面的代码，其中 MCMC 用于计算逻辑回归示例中的参数。

library(mcmc) 

data(logit)
foo <- logit

out <- glm(y ~ x1 + x2 + x3 + x4, data = foo, family = binomial())

x <- foo
x$y <- NULL 
x <- as.matrix(x)
x <- cbind(1, x)
dimnames(x) <- NULL
y <- foo$y

lupost <- function(beta, x, y){
  eta <- x %*% beta
  p <- 1/(1 + exp(-eta))

  logl <- sum(log(p[y == 1])) + sum(log(1-p[y == 0]))

  return(logl + sum(dnorm(beta, 0, 2, log = TRUE)))
}

set.seed(42)
beta.init <- as.numeric(coefficients(out))
out <- metrop(lupost, beta.init, 1000, x=x, y = y)

这个想法是计算标准 Monte Carlo 错误。因此

> apply(out$batch, 2, mean)
[1] 0.6531950 0.7920342 1.1701075 0.5077331 0.7488265
[6] 0.5145751 0.7560775 1.4973807 0.3913837 0.7244162

已使用。我的问题是，此命令输出中最后 5 列的含义是什么？ 本教程指出它是某种 E(X^2)。但是在哪一行生成了X呢？这里的 X 是什么？

Answer 1

如果我运行只是你在上面问题中发布的代码，那么我只会得到 5 个数字：

apply(out$batch, 2, mean)
# [1] 0.5990174 0.8027482 1.1539329 0.4547702 0.7172724

您关注的 "tutorial" 似乎是 mcmc::demo 小插图。小插图执行您在上面发布的所有步骤，然后它还计算

out <- metrop(out, nbatch = 100, blen = 100, outfun = function(z) c(z, z^2))

接着是

apply(out$batch, 2, mean)
# [1] 0.6531950 0.7920342 1.1701075 0.5077331 0.7488265
# [6] 0.5145751 0.7560775 1.4973807 0.3913837 0.7244162

在这里您可以看到 outfun 正在计算 X 和 X^2，并且当您使用 [=15= 取平均值时，您可以得到 E[X] 和 E[X^2] 的估计值].

这里的 X 似乎是从模型参数（beta）的后验分布中得出的。请注意，贝叶斯后验均值（前 5 个数字）与代码第四行中使用 glm 计算的非贝叶斯点估计非常接近。

了解 mcmc 包中使用的误差估计

Understanding error estimation used in mcmc package

r

bayesian

mcmc