## October 4, Maximum Likelihood 
rm(list = ls())

###################################################################################
## Poisson Example
salmon.count <- rpois(31, 10)  ## generate fake data for a month
y <- salmon.count

## we know the mle is the mean
mle.lambda.1 <- mean(y)
mle.lambda.1

## let's use optim to get the mle
log.likelihood.1 <- function(lambda){
LL <- 0
for(i in 1:length(y)){
LL <- LL + log(( exp(-lambda)*lambda^(y[i]) ) / (factorial(y[i]) ))
}
return(LL)
}

## or utilizing R functions for the poisson density
log.likelihood.2 <- function(lambda){
out <- sum(dpois(y, lambda, log=T))
return(out)
}

## test the functions
log.likelihood.1(lambda=2)
log.likelihood.2(lambda=2)

## now optimize wrt to lambda
optim(par=1, fn=log.likelihood.2,  control = list(fnscale=-1), method="BFGS")
optim(par=10, fn=log.likelihood.2,  control = list(fnscale=-1), method="BFGS")


#################################################################################
## Gamma example

light.bulbs <- rgamma(100, shape=2, scale=4)
hist(light.bulbs, xlab="life time of bulbs in months", cex.lab=1.5, col="yellow")

##
y <- light.bulbs

## let's take the log likelihood equation for alpha (shape) after substituting beta (scale) 

alpha.func <- function(alpha){
n <- length(y)
out <- -n*lgamma(alpha) - alpha*n*log(mean(y)*(1/alpha)) + (alpha-1) *sum(log(y)) - alpha*n
}

## now optimize wrt to alpha
optim(par=1, fn=alpha.func,  control = list(fnscale=-1), method="BFGS")
alpha.mle <-optim(par=1, fn=alpha.func,  control = list(fnscale=-1), method="BFGS")$par
beta.mle <- mean(y)/alpha.mle

alpha.mle
beta.mle


## now let's just use the dgamma() command in R
Log.Likelihood <- function(theta){
alpha <- theta[1]
beta <- theta[2]
out <- sum(dgamma(y, shape=alpha, scale=beta, log=T))
}

## now optimize wrt to alpha and beta
optim(par=c(1,1), fn=Log.Likelihood,  control = list(fnscale=-1), method="BFGS")
optim(par=c(10,10), fn=Log.Likelihood,  control = list(fnscale=-1), method="BFGS")