基于贝叶斯理论的R语言实例分析.docx

上海大学 2013 ~ 2014 学年春季学期研究生课程考试课程名称: 贝叶斯统计学课程编号: 01SAQ9009 论文题目: 基于贝叶斯理论的 R语言实例分析研究生姓名: 杨晓晓、李腾龙学号: 13720061 、13720067 研究生班级: 理学院统计系论文评语:成绩: 任课教师: 评阅日期: 基于贝叶斯理论的 R 语言实例分析杨晓晓( 13720061 ),李腾龙( 13720067 ) 摘要: Gibbs 抽样和 Metropolis-Hastings 算法是 MCMC 理论中最为重要的两种算法, Probit 模型也是二分类数据分析中非常重要的模型。本文的主要是通过小组两个人互相讨论的方式,应用 Gibbs 和 M-H 算法共同完成了 Probit 模型在贝叶斯理论框架下的估计问题, 深入学习并掌握 Probi t 模型、 Gibbs 抽样、 M-H 算法的相关知识, 并能够初步使用 R 语言进行编程。同时, 在文章第二部分我们俩还给出了多项式分布的 Gibbs 抽样的实现。关键词: Probit , Gibbs , Metropolis-Hastings ,多项式分布一、 Probit 模型介绍 Probit 模型的定义设y是一个二值的响应变量, 10或? iy 。y的值依赖于解释变量 x,通常我们可以认为 1? iy 的概率是关于 x的一个函数,即: )()|1( 1iixfxyP??假设存在潜在变量)1,(~)1,0(~???? iiiiiixNzNxz ????则: 其中, ??????00 01 i iiz zy)( )(1 )(1 )0()0()|1(?????? i i ii iiiiix x xP xPzPxyP?????????????????????是参数, iz 是潜变量。通常,我们称由上式决定的模型为 Probit 模型。二、 Probit 模型与 Gibbs 抽样 节,我们知道潜变量)1,(~??? iiiiixNzxz ????: 服从分布,由于对 iz 做了如下限制条件:0,0;1,0???? iiiiyzyz ,这暗示潜变量 iz 的分布是以 iy 为条件的截尾正态分布( truncated normal distribu tion,TN ):)1,(~|? iiix TN yz ?再者, iiixz?????,回归参数?和潜变量 iz 为简单线性关系,由实用多元统计分析[1] 第七章可知)1,(~? iixNz ?,所以: ?????????????????????????????????????????????????????????????????????npnp n pn ni iiiinz zZxx xxx xX ZXZX zxxyzzf????????? 11 1 1 11 1 1 2 12 )()( exp 2 )( exp ),,|,,(???????其中, 在先验分布 1)(???的条件下, ?的后验分布为: ????????????????????????????????????????????????????2 )()( exp 2 2 exp 2 )()( exp ),,|( 1 1ZXXXXXZXXX ZXXX ZXZXxzy iii?????????也即,?? 11)(,)(~,,| ?????XXZXXXNxzy piii?,其中 X为线性回归样本矩阵, 第一列元素为 1,Z为潜变量向量,最终得到回归参数?和潜变量 iz 的满条件分布: ?????????????1,~,| )(,)(~,| 11??? iii piix TN yz XXZXXXNyz Gibbs 算法实现由 我们得到了 Probit 模型的满条件分布,故 Probit 的 Gibbs 抽样可按下面过程执行: ( 1)选取合适的初始值)0(?; (2)从分布?? 1,? ix TN ?中抽取关于 iz 的样本,抽样过程中满足条件: .0,,1???? iiiizyzy ; ( 3)从分布?? 11)(,)( ?????XXZXXXN p中抽取关于?的样本; ( 4)重复过程( 2)和( 3),直到满足需要的样本量。 程序实现: 第一步,先给定参数)1,2( ???,生成 100 个样本: # 产生样本 x<-matrix() y<-vector() z<-vector() b<-c(2,1) x<-rnorm( 10 0,0,5) x<-cbind(rep(1, 10 0),x) for(i in 1: 10 0){ z[i]<-x[i,]%*%b+rnorm(1,0,1) if (z[i]>0){ y[i]=1 }els