R/eel.mean.R

In markean/bayesEL: Bayesian Empicial Likelihood

#' Extented empirical likelihood for mean
#'
#' Extented empirical likelihood for mean
#'
#' @param theta a vector of parameters to be tested.
#' @param x a matrix or vector of data. Each row is an observation vector.
#' @param abstol an optional value for the absolute convergence tolerance. Defaults to 1e-8.
#' @param maxit an optional value for the maximum number of iterations. Defaults to 50.
#'
#' @return An object with S3 class "eel"
#'
#' @examples
#' theta <- 10
#' x <- rnorm(100)
#' eel.mean(theta, x)
#'
#' @import Matrix
#' @importFrom stats uniroot dist
#' @export
eel.mean <- function(theta, x, abstol = 1e-8, maxit = 50) {
  ## Data
  x <- as.matrix(x)
  if (length(theta) != ncol(x))
    stop("length(theta) != ncol(x)")
  if (rankMatrix(x) != ncol(x))
    stop("model matrix must have full column rank")
  colnames(x) <- NULL
  n <- nrow(x)
  # distance between xbar and theta
  distance <- sapply(dist(rbind(theta, colMeans(x)), method = "euclidean"),
                     function(x) {attributes(x) <- NULL; x})
  # nothing to extend when theta == xbar
  if (distance == 0) {
    el_result <- el.mean(colMeans(x), x, abstol, maxit)
    result <- list()
    class(result) <- "eel"
    result$logLR <- el_result$logLR
    result$logL <- el_result$logL
    result$w <- el_result$w
    result$lambda <- el_result$lambda
    result$iterations <- el_result$iterations
    # result$message <- el_result$message
    return(result)
  }
  direction_unit_vec <- (theta - colMeans(x)) / distance


  ## EEL
  # simplified expansion factor
  f <- function(alpha) alpha * (1 - el2.mean(colMeans(x) + alpha * direction_unit_vec, x)$logLR / n) - distance
  # wrapper for f to remove Inf value
  g <- function(alpha) ifelse(f(alpha) == Inf, .Machine$double.xmax, f(alpha))
  # compute the expansion distance
  alpha <- uniroot(g, lower = 0, upper = distance, extendInt = "upX",
                     tol = .Machine$double.eps^.25, check.conv = F)$root
  # preimage of the expansion
  theta_prime <- colMeans(x) + alpha * direction_unit_vec


  ## Result
  el.result <- el.mean(theta_prime, x, abstol, maxit)
  result <- list()
  class(result) <- "eel"
  result$logLR <- el.result$logLR
  result$logL <- el.result$logL
  result$w <- el.result$w
  result$lambda <- el.result$lambda
  result$iterations <- el.result$iterations
  result$message <- ifelse(all.equal(sum(result$w), 1) == T,
                           "preimage is in the convex hull",
                           "something wrong"
  )
  result
}