原创力文档- 1、本文档共18页，可阅读全部内容。
- 2、原创力文档（book118）网站文档一经付费（服务费），不意味着购买了该文档的版权，仅供个人/单位学习、研究之用，不得用于商业用途，未经授权，严禁复制、发行、汇编、翻译或者网络传播等，侵权必究。
- 3、本站所有内容均由合作方或网友上传，本站不对文档的完整性、权威性及其观点立场正确性做任何保证或承诺！文档内容仅供研究参考，付费前请自行鉴别。如您付费，意味着您自己接受本站规则且自行承担风险，本站不退款、不进行额外附加服务；查看《如何避免下载的几个坑》。如果您已付费下载过本站文档，您可以点击 这里二次下载。
- 4、如文档侵犯商业秘密、侵犯著作权、侵犯人身权等，请点击“版权申诉”（推荐），也可以打举报电话：400-050-0827(电话支持时间：9:00-18:30)。
查看更多
标准文档
实用文案
上海大学2013～2014学年 春 季学期研究生课程考试
?
课程名称： 贝叶斯统计学 课程编号： 01SAQ9009
论文题目: 基于贝叶斯理论的R语言实例分析
研究生姓名: 杨晓晓、李腾龙 学号:13720067
研究生班级： 理学院统计系
论文评语:
?
??
?
成 绩: 任课教师:
?
评阅日期:
基于贝叶斯理论的R语言实例分析
杨晓晓，李腾龙
摘要：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){
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<-1
b2<-1
for(i in 2:5000){
for(j in 1:50){
if (y[j]==1){
z[j]<-rnormTrunc(1,x[j,]%*%c(b1[i-1],b2[i-1]),1,min=0,max=Inf)
#截尾正态
}else if(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)比较接近。但经过多次运行发型，Gib
文档评论（0）