资源描述：
《基于贝叶斯理论的r语言实例分析》由会员上传分享，免费在线阅读，更多相关内容在行业资料-天天文库。
1、上海大学2013～2014学年春季学期研究生课程考试 课程名称：贝叶斯统计学课程编号：01SAQ9009论文题目:基于贝叶斯理论的R语言实例分析研究生姓名:杨晓晓、李腾龙学号:13720061、13720067研究生班级：理学院统计系论文评语: 成绩:任课教师: 评阅日期:基于贝叶斯理论的R语言实例分析杨晓晓（13720061），李腾龙（13720067）摘要：Gibbs抽样和Metropolis-Hastings算法是MCMC理论中最为重要的两种算法，Probit模型也是二分类数据分析中非常重要的模型。本文的主要是通过小组两个人互相讨论的方
2、式，应用Gibbs和M-H算法共同完成了Probit模型在贝叶斯理论框架下的估计问题，深入学习并掌握Probit模型、Gibbs抽样、M-H算法的相关知识，并能够初步使用R语言进行编程。同时，在文章第二部分我们俩还给出了多项式分布的Gibbs抽样的实现。关键词：Probit，Gibbs，Metropolis-Hastings，多项式分布一、Probit模型介绍1.1Probit模型的定义设y是一个二值的响应变量，。y的值依赖于解释变量x，通常我们可以认为的概率是关于x的一个函数，即：假设存在潜在变量是参数，是潜变量。通常，我们称由上式决定的模型为P
3、robit模型。二、Probit模型与Gibbs抽样2.1满条件分布由1.1节，我们知道潜变量，由于对做了如下限制条件：，这暗示潜变量的分布是以为条件的截尾正态分布（truncatednormaldistribution,TN）：再者，，回归参数和潜变量为简单线性关系，由实用多元统计分析[1]第七章可知，所以：在先验分布的条件下，的后验分布为：也即，，其中X为线性回归样本矩阵，第一列元素为1，Z为潜变量向量，最终得到回归参数和潜变量的满条件分布：2.2Gibbs算法实现由2.1我们得到了Probit模型的满条件分布，故Probit的Gibbs抽样可
4、按下面过程执行：（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(iin1:100){z[i]<-x[i,]%*%b+rnorm(1,0,1)if(z[i]>0){y[i]=1}else{y[i]=0}}结果如下
5、（省略部分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
6、=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
7、Y,X)=1.9560，E(b2
8、Y,X)=1.0730，与真实值(2,1)比较接近。但经过多次运行发型，
9、Gibbs方法产生的结果不太稳定，每次运行的结果差别较大，这可能与迭代次数有关。三、Probit模型与Metropolis-Hastings算法3.1Probit模型在贝叶斯框架下的参数后验分布由1.1，我们知道，的概率为：设参数服从先验分布：，根据贝叶斯原理，我们得到的后验分布为：3.2M-H算法一般步骤一般的M-H算法是基于建议分布q(y
10、x)来产生Markov链的，令后验分布为平稳分布，由以下算法来产生Markov链：假设已由第t次迭代得到，为了得到，做以下步骤：（1）从一个建议分布取样：；（2）计算：（3）设定：3.3用于Probit模型的
11、M-H算法由3.1节，我们已经得到的后验分布，若取先验分布，得到：采用Metropolis-Hastings算法的迭代步骤
