Recursive function to quickly search for the parameter that maximizes a likelihood function

maximum_likelihood.R

# ============================================================================ #
# Quickly find the parameter that maximizes a likelihood function.
#
# This recursive function will take a range of n parameters within an interval
# and re-adjust the interval to border the maximum value and recursively call
# the function. Once the difference between the maximum value found and the 
# previous maximum value is less than the error term, the parameter value is
# returned. This is faster than the brute force method as the range is reduced
# every time the function is recursed.
#
# PARAMETERS:
#   LFUN     vectorized likelihood function
#   interval numeric inputs specifying the search space (inclusive)
#   n        the number of parameters to generate between the interval
#   max      set to zero. This will represent the current maximum likelihood 
#            value
#   e        an error term
#   quiet    if set to FALSE, a message with the parameter and likelihood values
#            will be printed to screen
# OUTPUT:
#   A numeric value representing the parameter that maximizes the likelihood
#   function. 
# ============================================================================ #
ml_parameter <- function(LFUN, interval = c(0, 1), n = 101, max = 0,
                         e = .Machine$double.eps^0.75, quiet = TRUE){
  interval <- seq(from = interval[1], to = interval[2], length = n)
  lik <- LFUN(interval)
  max.val <- which.max(lik)
  max.lik <- lik[max.val]
  max.int <- interval[max.val]
  if (abs(log1p(max) - log1p(max.lik)) < log1p(e) && log1p(max) <= log1p(max.lik)){
    if (!quiet){
      cat("PARAMETER:\t", max.int, "\tLIKELIHOOD:\t", max.lik, "\n")
    }
    return(max.int)
  } else if (log1p(max) > log1p(max.lik)) {
    return(NULL)
  } else if (max.val == 1){
    return(ml_parameter(LFUN, c(interval[1], interval[2]), n, max.lik, e, quiet))
  } else {
    return(ml_parameter(LFUN, interval[max.val + c(-1, 1)], n, max.lik, e, quiet))
  }
}

# ============================================================================ #
# # Uncomment this to run an example from the R package adegenet:
# 
# library(adegenet)
# data(microbov)
# sal <- seppop(microbov)$Salers
# Fsal <- inbreeding(sal, res.type = "function")
# Fest <- sapply(Fsal, ml_parameter)
# plot(Fest)
# ============================================================================ #

# ============================================================================ #
# Benchmarking
#
# This benchmark is based on the function below, which will find the parameter
# that maximizes the likelihood function among n parameters.
# ============================================================================ #
seq_max <- function(LFUN, n = 1001, interval = c(0, 1)){
  ran <- seq(from = interval[1], to = interval[2], length= n)
  ran[which.max(LFUN(ran))]
}

# ============================================================================ #
# # Uncomment this section to perform the benchmark
#
# F37 <- Fsal[[37]] ## Maximum inbreeding coefficient
# (mlp.res <- ml_parameter(F37))
# (sm.res <- seq_max(F37)) # Note, this is with 1001 parameters
# (p1k <- F37(c(mlp.res, sm.res)))
# which.max(p1k)
# diff(p1k)
# sm.res <- seq_max(F37, n = 10001) # increasing the number of iterations
# (p10k <- F37(c(mlp.res, sm.res)))
# which.max(p10k)
# diff(p10k)
# library(microbenchmark)
# library(ggplot2)
# (mb.res <- microbenchmark(ml_parameter(F37), seq_max(F37)))
# autoplot(mb.res)
# ============================================================================ #

Note: this function bows before the stats function, optimize, which is faster and more rigorous.

library(adegenet)
data(microbov)
sal  <- seppop(microbov)$Salers
# Generate likelihood functions
Fsal <- inbreeding(sal, res.type = "function")

## Calculating maximum likelihood values
# Method above
Fest <- vapply(Fsal, ml_parameter, numeric(1))
# stats::optimize
op   <- vapply(lapply(Fsal, optimize, interval = c(0, 1), maximum = TRUE, 
                    tol = .Machine$double.eps^0.75), unlist, numeric(2))
# Getting likelihood values for method above
Festml <- vapply(1:length(Fsal), function(x) Fsal[[x]](Fest[x]), numeric(1))

## PLOTTING
mlcol <- ifelse(op[2, ] - Festml > 0, "red", 
                ifelse(op[2, ] - Festml < 0, "blue", "black"))
yl <- "Likelihood difference"
plot(abs(op[2, ] - Festml), col = mlcol, log = "y", ylab = yl)
abline(h = .Machine$double.eps, lty = 2)
text(x = 0, y = .Machine$double.eps, labels = "epsilon", adj = c(0, -0.5))
legend("bottomleft", legend = c("optimize", "ml_parameter"),bty = "n", 
       pch = c(1, 1), col = c("red", "blue"), cex = 0.75,
       title = "Greater likelihood")