####### Maximum likelihood - from handout #######
# likelihood for poisson
pois.lhood <- function(theta, y){
  prod(theta ^ y * exp(-theta) / factorial(y))
}

# some fictional data
y <- c(5, 0, 6, 3, 7, 5, 5, 7, 3, 5)

# want to plot likelihood for theta
theta <- seq(0, 10, 0.1)
l.at.theta <- sapply(theta, pois.lhood, y = y)

par(mfrow = c(2,1))
plot(theta, l.at.theta, type = "l", ylab = "Likelihood")

# log-likelihood for poisson
pois.llhood <- function(theta, y){
  sum(log(theta) * y - theta - log(factorial(y)))
}

ll.at.theta <- sapply(theta, pois.llhood, y = y)
plot(theta, ll.at.theta, type = "l", ylab = "Log likelihood")

# maximise by numerical methods 
# use function nlm - but this is a minimiser so we minimise negative log-likelihood
neg.pois.llhood <- function(theta, y){
  -pois.llhood(theta, y)
}

# nlm - The first argument should be the name of the function that needs 
# minimizing. This function (in our case neg.pois.llhood) should have as 
# its first argument the parameter we wish to minimize over (in our case theta).
# The second argument is a starting value for the parameter of interest.
# Any remaining arguments are passed on to the function we are minimizing.

nlm(neg.pois.llhood, p = 5, y = y)
# MLE estiamte of theta is 4.6

abline(v = nlm(neg.pois.llhood, p = 5, y = y)$estimate)

# We can actually calculate the maximum likelihood estimate of theta
# analytically. Do it for practice and check we get the same answer

####### 2d function minimization ######
# If we want to minimize a 2d function, we need to write the function 
# to take a vector argument.

# For example, find the minimum of the quadratic x^2 + y^2 
quad.func <- function(z){
    x <- z[1]
    y <- z[2]
    x^2 + y^2
}

# ok and minimizing it
nlm(quad.func, p = c(0.1, 0.2))

x<- seq(0,10,0.1)
y<- seq(0,10,0.1)
z <- expand.grid(x = x, y = y)
height <- matrix(apply(z, 1, quad.func), ncol = length(y))

contour(x, y, height)
image(x, y, height)
persp(x, y, height)