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

_上海大学20132014学年 春 季学期研究生课程考试课程名称： 贝叶斯统计学 课程编号： 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.1 Probit模型的定义设y是一个二值的响应变量，。y的值依赖于解释变量x，通常我们可以认为的概率是关于x的一个函数，即：假设存在潜在变量是参数，是潜变量。 通常，我们称由上式决定的模型为Probit模型。二、Probit模型与Gibbs抽样2.1 满条件分布由1.1节，我们知道潜变量，由于对做了如下限制条件：，这暗示潜变量的分布是以为条件的截尾正态分布（truncated normal distribution,TN）：再者，回归参数和潜变量为简单线性关系，由实用多元统计分析1第七章可知，所以：在先验分布的条件下，的后验分布为：也即，其中X为线性回归样本矩阵，第一列元素为1，Z为潜变量向量，最终得到回归参数和潜变量的满条件分布：2.2 Gibbs算法实现由2.1 我们得到了Probit模型的满条件分布，故Probit的Gibbs抽样可按下面过程执行：（1）选取合适的初始值；（2）从分布中抽取关于的样本，抽样过程中满足条件：；（3）从分布中抽取关于的样本；（4）重复过程（2）和（3），直到满足需要的样本量。2.3 程序实现：第一步，先给定参数，生成100个样本： #产生样本x-matrix()y-vector()z-vector()b-c(2,1)x-rnorm(100,0,5)x-cbind(rep(1,100),x)for(i in 1:100) zi0) yi=1 else yi=0 结果如下（省略部分X值）：第二步，用Gibbs抽样的方法，对参数进行估计： #Gibbs抽样library(EnvStats) #截尾正态包library(MASS) #多元正态包b1-1b2-1for(i in 2:5000) for(j in 1:50) if (yj=1) zj-rnormTrunc(1,xj,%*%c(b1i-1,b2i-1),1,min=0,max=Inf) #截尾正态 else if(yj=0) zj-rnormTrunc(1,xj,%*%c(b1i-1,b2i-1),1,min=-Inf,max=0) b-mvrnorm(1,solve(t(x)%*%x)%*%t(x)%*%z,solve(t(x)%*%x) #多元正态 # solve()求逆,t()求转置b1i-b1b2i-b2summary(b12000:5000)summary(b22000:5000)去掉前2000次收敛迭代过程，计算后3000次迭代值得到两个参数的后验均值估计结果如下（大约需要运行2分零38秒）：我们可以看到，估计结果E(b1|Y,X)=1.9560，E(b2|Y,X)=1.0730，与真实值(2,1)比较接近。但经过多次运行发型，Gibbs方法产生的结果不太稳定，每次运行的结果差别较大，这可能与迭代次数有关。三、Probit模型与Metropolis-Hastings算法3.1 Probit模型在贝叶斯框架下的参数后验分布由1.1，我们知道，的概率为：设参数服从先验分布：，根据贝叶斯原理，我们得到的后验分布为：3.2 M-H算法一般步骤一般的M-H算法是基于建议分布q(y|x) 来产生Markov链的，令后验分布为平稳分布，由以下算法来产生Markov链：假设已由第t次迭代得到，为了得到，做以下步骤：（1）从一个建议分布取样：；（2）计算：（3）设定：3.3 用于Probit模型的M-H算法由3.1节，我们已经得到的后验分布，若取先验分布，得到：采用Metropolis-Hastings算法的迭代步骤为：任意选取，假设已经产生了，为了产生（1）从建议分布产生一个备选值则，概率密度函数 （2）计算（3）设定：3.4 程序实现对于2.3 Gibbs算法中产生的同一组样本，我们实现M-H算法 #产生样本x-matrix()y-vector()z-vector()b-c(2,1)x-rnorm(100,0,5)x-cbind(rep(1,100),x)for(i in 1:100) zi0) yi-1 else yi-0 #M-H算法实现library(MASS) #多元正态包b-c(10,10) #取初始值phiz-vector()phib-vector()b1-vector()b2-vector()for(i in 1:10000) z-mvrnorm(1,b,diag(1,2,2) f-1 for(j in 1:100) phizj-pnorm(xj,%*%z) #标准正态概率密度积分值 phibj-pnorm(xj,%*%b)f-f*(phizj)yj)*(1-phizj)(1-yj)/(phibj)yj)*(1-phibj)(1-yj) r-min(f,1)u=runif(1,0,1)if(ur) b=z else b=bb1i-b1b2i-b2summary(b15000:10000)summary(b25000: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算法，我们采用了二元正态分布作为建议分布，其中直接采用了单位矩阵作为其协方差矩阵。而在文献10中，对M-H算法进行了改进，给出了建议分布协方差矩阵的另一种取法。在第t+1次迭代中，我们取建议分布为：，其中4.2 程序实现 #产生样本x-matrix()y-vector()z-vector()b-c(2,1)x-rnorm(50,0,1)x-cbind(rep(1,50),x)for(i in 1:50) zi0) yi-1 else yi-0 #M-H算法实现library(MASS) #多元正态包b-c(0,0) #取初始值phiz-vector()phib-vector()b1-vector()b2-vector()for(i in 1:10000) a1-0 a2-0 a4-0# 建议分布协方差矩阵 for(t in 1:50) phibt-pnorm(xt,%*%b) a1-a1+(yt/(2*pi*phibt2)+(1-yt)/(2*pi*(1-phibt)2)*exp(-(xt,%*%b)2)+(yt-phibt)*(xt,%*%b)/(sqrt(2*pi)*phibt*(1-phibt)*exp(-(xt,%*%b)2/2) a2-a2+(yt*xt,2/(2*pi*phibt2)+(1-yt)*xt,2/(2*pi*(1-phibt)2)*exp(-(xt,%*%b)2)+(xt,2*(yi-phibt)*(xt,%*%b)/(sqrt(2*pi)*phibt*(1-phibt)*exp(-(xt,%*%b)2/2) a4-a4+(yt*xt,22/(2*pi*phibt2)+(1-yt)*xt,22/(2*pi*(1-phibt)2)*exp(-(xt,%*%b)2)+(xt,22*(yt-phibt)*(xt,%*%b)/(sqrt(2*pi)*phibt*(1-phibt)*exp(-(xt,%*%b)2/2) sigma-solve(cbind(c(a1,a2),c(a2,a4) z-mvrnorm(1,b,sigma) f-1 for(j in 1:50) phizj-pnorm(xj,%*%z) #标准正态概率密度积分值 f-f*(phizj)yj)*(1-phizj)(1-yj)/(phibj)yj)*(1-phibj)(1-yj) r-min(f,1)u-runif(1,0,1)if(ur) b=z else b=b b1i-b1b2i-b2summary(b15000:10000)summary(b25000: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()4.3 运行结果很遗憾，因为算法本身存，或者我们程序编写存在问题，不管初始值如何变换，在迭代的过程中，始终会出现非正定的情况，使得迭代中断。错误提醒：）但即便如此，我们也产生了29个b的值如下：迭代过程图如下：从图中我们也可以看到，在迭代过程中，b1和b2也是朝着真实值在慢慢接近，这也可初步说明理论或许是可行的，但在实际操作中由于是数值计算求逆，从而产生了一些问题。五、Gibbs sampling (Multinomial posterior)5.1 多项式分布考虑一个多项式分布：用Gibbs抽样方法进行Bayes估计：可以将原模型看成：在这个条件下，的后验分布可以容易得到：更多的，考察原多项式分布，第一种情况共出现了次，其中出现了次，出现了，而出现的概率为。所以，可求得的条件概率密度为：，同理：5.2 Gibbs 抽样5.2.1 采用Gibbs抽样的迭代步骤为：给定初始值，假设经过前t步迭代，已经产生了第t+1步迭代：1. 2. 5.2.2 程序实现：第一步，产生样本：对于模型，我们给定真实值，a=(0.2,0.4,0.1,0.5)，b=(0.05,0.1,0.1,0.15)，则根据建议，可以取R程序如下：y-rmultinom(1, size=800, prob=c(0.12,0.24,0.125,0.275,0.24)结果：y=( 94,216,82,224,184)第二步，进行Gibbs抽样：R程序如下：a-c(0.2,0.4,0.1,0.5)b-c(0.05,0.1,0.1,0.15)y-c(94,216,82,224,184)alpha-c(0.5,0.5,0.5)z-c(1,1,1,1)mu-vector()eta-vector()library(BGSIMD) #Dirichlet函数包for(i in 1:10000)mueta-rdirichlet(1,c(z1+z2+alpha1,z3+z4+alpha2,y5+alpha3) z1-rbinom(1,y1,(a1*mueta1/(a1*mueta1+b1) z2-rbinom(1,y2,(a2*mueta1/(a2*mueta1+b2) z3-rbinom(1,y3,(a3*mueta2/(a3*mueta2+b3) z4-rbinom(1,y4,(a4*mueta2/(a4*mueta2+b4) mui-mueta1 etai-mueta2summary(mu2000:10000)summary(eta2000:10000)plot(mu,eta,pch=20,col=blue)5.2.3 结果从结果中我们看到，Gibbs估计的精确度并不非常高，并且经过实验，每次估计所得到的结果相差较大。六、总结从上面几个例子中，我们可以看到，使用Gibbs抽样和M-H算法最终得到的结果，有时并不十分精确，并且每次结果的差异也较大。并且文献中改进的MH算法，迭代过程也未能完全实现。我们后续研究的方向，首先还是要把理论证明弄清楚，这样才能容易找到问题的原因。其次，或许会尝试用广义逆的理论，解决在改进MH算法中遇到的问题。参考文献1 Richard A.Johnson/Dean W.Wichern. Applied Multivariate Statistical Analysis. Pearson Prentice Hall.20082 Christian P.Robert. The Bayesian Choice. New York. Springer Verlag, 20073