par(mfrow=c(2,2)) # 2 by 2 array of plots
m <- 10000 # number of MLE's generated for each plot
theta.hat <- rep(0,m) # vector of MLE'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] <- -1/mean(log(x[j,])) # creates list of m MLE's
 hist(theta.hat,main=paste("MLE: n=",n))
}