tao和w是用rgamma()生成的样本,数字是不一样,是正常的。u、miu和beita是用rnorm()生成的,数字都一样。是rnorm()调用有问题,还是其他原因呢?我实在找不出原因了,求大神指教。代码如下:
y=c(0.1169,1.7311,0.0144,0.4881,2.3877,-0.0853,0.0522,0.7226,0.8718,3.0415,0.6363,1.359,1.5958,0.4567,
0.9885,1.7001,1.038,1.0195,0.3528,1.4249,0.317,1.3276,-0.1032,1.1568,0.0699,1.6802,0.247,0.4147,
1.6882,0.7256,1.3568,1.1091,0.05,2.2886,1.4985,2.7261,1.8443,-0.0687,0.9441,0.6872,0.5258,0.4743,
0.424,0.7349,0.4428,0.188,0.4642,0.2786,0.2742,0.5678)
Gibbs.sn=function(N,y,w.start,miu.start,beita.start,tao.start,b,d,burnin,thin=1){
w.draws=c()
miu.draws=c()
beita.draws=c()
tao.draws=c()
u.draws=matrix(NA,nrow = N-burnin,ncol = length(y))
w.cur=w.start
miu.cur=w.start
beita.cur=w.start
tao.cur=w.start
u.update=function(y,miu,beita,tao){
rnorm(length(y),(y-miu)*beita/((1/tao)+beita^2),sqrt((1/tao)/((1/tao)+beita^2)))
}
w.update=function(b,d,beita,tao){
rgamma(1,(d+1)/2,((beita^2)/((1/tao)*b)+d)/2)
}
miu.update=function(y,beita,u,tao){
rnorm(1,sum(y-beita*u)/length(y),sqrt((1/tao)/length(y)))
}
beita.update=function(b,y,miu,u,tao,w){
rnorm(1,sum((y-miu)*u)/((w/b)+sum(u^2)),sqrt((1/tao)/((w/b)+sum(u^2))))
}
tao.update=function(b,y,miu,beita,u,w){
rgamma(1,(length(y)+1)/2,(((w*beita^2)/b)+sum((y-miu)-beita*u)^2)/2)
}
for(i in 1:5000){
u.cur=u.update(y=y,miu=miu.cur,beita=beita.cur,tao=tao.cur)
w.cur=w.update(b=b,d=d,beita = beita.cur,tao = tao.cur)
miu.cur=miu.update(y=y,beita = beita.cur,u=u.cur,tao = tao.cur)
beita.cur=beita.update(b=b,y=y,miu=miu.cur,u=u.cur,tao=tao.cur,w=w.cur)
tao.cur=tao.update(b=b,y=y,miu=miu.cur,beita=beita.cur,u=u.cur,w=w.cur)
if(i>burnin&&(i-burnin)%%thin==0){
u.draws[(i-burnin)/thin,]=u.cur
w.draws[(i-burnin)/thin]=w.cur
miu.draws[(i-burnin)/thin]=miu.cur
beita.draws[(i-burnin)/thin]=beita.cur
tao.draws[(i-burnin)/thin]=tao.cur
}
}
return(list(u.draws=u.draws,w.draws= w.draws,miu.draws= miu.draws,beita.draws=beita.draws,tao.draws=tao.draws))
}
result=Gibbs.sn(N=5000,y=y,w.start=1,miu.start = 0.8,beita.start =1.2,tao.start = 26,b=1/2,d=2,burnin = 2000)