文档介绍：评阅日期：
上海大学2013〜2014学年春季学期研究生课程考试
课程名称：贝叶斯统计学课程编号：_0〔SAQ9009,
论文题目:基于贝叶斯理论的R语言实例分析
研究生姓名：杨晓晓、李腾龙学号：13720061、1372006
0
0
1
1
0
1
Q
0
1
Q
1
0
1
0
1
1
0
1
1
0
0
1
0
1
0
1
Q
0
1
[73]
X
1
1
1
1
0
1
1
0
1
1
0
Q
1
1
1
0
-
0
0
0
1
1
1
0
1
0
1
1
01
1o
[
1
-5・72037102
[2rl
1
-4・07343339
[3.]
1
£・28623506
[£】
1
-6533
【乩］
1
-乳82013333
【巧］
1
7・34048503
[
1
[Sri
1
-1・3S405771
【即］
1
第二步，用Gibbs抽样的方法，对参数进行估计：
#Gibbs抽样
library(EnvStats)#截尾正态包
library(MASS)#多元正态包
b1<-1
b2<-1
for(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=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秒)：
aummazsy(LL:5000；)
--\.
。
swuaa^y(b2)
Min・151;Qj・\i・Ma.)t・
.-^
我们可以看到，估计结果E(bllY,X)=,E(b2IY,X)=,与真实值(2,1)比较接近。但经过多次运行发型，Gibbs方法产生的结果不太稳定，每次运行的结果差别较大,这可能与迭代次数有关。
三、Probit模型与Metropolis-Hastings算法
Probit模型在贝叶斯框架下的参数后验分布
，我们知道，y二1的概率为：P(y二1IX)二①(x'0)
iiii
所以：yIx〜BernoulliG,①(x'B))
iii
得到似然函数：
/(Y，…,YIx,…,x,P)二肝［①(x'B)］几［1—①(x'B)»-几
1n1nii
i=1
设参数p服从先验分布：p〜兀(p)，根据贝叶斯原理，我们得到B的后验
分布为：
兀(pIY,X)x/(Y,…，YIx,…，x,p)兀(p)
1n1n
x
i
i=1
njxip
乂
i=1
—82兀
P2
2Dp
Yi
1—Jx；p丄e
—82兀
2Dp
n［①(X’p)］Y1［1—0(X甲)］1-Y1兀(p)
M—H算法一般步骤
一般的M-H算法是基于建议分布q(yIx)来产生Markov链的，令后验分布为平稳分布，由以下算法来产生Markov链：
假设已由第t次迭代得到x(t)，为了得到x(t+1)，做以下步骤:
1)从一个建议分布取样：y~q(y|x(t))；t
2)计算：
r(x,y)二
min
f(y)q(x(t)|y)
tt
f(x(t))q(y|x(t))
t
,1>
(3)设定：x(t+1)=<儿
x(t)
以概率r
以概率1-r
-H算法
，我们已经得到0的后验分布，若取先验分布兀(P)