# Maximum Likelihood Estimation for the Exponential Distribution
# 
# ADM 6/5/2006

# simulate some fake data
lambda <- 2
n <- 1000
y <- rexp(n, 1/lambda)

# log-likelihood function
log.likelihood <- function(lambda, data) {
  contributions <-  dexp(data, 1/lambda, log=TRUE)
  return(sum(contributions))
}

# plot log-likelihood
ruler <- seq(1,3,0.01)
ll <- sapply(ruler, log.likelihood, data=y)
pdf(file="exponential.pdf", width=11, height=8.5)
plot(ruler, ll, type="l", lwd=2, col="red", 
   xlab="lambda", ylab="L(lambda)",
   main="Log-Likelihood for Exponential Distribution")
dev.off()

# numerically optimize using optim()
mle <- optim(c(1),
   log.likelihood,              # the likelihood function   
   method="BFGS",               # optimization method
   hessian=TRUE,                # return numerical Hessian
   control=list(fnscale=-1),    # maximize function instead of minimize  
   data=y)                      # the data
cat("@@@@@ MLE Summary for Exponential Example @@@@@ \n",
    "MLE of lambda (numerical) = ", mle$par,
    "\n Standard Error (numerical) = ", sqrt(-1 * 1/mle$hessian), "\n")
cat("\n MLE of lambda (analytic) = ", mean(y),
    "\n Standard Error (analytic) = ", sqrt(mean(y)^2 / n), "\n")
cat("\n")

# Monte Carolo analysis of coverage properties
G <- 5000
cover <- matrix(NA, G, 1)
for(g in 1:G) {
   y <- rexp(n, 1/lambda)
   mle <- mean(y)
   se <- sqrt(mean(y)^2 / n)
   if(lambda > mle-1.96*se & lambda < mle+1.96*se) cover[g] <- 1
   else cover[g] <- 0
}
cat(" Covered", sum(cover), "out of", G, "trials.\n")
cat(" Coverage Rate:", sum(cover)/G, "\n\n")