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(rexp(n*m,10),m) # creates m by n matrix of exp(10) values
 for(j in 1:m) theta.hat[j] <- 1/mean(x[j,]) # creates list of m MLE's 
 hist(theta.hat,main=paste("n=",n))
}