标签归档:参数估计

MCMC案例学习

本文是R中mcmc包的一篇帮助文档,作者为Charles J.Geyer。经过knitr编译后的pdf文档可见此处,提供中文译稿的作者:
闫超,天津财经大学统计系2011级研究生,方向:非寿险准备金评估。
高磊,天津财经大学统计系2011级研究生,方向:非寿险准备金评估。

这个案例,我们不关心题目的具体意义,重点放在利用贝叶斯的观点来解决问题时,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, 
    x = T)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-1.746  -0.691   0.154   0.704   2.194  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept)    0.633      0.301    2.10   0.0354 * 
x1             0.739      0.362    2.04   0.0410 * 
x2             1.114      0.363    3.07   0.0021 **
x3             0.478      0.354    1.35   0.1766   
x4             0.694      0.399    1.74   0.0817 . 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 137.628  on 99  degrees of freedom
Residual deviance:  87.668  on 95  degrees of freedom
AIC: 97.67

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 <- as.numeric(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])
    return(logl - sum(beta^2)/8)
}

继续阅读MCMC案例学习

闲扯“自由度”

        “闲扯”是一个四川方言词汇,指的就是大家在没事的时候坐下来吹吹牛,聊聊天。现在正是夜深人静的时候,找人聊聊天是不合适的,就由我一个人来自言自语下什么是自由度。

        我们进行统计分析,就像一个摄影师在拿着镜头在记录世界。但这个摄影师如果用的是广角镜头,那么他将面临一个问题:几何失真。特别是拍近景的时,拍出来的直线是弯曲的。这样就没有真失地反映客观事物的图像。所以这个时候他的反映真实客观现实的“自由”被限制了。虽然他的自由被限制了,但摄影师还是有办法矫正照的几何失真的:他可以尽量避免用广角镜头拍近景;他可以将照片交给专业的图像处理软件修复。所以,这个摄影师是有很多“自由”的手段来矫正照片失真的问题。这就可以当作是自由度的一个不恰当的类比。

继续阅读闲扯“自由度”