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

上海大学2013～2014学年春季学期研究生课程考试课程名称：贝叶斯统计学课程编号：01SAQ9009论文题目:基于贝叶斯理论的R语言实例分析研究生姓名:杨晓晓、李腾龙学号:13720061、13720067研究生班级：理学院统计系论文评语:成绩:任课教师:评阅日期:基于贝叶斯理论的R语言实例分析杨晓晓（13720061），李腾龙（13720067）摘要：Gibbs抽样和Metropolis-Hastings算法是MCMC理论中最为重要的两种算法，Probit模型也是二分类数据分析中非常重要的模型。本文的主要是通过小组两个人互相讨论的方式，应用Gibbs和M-H算法共同完成了Probit模型在贝叶斯理论框架下的估计问题，深入学习并掌握Probit模型、Gibbs抽样、M-H算法的相关知识，并能够初步使用R语言进行编程。同时，在文章第二部分我们俩还给出了多项式分布的Gibbs抽样的实现。关键词：Probit，Gibbs，Metropolis-Hastings，多项式分布一、Probit模型介绍1.1Probit模型的定义设y是一个二值的响应变量，10或iy。y的值依赖于解释变量x，通常我们可以认为1iy的概率是关于x的一个函数，即：)()|1(1iixfxyP假设存在潜在变量)1,(~)1,0(~iiiiiixNzNxz则：其中，0001iiizzy)()(1)(1)0()0()|1(iiiiiiiiixxxPxPzPxyP是参数，iz是潜变量。通常，我们称由上式决定的模型为Probit模型。二、Probit模型与Gibbs抽样2.1满条件分布由1.1节，我们知道潜变量)1,(~iiiiixNzxz：服从分布，由于对iz做了如下限制条件：0,0;1,0iiiiyzyz，这暗示潜变量iz的分布是以iy为条件的截尾正态分布（truncatednormaldistribution,TN）：)1,(~|iiixTNyz再者，iiixz，回归参数和潜变量iz为简单线性关系，由实用多元统计分析[1]第七章可知)1,(~iixNz，所以：npnpnpnniiiiinzzZxxxxxxXZXZXzxxyzzf11111111212)()(exp2)(exp),,|,,(其中，在先验分布1)(的条件下，的后验分布为：2)()(exp22exp2)()(exp),,|(11ZXXXXXZXXXZXXXZXZXxzyiii也即，11)(,)(~,,|XXZXXXNxzypiii，其中X为线性回归样本矩阵，第一列元素为1，Z为潜变量向量，最终得到回归参数和潜变量iz的满条件分布：1,~,|)(,)(~,|11iiipiixTNyzXXZXXXNyz2.2Gibbs算法实现由2.1我们得到了Probit模型的满条件分布，故Probit的Gibbs抽样可按下面过程执行：（1）选取合适的初始值)0(；（2）从分布1,ixTN中抽取关于iz的样本，抽样过程中满足条件：.0,0.0,1iiiizyzy；（3）从分布11)(,)(XXZXXXNp中抽取关于的样本；（4）重复过程（2）和（3），直到满足需要的样本量。2.3程序实现：第一步，先给定参数)1,2(，生成100个样本：#产生样本x-matrix()y-vector()z-vector()b-c(2,1)x-rnorm(100,0,5)x-cbind(rep(1,100),x)for(iin1:100){z[i]-x[i,]%*%b+rnorm(1,0,1)if(z[i]0){y[i]=1}else{y[i]=0}}结果如下（省略部分X值）：第二步，用Gibbs抽样的方法，对参数进行估计：#Gibbs抽样library(EnvStats)#截尾正态包library(MASS)#多元正态包b1-1b2-1for(iin2:5000){for(jin1:50){if(y[j]==1){z[j]-rnormTrunc(1,x[j,]%*%c(b1[i-1],b2[i-1]),1,min=0,max=Inf)#截尾正态}elseif(y[j]==0){z[j]-rnormTrunc(1,x[j,]%*%c(b1[i-1],b2[i-1]),1,min=-Inf,max=0)}}b-mvrnorm(1,solve(t(x)%*%x)%*%t(x)%*%z,solve(t(x)%*%x))#多元正态#solve()求逆,t()求转置b1[i]-b[1]b2[i]-b[2]}summary(b1[2000:5000])summary(b2[2000:5000])去掉前2000次收敛迭代过程，计算后3000次迭代值得到两个参数的后验均值估计结果如下（大约需要运行2分零38秒）：我们可以看到，估计结果E(b1|Y,X)=1.9560，E(b2|Y,X)=1.0730，与真实值(2,1)比较接近。但经过多次运行发型，Gibbs方法产生的结果不太稳定，每次运行的结果差别较大，这可能与迭代次数有关。三、Probit模型与Metropolis-Hastings算法3.1Probit模型在贝叶斯框架下的参数后验分布由1.1，我们知道，1iy的概率为：)()|1(iiixxyP)(,1~|iiixBernoullixy所以：niyiyinnxxxxyyf111111)](1[)]([),,,|,,(得到似然函数：设参数服从先验分布：)(~，根据贝叶斯原理，我们得到的后验分布为：niyxyxniyiyinniiiidedexxxxyyfXY11`221111)(21121)()](1[)]([)(),,,|,,(),|(22113.2M-H算法一般步骤一般的M-H算法是基于建议分布q(y|x)来产生Markov链的，令后验分布为平稳分布，由以下算法来产生Markov链：假设已由第t次迭代得到)(tx，为了得到)1(tx，做以下步骤：（1）从一个建议分布取样：)|(~)(ttxyqy；（2）计算：1,)|()()|()(min),()()()(ttttttxyqxfyxqyfyxr（3）设定：rxryxttt1)()1(以概率以概率3.3用于Probit模型的M-H算法由3.1节，我们已经得到的后验分布，若取先验分布1)(，得到：niyiyinnxxxxyyfXY111111)](1[)]([)(),,,|,,(),|(采用Metropolis-Hastings算法的迭代步骤为：任意选取),()0()0()0(ba，假设已经产生了),(,,),()()()0()0(ttbaba，为了产生),()1()1()1(tttba（1）从建议分布产生一个备选值IbaNzzzttt,),(~),()()(221则，概率密度函数2exp)|()|()()()()(ttttttttzzzqzq（2）计算1,)](1[)]([)](1[)]([min1,)|()()|()(min),(11)()(11)()()(niytiytiniytiytittttttiiiixxzxzxzqzqzyxr（3）设定：rrzttt1)()1(以概率以概率3.4程序实现对于2.3Gibbs算法中产生的同一组样本，我们实现M-H算法#产生样本x-matrix()y-vector()z-vector()b-c(2,1)x-rnorm(100,0,5)x-cbind(rep(1,100),x)for(iin1:100){z[i]-x[i,]%*%b+rnorm(1,0,1)if(z[i]0){y[i]-1}else{y[i]-0}}#M-H算法实现library(MASS)#多元正态包b-c(10,10)#取初始值phiz-vector()phib-vector()b1-vector()b2-vector()for(iin1:10000){z-mvrnorm(1,b,diag(1,2,2))f-1for(jin1:100){phiz[j]-pnorm(x[j,]%*%z)#标准正态概率密度积分值phib[j]-pnorm(x[j,]%*%b)f-f*((((phiz[j])^y[j])*((1-phiz[j])^(1-y[j])))/(((phib[j])^y[j])*((1-phib[j])^(1-y[j]))))}r-min(f,1)u=runif(1,0,1)if(ur){b=z}else{b=b}b1[i]-b[1]b2[i]-b[2]}summary(b1[5000:10000])summary(b2[5000:10000])png(4.png,height=480,width=1440)par(mfrow=c(1,3))plot(b1,b2,main='b1&b2',type='l',pch=20,col='blue')plot(b1,main='b1',type='l',pch=20,col='blue')plot(b2,main='b2',type='l',pch=20,col='blue')dev.off()取初始值b=(10,10)'去掉前5000次收敛迭代过程，计算后5000次迭代值得到两参数的后验均值估计结果如下（约运行20秒）：b的真实值为(2,1)'迭代过程图如下：对于MH算法，我们产生的后验均值估计与真实值比较接近，并且从迭代过程图中我们也可以看到，算法收敛性质不错。四、改进的MH算法4.1算法介绍在第三节中介绍的MH算法，我们采用了二元正态分布IbaNtt,),()()(2作为建议分布，其中直接采用了单位矩阵作为其协方差矩阵。而在文献[10]中，对M-H算法进行了改进，给出了建议分布协方差矩阵的另一种取法。在第t+1次迭代中，我们取建议分布为：)()()(2,),(tttbaN，其中nixbaittiixbaiiiinixbaittiixbaiiiinixbaittixbaiiitttittittittittittittexbayxexyxyaexbayxexyxyaaexbayeyyaxbaaaaa12)()()(2)(2222412)()()()(223112)()()()(221)()(14321)(2)()(2)()(2)()(2)()(2)()(2)()()1(2))(()1(2)1(2)1(2))(()1(2)1(2)1(2))(()1(212)(令：4.2程序实现#产生样本x-matrix()y-vector()z-