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

大学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 模型介绍
Probit模型的定义
设y是一个二值的响应变量,。y的值依赖于解释变量x,通常我们可以认为的概率是关于x的一个函数,即:
假设存在潜在变量
是参数,是潜变量。
通常,我们称由上式决定的模型为Probit模型。
二、Probit模型与Gibbs抽样
满条件分布
,我们知道潜变量,由于对做了如下限制条件:,这暗示潜变量的分布是以为条件的截尾正态分布〔truncated normal distribution,TN〕:
再者,,回归参数和潜变量为简单线性关系,由实用多元统计分析[1]第七章可知,所以:
在先验分布的条件下,的后验分布为:
也即,,其中X为线性回归样本矩阵,第一列元素为1,Z为潜变量向量,最终得到回归参数和潜变量的满条件分布:
Gibbs算法实现
我们得到了Probit模型的满条件分布,故Probit的Gibbs抽样可按下面过程执行:
〔1〕选取适宜的初始值;
〔2〕从分布中抽取关于的样本,抽样过程中满足条件:;
〔3〕从分布中抽取关于的样本;
〔4〕重复过程〔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)=,E(b2|Y,X)=,与真实值(2,1)比拟接近。但经过屡次运行发型,Gibbs方法产生的结果不太稳定,每次运行的结果差异较大,这可能与迭代次数有关。
三、Probit模型与Metropolis-Hastings算法
Probit模型在贝叶斯框架下的参数后验分布
,我们知道,的概率为:
设参数服从先验分布:,根据贝叶斯原理,我们得到的后验分布为:
M-H算法一般步骤