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

Question

我正在看the MCMC Example tutorial的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)

想法是计算标准蒙特卡罗误差。因此

> 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

看起来你所遵循的“教程”是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，并且当您使用apply()取平均值时，您可以得到E [X]和E [X ^ 2]的估计值。

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