### R DEMO
library(MASS)
library(ggplot2)

# Negative log likelihood for a Weibull sample
nll_weibull <- function(pars, y) {
  # Handle the case of negative parameter values
  if (isTRUE(any(pars <= 0))) {
    # parameters must be positive
    return(1e10) # large value (not infinite, to avoid warning messages)
  }
  - sum(dweibull(
    x = y,
    scale = pars[1],
    shape = pars[2],
    log = TRUE
  ))
}
# Gradient of the Weibull negative log likelihood
gr_nll_weibull <- function(pars, y) {
  scale <- pars[1]
  shape <- pars[2]
  n <- length(y)
  grad_ll <- c(
    scale = -n * shape / scale + shape * scale^(-shape - 1) * sum(y^shape),
    shape = n / shape - n * log(scale) + sum(log(y)) -
      sum(log(y / scale) * (y / scale)^shape)
  )
  return(-grad_ll)
}

# load data from the package MASS
data(geyser, package="MASS")
waiting <- geyser$waiting #waiting time between eruptions of the geyser
# Use exponential submodel MLE as starting parameters
start <- c(mean(waiting), 1) #with alpha=1, E(Y)= lambda
# Check gradient function is correctly coded!
# Returns TRUE if numerically equal to tolerance
isTRUE(all.equal(
  numDeriv::grad(nll_weibull, x = start, y = waiting),
  gr_nll_weibull(pars = start, y = waiting),
  check.attributes = FALSE
))
# Numerical minimization using optim
opt_weibull <- optim(
  par = start,
  # starting values
  fn = nll_weibull,
  # pass function, whose first argument is the parameter vector
  gr = gr_nll_weibull,
  # optional (if missing, numerical derivative)
  method = "BFGS",
  # gradient-based algorithm, common alternative is "Nelder"
  y = waiting,
  # vector of observations, passed as additional argument to fn
  hessian = TRUE
) # return matrix of second derivatives evaluated at MLE

# Parameter estimates - MLE
mle_weibull <- opt_weibull$par
# Check gradient for convergence
gr_nll_weibull(mle_weibull, y = waiting)
# Is the Hessian matrix positive definite (all eigenvalues are positive)
# If so, we found a maximum and the matrix is invertible
isTRUE(all(eigen(opt_weibull$hessian)$values > 0))

# The Hessian matrix of the negative log likelihood
# evaluated at the MLE (observed information matrix)
obsinfo_weibull <- opt_weibull$hessian
vmat_weibull    <- solve(obsinfo_weibull)
# Standard errors
se_weibull <- sqrt(diag(vmat_weibull))

nll <- matrix(nrow = 101, ncol = 100)
alpha <- seq(mle_weibull[2] - 2.5 * se_weibull[2],
             mle_weibull[2] + 2.5 * se_weibull[2],
             length.out = 100)
lambda <- seq(mle_weibull[1] - 2.5 * se_weibull[1],
              mle_weibull[1] + 2.5 * se_weibull[1],
              length.out = 101)
z <- rep(NA, length(nll))
for (i in seq_along(lambda)) {
  for (j in seq_along(alpha)) {
    z[(i - 1) * 100 + j] <- nll[i, j] <-
      nll_weibull(pars = c(lambda[i], alpha[j]), y = waiting)
  }
}


# Plot likelihood surface for the Weibull model with 10%, 20%, ..., 90% likelihood ratio confidence regions (contour curves). 
# Higher log likelihood values are indicated by darker colours.
ggplot() +
  geom_raster(data = data.frame(
    x = rep(lambda, each = length(alpha)),
    y = rep(alpha, length.out = length(alpha) *
              length(lambda)),
    z = c(-z + opt_weibull$value)
  ),
  mapping = aes(
    x = x,
    y = y,
    fill = pchisq(-2 * z, df = 2)
  )) +
  geom_contour(
    data = data.frame(
      x = rep(lambda, each = length(alpha)),
      y = rep(alpha, length.out = length(alpha) *
                length(lambda)),
      z = c(-z + opt_weibull$value)
    ),
    mapping = aes(x = x, y = y, z = z),
    col = "white",
    breaks = -qchisq(seq(0.1, 0.9, by = 0.1), df = 2) / 2
  ) +
  # Add cross for MLE
  geom_point(
    shape = 4,
    color = "white",
    data = data.frame(x = mle_weibull[1], y = mle_weibull[2]),
    mapping = aes(x = x, y = y)
  ) +
  # Use a different color palette
  scale_fill_viridis_c(direction = 1, option = "viridis") +
  labs(x = expression(paste("scale ", lambda)),
       y = expression(paste("shape ", alpha)),
       fill = "probability level") +
  scale_y_continuous(expand = expansion(), limits = range(alpha)) +
  scale_x_continuous(expand = expansion(), limits = range(lambda)) +
  theme(legend.position = "none")