Gibbs（R）

Gibbs<-function(m,n,a,b){
  th<-matrix(0,ncol=2,nrow=m)
  theta0<-rbeta(1,a,b)
  x0<-rbinom(1,n,theta0)
  th[1,]<-c(theta0,x0)
  for(i in 2:m){
  	th[i,1]<-rbeta(1,th[i-1,2]+a,n-th[i-1,2]+b)
	th[i,2]<-rbinom(1,n,th[i,1])
      }
   list(x<-th[,2],theta<-th[,1])   
}
set.seed(123)
result<-Gibbs(20000,15,3,7)
result[[1]][1:10]
result[[2]][1:10]

hist(result[[1]],xlab="x",main="直方图")
d1<-density(result[[1]])
plot(d1,main="密度曲线")
hist(result[[2]],xlab="theta",main="直方图")
d2<-density(result[[2]])
plot(d2,main="密度曲线")