R/vieu.R

In modeest: Mode Estimation

#' @title
#' Vieu's mode estimator
#'
#' @description
#' Vieu's mode estimator is the value at which the kernel density derivative
#' estimate is null.
#'
#' @note
#' The user may call \code{vieu} through
#' \code{mlv(x, method = "vieu", ...)}.
#'
#' Presently, \code{vieu} is quite slow.
#'
#' @references
#' \itemize{
#'   \item Vieu P. (1996). A note on density mode estimation.
#'   \emph{Statistics & Probability Letters}, \bold{26}:297--307.
#' }
#'
#' @param x
#' numeric. Vector of observations.
#'
#' @param bw
#' numeric. The smoothing bandwidth to be used.
#'
#' @param kernel
#' character. The kernel to be used. Available kernels are \code{"biweight"},
#' \code{"cosine"}, \code{"eddy"}, \code{"epanechnikov"}, \code{"gaussian"},
#' \code{"optcosine"}, \code{"rectangular"}, \code{"triangular"},
#' \code{"uniform"}. See \code{\link[stats]{density}} for more details on some
#' of these kernels.
#'
#' @param abc
#' logical. If \code{FALSE} (the default), the root of the density derivate
#' estimate is searched with \code{\link[stats]{uniroot}}.
#'
#' @param ...
#' If \code{abc = FALSE}, further arguments to be passed to
#' \code{\link[stats]{uniroot}}.
#'
#' @return
#' \code{vieu} returns a numeric value, the mode estimate. If \code{abc = TRUE},
#' the \code{x} value at which the density derivative estimate is null is
#' returned. Otherwise, the \code{\link[stats]{uniroot}} method is used.
#'
#' @importFrom stats uniroot
#' @importFrom statip .kernelsList bandwidth kernelfun
#' @export
#' @aliases Vieu
#'
#' @seealso
#' \code{\link[modeest]{mlv}}, \code{\link[modeest]{parzen}}.
#'
#' @examples
#' # Unimodal distribution
#' x <- rlnorm(10000, meanlog = 3.4, sdlog = 0.2)
#'
#' ## True mode
#' lnormMode(meanlog = 3.4, sdlog = 0.2)
#'
#' ## Estimate of the mode
#' mlv(x, method = "vieu", kernel = "gaussian")
#'
vieu <- function(x, bw = NULL, kernel = "gaussian", abc = FALSE, ...) {
  if (pmatch(tolower(kernel), "normal", nomatch = 0L)) {
    kernel <- "gaussian"
  } #else {
    #kernel <- match.arg(tolower(kernel), c(statip::.kernelsList(), "uniform"))
  #}

  if (is.null(bw)) bw <- "nrd0"
  f <- statip::densityfun(x, bw = bw, kernel = kernel, ...)
  bw <- statip::bandwidth(x, bw)

  fn <- function(z) {
    k <- statip::kernelfun(kernel, derivative = TRUE)((z - x) / bw)
    sum(k)
  }

  if (!abc) {
    r <- stats::uniroot(f = fn, interval = range(x), ...)
    M <- r$root
  } else {
    .NotYetImplemented()
    #FN <- Vectorize(fn)
    #f <- abs(FN(x))
    #M <- x[f == min(f)]
    # TODO:
  }
  M
}