# Maximum Likelihood Estimation for the Battle of the Sexes Example
#
# ADM 6/5/2006

bos.ll <- function(theta, y, n) {
   pi <- 2 * theta / (1 + theta)^2
   output <- y*log(pi) + (n-y)*log(1-pi)
   return(output)
}

# plot the likelihood function for Case 1
ruler <- seq(1,5,0.01)
loglikelihood <- bos.ll(ruler, y=39, n=100)
pdf(file="bosll.pdf", width=11, height=8.5)
plot(ruler, loglikelihood, type="l", lwd=2, col="green", 
   xlab="theta", ylab="L(theta)",
   main="Log-Likelihood for Battle of the Sexes Example")
dev.off()

# optimize the log-likelihood
case1 <- optim(c(5),
   bos.ll,                      # the likelihood function   
   method="BFGS",               # optimization method
   hessian=TRUE,                # return numerical Hessian
   control=list(fnscale=-1),    # maximize function instead of minimize  
   y=39, n=100)                 # the data

case2 <- optim(c(5),
   bos.ll,                      # the likelihood function   
   method="BFGS",               # optimization method
   hessian=TRUE,                # return numerical Hessian
   control=list(fnscale=-1),    # maximize function instead of minimize  
   y=26, n=100)                 # the data

cat("Case1 Summary", "\n n = ", 100, "\n y = ", 39,
   "\n MLE of theta = ", case1$par, "\n Standard Error = ",
   sqrt(-1 * 1/case1$hessian), "\n")
cat("\n")
cat("Case2 Summary", "\n n = ", 100, "\n y = ", 26,
   "\n MLE of theta = ", case2$par, "\n Standard Error = ",
   sqrt(-1 * 1/case2$hessian), "\n")
cat("\n")