文档介绍：MCMC 案例学习
Charles
这个案例,我们不关心题目的具体意义,重点放在利用贝叶斯的观点来解决问题
时,MCMC 在后续的计算中所发挥的巨大作用。我们知道,贝叶斯的结果往往是一
个后验分布。这个后验分布往往很复杂,我们难以用经典的方法求解其期望与方差等
一系列的数据特征,这时 MCMC 来了,将这一系列问题通过模拟来解决。从这个意
义上说,MCMC 是一种计算手段。依频率学派看来,题目利用广义线性模型可以解
决,在贝叶斯看来同样可以解决,但是遇到了一个问题,就是我们得到的非标准后验
分布很复杂。我们正是利用 MCMC 来解决了这个分布的处理问题。本文的重点也在
于此。
在使用 MCMC 时作者遵循了这样的思路,首先依照贝叶斯解决问题的套路,构
建了非标准后验分布函数。然后初步运行 MCMC,确定合适的 scale。继而,确定适
当的模拟批次和每批长度(以克服模拟取样的相关性)。最后,估计参数并利用 delta
方法估计标准误。
1 问题的提出
这是一个关于 R 软件中 mcmc 包的应用案例。问题出自明尼苏达大学统计系博
士入学考试试题。这个问题所需要的数据存放在 logit 数据集中。在这个数据集中
有五个变量,其中四个自变量 x1、x2、x3、x4,一个响应变量 y。
对于这个问题,频率学派的处理方法是利用广义线性模型进行参数估计,下面是
相应的 R 代码以及结果:
> library(mcmc)
> data(logit)
> out <- glm(y ~ x1 + x2 + x3 + x4, data = logit, family = binomial(), x = T)
> summary(out)
Call:
glm(formula = y ~ x1 + x2 + x3 + x4, family = binomial(), data = logit,
1
1 问题的提出 2
x = T)
Deviance Residuals:
Min 1Q Median 3Q Max
- -
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) *
x1 *
x2 **
x3
x4 .
---
Signif. codes: 0 '***' '**' '*' '.' ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: on 99 degrees of freedom
Residual deviance: on 95 degrees of freedom
AIC:
Number of Fisher Scoring iterations: 6
但是,使用频率学派的分析方法解决这个问题并不是我们想要的,我们希望使用
Bayesian 分析方法。对于 Bayesian 分析而言,我们假定数据模型(即广义线性模型)
与频率学派一致。同时假定,五个参数(回归系数) 相互独立,并服从均值为 0、标
准差为 2 的先验正态分布。
定义下面的 R 函数来计算非标准的对数后验分布概率密度(先验密度与似然函
数相乘)。我们为什么要定义成密度函数的对数形式?因为虽然我们是从分布中采样,
但是 MCMC 算法的执行函数 metrop() 需要的一个参数正是这个分布的密度函数的
对数形式。
> x <- out$x
> y <- out$y
> lupost <- function(beta, x, y) {
+ eta <- (x %*% beta)
+ logp <- ifelse(eta < 0, eta - log1p(exp(eta)), -log1p(exp(-eta)))
+ logq <- ifelse(eta < 0, -log1p(exp(eta)), -eta - log1p(exp(-eta)))
+ logl <- sum(logp[y == 1]) + sum(logq[y == 0])
2 开始