numIter <- 1000;
n <- 10;
theta <- 3;
mle <- rep(0,numIter);

for (i in 1:numIter) {
  x <- rexp(n,theta);
  mle[i] <- 1/mean(x);
}

par(mfrow=c(2,1))
hist(mle,nclass=30,xlim=c(mean(mle)-3*sqrt(var(mle)),mean(mle)+3*sqrt(var(mle))));
par(new=TRUE);
mle<-sort(mle);
plot(mle,dnorm(mle,mean(mle),sqrt(var(mle))),xlim=c(mean(mle)-3*sqrt(var(mle)),mean(mle)+3*sqrt(var(mle))),type="l");

qqnorm(mle);
cat("Mean of the MLE: ",mean(mle),"\n");
cat("Variance of the MLE: ",var(mle),"\n");
cat("Fisher Variance of the MLE: ",theta^2/n,"\n");