n <- 16 a <- 2 b <- 4 S <- matrix(0,2,10) S[2,1] <- 0.5 for(i in 2:ncol(S)) { S[1,i] <- rbinom(1,n,S[2,i-1]) S[2,i] <- rbeta(1,a+S[1,i],b+n-S[1,i]) } gibbs <- function(n, a, b, iter = 1000) { S <- matrix(0,iter+1,2) Y <- 0.5 S[1,2] <- Y for (i in 1:iter) { X <- rbinom(1,n,Y) Y <- rbeta(1,a+X,b+n-X) S[i+1,] <- c(X,Y) } return(S) } S <- gibbs(16,2,4) X <- S[-(1:101),1] mean(X==2) choose(n,2)*beta(2+a,n-2+b)/beta(a,b) X <- S[-(1:501),1] mean(X==2) S <- gibbs(16,2,4,10^5) X <- S[-(1:1001),1] mean(X==2)