gibbs=function (n=10000,k=10, mean=0) { mat=matrix(ncol = k+1, nrow = n) x=rep(1,k) y=1 mat[1, ]=c(x, y) for (i in 2:n) { x=rnorm(n=k, mean=0, sd=sqrt(y)) s=sum(x*x) y=1/rchisq(n=1, df=max(x+s,1)) #print(c(i,s,x,y)) mat[i, ]=c(x, y) } mat } bvn<-gibbs(100000,k=10,mean=0) plot(bvn[,c(1,k)]) plot(bvn[,1]) plot(bvn[,2],ylim=c(-100,100)) plot(bvn[,k+1],ylim=c(0,1000)) #cor(bvn)