《基于贝叶斯理论的R语言实例分析报告》由会员分享，可在线阅读，更多相关《基于贝叶斯理论的R语言实例分析报告（13页珍藏版）》请在金锄头文库上搜索。
1、. . .上海大学20132014学年春季学期研究生课程考试课程名称：贝叶斯统计学课程编号：01SAQ9009论文题目:基于贝叶斯理论的R语言实例分析研究生姓名: 杨晓晓、李腾龙 学号:13720061、13720067研究生班级： 理学院统计系 论文评语:成 绩: 任课教师: 评阅日期: 基于贝叶斯理论的R语言实例分析 杨晓晓13720061,李腾龙13720067摘要：Gibbs抽样和Metropolis-Hastings算法是MCMC理论中最为重要的两种算法,Probit模型也是二分类数据分析中非常重要的模型。本文的主要是通过小组两个人互相讨论的方式,应用Gibbs和M-H算法共同完成了
2、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节
3、,我们知道潜变量,由于对做了如下限制条件：,这暗示潜变量的分布是以为条件的截尾正态分布truncated normal distribution,TN：再者,回归参数和潜变量为简单线性关系,由实用多元统计分析1第七章可知,所以：在先验分布的条件下,的后验分布为：也即,其中X为线性回归样本矩阵,第一列元素为1,Z为潜变量向量,最终得到回归参数和潜变量的满条件分布：2.2 Gibbs算法实现由2.1 我们得到了Probit模型的满条件分布,故Probit的Gibbs抽样可按下面过程执行：1选取合适的初始值；2从分布中抽取关于的样本,抽样过程中满足条件：；3从分布中抽取关于的样本；4重复过程2和3,
4、直到满足需要的样本量。2.3 程序实现：第一步,先给定参数,生成100个样本： #产生样本x-matrixy-vectorz-vectorb-cx-rnormx-cbindrep,xfor zi-xi,%*%b+rnorm if 0 yi=1 else yi=0 结果如下省略部分X值：第二步,用Gibbs抽样的方法,对参数进行估计：#Gibbs抽样library #截尾正态包library #多元正态包b1-1b2-1for for if zj-rnormTrunc1,xj,%*%c,1,min=0,max=Inf #截尾正态 else if zj-rnormTrunc1,xj,%*%c,1,
5、min=-Inf,max=0 b-mvrnorm1,solvet%*%x%*%t%*%z,solvet%*%x #多元正态 # solve求逆,t求转置b1i-b1b2i-b2summarysummary去掉前2000次收敛迭代过程,计算后3000次迭代值得到两个参数的后验均值估计结果如下大约需要运行2分零38秒：我们可以看到,估计结果E=1.9560,E=1.0730,与真实值比较接近。但经过多次运行发型,Gibbs方法产生的结果不太稳定,每次运行的结果差别较大,这可能与迭代次数有关。三、Probit模型与Metropolis-Hastings算法3.1 Probit模型在贝叶斯框架下的参数
6、后验分布由1.1,我们知道,的概率为：设参数服从先验分布：,根据贝叶斯原理,我们得到的后验分布为：3.2 M-H算法一般步骤一般的M-H算法是基于建议分布q 来产生Markov链的,令后验分布为平稳分布,由以下算法来产生Markov链：假设已由第t次迭代得到,为了得到,做以下步骤：1从一个建议分布取样：；2计算：3设定：3.3 用于Probit模型的M-H算法由3.1节,我们已经得到的后验分布,若取先验分布,得到：采用Metropolis-Hastings算法的迭代步骤为：任意选取,假设已经产生了,为了产生1从建议分布产生一个备选值则,概率密度函数 2计算3设定：3.4 程序实现对于2.3 G
7、ibbs算法中产生的同一组样本,我们实现M-H算法 #产生样本x-matrixy-vectorz-vectorb-cx-rnormx-cbindrep,xfor zi-xi,%*%b+rnorm if 0 yi-1 else yi-0 #M-H算法实现library #多元正态包b-c #取初始值phiz-vectorphib-vectorb1-vectorb2-vectorfor z-mvrnorm1,b,diag f-1 for phizj-pnorm #标准正态概率密度积分值 phibj-pnormf-f*yj*/yj* r-minu=runififu b=z else b=bb1i-b
8、1b2i-b2summarysummarypngparmfrow=cplotplotplotdev.off取初始值b=去掉前5000次收敛迭代过程,计算后5000次迭代值得到两参数的后验均值估计结果如下约运行20秒：b的真实值为迭代过程图如下：对于MH算法,我们产生的后验均值估计与真实值比较接近,并且从迭代过程图中我们也可以看到,算法收敛性质不错。四、改进的MH算法4.1 算法介绍在第三节中介绍的MH算法,我们采用了二元正态分布作为建议分布,其中直接采用了单位矩阵作为其协方差矩阵。而在文献10中,对M-H算法进行了改进,给出了建议分布协方差矩阵的另一种取法。在第t+1次迭代中,我们取建议分布为
9、：,其中4.2 程序实现 #产生样本x-matrixy-vectorz-vectorb-cx-rnormx-cbindrep,xfor zi-xi,%*%b+rnorm if 0 yi-1 else yi-0 #M-H算法实现library #多元正态包b-c #取初始值phiz-vectorphib-vectorb1-vectorb2-vectorfor a1-0 a2-0 a4-0# 建议分布协方差矩阵for phibt-pnorm a1-a1+yt/+/2*pi*2*exp-2+*/sqrt*phibt*exp-2/2 a2-a2+yt*xt,2/+*xt,2/2*pi*2*exp-2+xt,2*/sqrt*phibt*exp-2/2 a4-a4+yt*xt,22/+*xt,22/2*pi*2*exp-2+xt,22*/sqrt*phibt*exp-2/2 sigma-solvecbindc,c z-mvrnorm f-1 for phizj-pnorm #标准正态概率密度积分值 f-f*yj*/yj* r-minu-runififu b=z else b=b b1i-b1b2i-b2summaryb15000:10000
