bold<-1
clambda<-rep(0,(nsims+kstart))
cb<-rep(0,(nsims+kstart))
for(i in 1: (nsims+kstart))
{clambda[i]<-rgamma(1,shape=(x+1),scale=(bold/(bold+1)))
newy<-rgamma(1,shape=2,scale=(tau/(clambda[i]*tau+1)))
cb[i] <-1/newy
bold<-1/newy
}
gibbslambda<-clambda[(kstart+1):(nsims+kstart)]
gibbsb<-cb[(kstart+1):(nsims+kstart)]
out<-list(clambda=clambda, cb=cb,gibbslambda=gibbslambda,gibbsb= gibbsb)
return(out)
}
mean(hieracrh1(4000,6,0.04,1000)$gibbslambda)
mean(hieracrh1(4000,6,0.045,1000)$gibbslambda)
mean(hieracrh1(4000,6,0.05,1000)$gibbslambda)
mean(hieracrh1(4000,6,0.055,1000)$gibbslambda)
mean(hieracrh1(4000,6,0.06,1000)$gibbslambda)
