par(mfrow=c(2,2)) # 2 by 2 array of plots
m <- 10000 # number of MOM's generated for each plot
theta.hat <- rep(0,m) # vector of MOM's for each plot
for(i in 1:4){
 n <- 5^i # sample sizes are n=5, 25, 125, 625
 x <- matrix(rbeta(n*m,5,1),m) # creates m by n matrix of beta(5,1) values
 for(j in 1:m) theta.hat[j] <- mean(x[j,])/(1-mean(x[j,])) # creates list of m MOM's
 hist(theta.hat,main=paste("MOM: n=",n))
}