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

###################################################################################
## Multivariate Normal example with p=2, n=100
library(mvtnorm)

Mu <- matrix(c(3,4),2,1)
Sigma <- matrix(c(2,1,1,3), 2,2)

y <- as.matrix(rmvnorm(100, Mu, Sigma))
plot(y[,1], y[,2], pch=16)
n <- nrow(y)

Mu.hat <- apply(y,2,mean)
Sigma.hat <- ((n-1)/n)*cov(y) 


## or utilizing R functions for the poisson density
Log.Likelihood <- function(theta){
Mu <- matrix(c(theta[1],theta[2]),2,1)
Sigma <-  matrix(c(theta[3],theta[4],theta[4],theta[5]), 2,2)

out <- sum(dmvnorm(y, Mu, Sigma,log=T))
return(out)
}


## test the function
theta.test <- c(Mu.hat, Sigma.hat)
theta.test

## remove the one of the covariances for the optimization
theta.test <- theta.test[-4]
theta.test
Log.Likelihood(theta.test)


## use optim()
optim(par=c(0,0,1,0,1), fn=Log.Likelihood,  control = list(fnscale=-1), method="BFGS")