R/kdensity.R

In JonasMoss/kdensity: Kernel Density Estimation with Parametric Starts and Asymmetric Kernels

#' Parametrically guided kernel density estimation
#'
#' `kdensity` computes a parametrically guided kernel density estimate
#' for univariate data. It supports asymmetric kernels and parametric starts
#' through the `kernel` and `start` arguments.
#'
#' The default values for `bw`, `kernel`, `start`, and
#' `support` are interdependent, and are chosen to make sense. E.g.,
#' the default value for `support` when `start = beta` is
#' `c(0, 1)`.
#'
#' The `start` argument defaults to `uniform`, which corresponds
#' to ordinary kernel density estimation. The typical default value for
#' `kernel` is `gaussian`.
#'
#' @export
#'
#' @param x Numeric vector containing the data.
#'
#' @param bw A bandwidth function. Can be either a string, a custom-made
#' function, or a double. The supported bandwidth functions are documented
#' in [bandwidths()].
#'
#' @param adjust An adjustment constant, so that `h = adjust*bw*sd`, where `sd`
#' varies with the chosen kernel.
#'
#' @param kernel The kernel function. Can be chosen from the list of built-in
#' kernels or be custom-made. See [kernels()] for details.
#'
#' @param start Parametric start. Can be chosen from the list of built-in
#' parametric starts or be custom-made. See [parametric_starts()] for
#' details.
#'
#' @param support The support of the data. Must be compatible with the supplied
#' `x` and the supplied `start` and `kernel`. Is used to find the
#' normalization constant, see `normalized`.
#'
#' @param na.rm Logical; if `TRUE`, `NA`s will be removed from `x`.
#'
#' @param normalized Logical; if `TRUE`, the density is normalized.
#'
#' @param tolerance Numeric; the relative error to tolerate in normalization.
#'
#' @return `kdensity` returns an S3 function object of
#' [base::class()] "kdensity". This is a callable function with the
#' following elements, accessible by '$':
#' \describe{
#'   \item{`x`}{The data supplied in `x`.}
#'   \item{`bw_str, bw, adjust, h`}{The bandwidth function, the resulting
#'               bandwidth, the `adjust` argument, and the adjusted
#'               bandwidth.}
#'   \item{`kernel_str, kernel, start, start_str, support`}{Name of the kernel,
#'               the kernel object, name of the parametric start, the start object,
#'               and the support of the density.}
#'   \item{`data.name, n, range, has.na, na.rm, normalized`}{Name of the data, number of
#'               observations, the range of the data, whether the data
#'               `x` contained `NA` values, whether na.rm is `TRUE`
#'               or not, and whether the density is normalized.}
#'   \item{`call`}{The `call` to `kdensity`.}
#'   \item{`estimates`}{Named numeric vector containing the parameter
#'               estimates from the parametric start.}
#'   \item{`logLik`}{The log-likelihood of the parametric starts. Is `NA`
#'               for the uniform start.}
#'
#' }
#'
#'
#' @details If `normalized` is `FALSE` and `start != "uniform"`, the resulting
#' density will not integrate to 1 in general.
#'
#' @seealso The `stats` package function [stats::density()].
#' @references
#'   Hjort, Nils Lid, and Ingrid K. Glad. "Nonparametric density estimation with a parametric start." The Annals of Statistics (1995): 882-904.
#'
#'   Jones, M. C., and D. A. Henderson. "Miscellanea kernel-type density estimation on the unit interval." Biometrika 94.4 (2007): 977-984.
#'
#'   Chen, Song Xi. "Probability density function estimation using gamma kernels." Annals of the Institute of Statistical Mathematics 52.3 (2000): 471-480.
#'
#'   Silverman, Bernard W. Density estimation for statistics and data analysis. Vol. 26. CRC press, 1986.
#'
#' @examples
#' ## Use gamma kernels to model positive data, the concentration of
#' ## theophylline
#'
#' concentration <- Theoph$conc + 0.001
#' plot(kdensity(concentration, start = "gamma", kernel = "gamma", adjust = 1 / 3),
#'   ylim = c(0, 0.15), lwd = 2, main = "Concentration of theophylline"
#' )
#' lines(kdensity(concentration, start = "gamma", kernel = "gaussian"),
#'   lty = 2, col = "grey", lwd = 2
#' )
#' lines(kdensity(concentration, start = "gaussian", kernel = "gaussian"),
#'   lty = 3, col = "blue", lwd = 2
#' )
#' lines(kdensity(concentration, start = "gaussian", kernel = "gamma", adjust = 1 / 3),
#'   lty = 4, col = "red", lwd = 2
#' )
#' rug(concentration)
#'
#' ## Using a density and and estimator from another package.
#'
#' skew_hyperbolic <- list(
#'   density   = SkewHyperbolic::dskewhyp,
#'   estimator = function(x) SkewHyperbolic::skewhypFit(x, printOut = FALSE)$param,
#'   support   = c(-Inf, Inf)
#' )
#'
#' kde <- kdensity(diff(LakeHuron), start = skew_hyperbolic)
#' plot(kde,
#'   lwd = 2, col = "blue",
#'   main = "Annual differences in water level (ft) of Lake Huron, 1875 - 1972"
#' )
#' lines(kde, plot_start = TRUE, lty = 2, lwd = 2) # Plots the skew hyperbolic density.
#' rug(diff(LakeHuron))
#'
#' kde$estimates # Also: coef(kde)
#' # Displays the parameter estimates:
#' #        mu     delta      beta        nu
#' # -1.140713  3.301112  2.551657 26.462469
#'
kdensity <- function(x, bw = NULL, adjust = 1, kernel = NULL, start = NULL,
                     support = NULL, na.rm = FALSE, normalized = TRUE,
                     tolerance = 0.01) {
  ## These tests are based purely on the data, and should be run no matter what.
  data.name <- deparse(substitute(x))
  has.na <- anyNA(x)

  assertthat::assert_that(!(has.na & !na.rm),
    msg = "x contains NAs and na.rm = FALSE."
  )

  x <- x[!is.na(x)] # This line is reached only if (has.na & !na.rm) is FALSE.

  ## 'kernel', 'start' and 'bw' can be custom made: In this case, they must
  ## be added to their environments.

  if (!is.null(start)) {
    if (!is.character(start)) {
      start_str <- deparse(substitute(start))
      add_start(start_str = start_str, start = start)
      start <- start_str
    }
  }

  if (!is.null(kernel)) {
    if (!is.character(kernel)) {
      kernel_str <- deparse(substitute(kernel))
      add_kernel(kernel_str = kernel_str, kernel = kernel)
      kernel <- kernel_str
    }
  }

  ## The case of bw == Inf is special! In this case, we return the parametric
  ## start itself.

  ## Now we handle the combinations of kernel, start and support.
  ## This is fancy defaults management.
  kss_list <- get_kernel_start_support(kernel, start, support)

  ## The strings are used for reporting and inside bandwidth functions,
  ## at must be preserved.
  start_str <- kss_list$start_str
  kernel_str <- kss_list$kernel_str

  ## We overwrite the kernel, start, and support with what we obtained from kss.
  kernel <- kss_list$kernel
  start <- kss_list$start
  support <- kss_list$support

  ## Tests for incompabibilities in the supplied values.
  support_compatible(kernel, start, support)

  ## We continue by checking if the supplied values make sense.
  if (!all(x <= support[2] & x >= support[1])) {
    stop(
      "The supplied data x is not contained in the support: (",
      support[1], ", ", support[2], ")."
    )
  }

  estimates <- start$estimator(x)
  parametric_start <- start$density

  # Name of the variable where the density is evaluated. Typically x.
  x_name <- names(formals(start$density))[1]

  parametric_start_vector <- function(x) {
    arguments <- list()
    arguments[[1]] <- x
    names(arguments)[1] <- x_name
    arguments <- append(arguments, as.list(estimates))
    do.call(parametric_start, arguments)
  }

  ## Takes care of the bandwidth. Can be either a double, a string, or a
  ## function taking the arguments data, kernel, start support.

  if (is.null(bw)) {
    bw <- get_standard_bw(kernel_str, start_str, support)
  }

  if (!is.numeric(bw)) {
    if (is.character(bw)) {
      bw_str <- bw
      bw <- get_bw(bw)(x, kernel_str, start_str, support)
    } else {
      bw_str <- deparse(substitute(bw))
      bw <- bw(x, kernel_str, start_str, support)
    }
  } else {
    bw_str <- "user supplied"
    if (bw == Inf) {
      msg <- "bw = Inf does not work with a uniform start."
      assertthat::assert_that(start_str != "uniform", start_str != "constant", msg = msg)
    }
  }

  ## The parameter h is computed. The basic bandwidth is h = bw*adjust for the
  ## normal kernel, and is adjusted for all the other kernels so that the sd
  ## of the kernel equals h.
  sd <- ifelse(!is.null(kss_list$kernel$sd), kss_list$kernel$sd, 1)
  h <- bw * adjust * sd

  ## The denominator can be computed.
  parametric_start_data <- parametric_start_vector(x)

  ## Now we handle normalization.

  normalization <- 1

  if (normalized & bw != Inf) {
    integrand <- function(y) {
      sapply(y, function(y) {
        mean(1 / h * kernel$kernel(y, x, h) /
          parametric_start_data) *
          parametric_start_vector(y)
      })
    }

    integral <- tryCatch(
      stats::integrate(integrand,
        lower = support[1],
        upper = support[2]
      ),
      error = function(e) {
        stop(paste0(
          "Normalization error: The function will not integrate.",
          "Two common causes are: 1.) The kernel is non-smooth, ",
          "try a smooth kernel if possible. 2.) The supplied ",
          "support is incorrect."
        ))
      }
    )

    msg <- "The normalization constant has too large relative error.
   1.) try to normalize the data, or
   2.) change kernel, or
   3.) increase the 'tolerance' parameter for kdensity"

    assertthat::assert_that(!is.nan(integral$abs.error / integral$value),
      msg = msg
    )
    assertthat::assert_that(integral$abs.error / integral$value < tolerance,
      msg = msg
    )

    normalization <- integral$value
  }


  ## This is the main part of the returned object.

  return_function <- function(y) {
    if (missing(y)) {
      stop("A kdensity object is a density function. Call it with numerical input.")
    }

    if (h != Inf) {
      if (length(y) == 1) {
        mean(1 / h * kernel$kernel(y, x, h) * parametric_start_vector(y) /
          parametric_start_data) / normalization
      } else {
        sapply(y, function(y) {
          parametric_start_vector_y <- parametric_start_vector(y)
          1 / h * mean(kernel$kernel(y, x, h) * parametric_start_vector_y /
            parametric_start_data) / normalization
        })
      }
    } else {
      ## If bw = Inf, the parametric start is returned.
      arguments <- list()
      arguments[[1]] <- y
      names(arguments)[1] <- x_name
      arguments <- append(arguments, as.list(estimates))
      do.call(parametric_start, arguments)
    }
  }

  ## 'R6'-ish attributes:
  call <- match.call()
  range <- c(min(x), max(x))
  n <- length(x)
  logLik <- ifelse(start_str == "uniform" | start_str == "constant", NA,
    sum(parametric_start_vector(x))
  )

  ## S3 attributes:
  class(return_function) <- "kdensity"
  attr(return_function, "srcref") <- NULL
  return_function
}