R/generics.R

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

### ===========================================================================
### GENERICS
### This package supports the following generics:
### -- plot, points, lines
### -- summary, print
### ===========================================================================

#' @export
`[[.kdensity` = function(x, i) {
  allowed_arg = c("x", "bw_str", "bw", "adjust", "h", "kernel_str", "kernel",
                  "start", "start_str", "support", "data.name", "n", "range",
                  "has.na", "na.rm", "normalized", "call", "estimates", "logLik")

  i = match.arg(i, allowed_arg)
  environment(x)[[i]]
}

#' @export
`[[<-.kdensity` = function(x, i, value) {
  allowed_arg = c("x", "bw", "adjust", "kernel", "start", "support", "na.rm",
                  "normalized")

  i = match.arg(i, allowed_arg)
  environment(x)$obj_name = "x"
  args = list(object = x)
  args[[i]] = value
  do.call(update.kdensity, args)
  x
}


#' @export
`$<-.kdensity` = function(x, name, value) {
  x[[name]] = value
  x
}


#' @export
`$.kdensity` = function(x, name) {
  x[[name]]
}

#' @export
update.kdensity = function(object, ...) {


  current = list(x = object$x,
                 bw = object$bw_str,
                 adjust = object$adjust,
                 kernel = object$kernel_str,
                 start = object$start_str,
                 support = object$support,
                 na.rm = object$na.rm,
                 normalized = object$normalized)

  passed = list(...)

  ## Part of a hack to make $<- and [[<- work.
  if(!is.null(environment(object)$obj_name)) {
    obj_name = environment(object)$obj_name
  } else {
    obj_name = deparse(substitute(object))
  }

  arg_names = lapply(match.call(expand.dots=TRUE)[-1], deparse)
  args = listmerge(current, passed, type = "template")
  new_object = do.call(kdensity, args)

  if("x" %in% names(arg_names)) {
    data.name = arg_names$x
  } else {
    data.name = object$data.name
  }

  call = call("kdensity", x = data.name, adjust = args$adjust,
                kernel = args$kernel, start = args$start, support = args$support)

  environment(new_object)$call = call

  environment(new_object)$data.name = data.name
  assign(obj_name, new_object, envir = parent.frame())
}

#' @export
coef.kdensity = function(object, ...) object$estimates

#' @export
logLik.kdensity = function(object, ...) {
  msg = "'logLik' only makes sense for kdensity objects with a non-uniform parametric start."
  assertthat::assert_that(object$start_str != "uniform" & object$start_str != "constant", msg = msg)
  val = object$logLik
  attr(val, "nobs") = length(object$n)
  attr(val, "df")   = length(stats::coef(object))
  class(val) = "logLik"
  val
}

confint.kdensity = function(object, parm, level = 0.95, ...) {
  # Implement pointwise confidence intervals in some way.
}

#' Supplies a plotting range from a kdensity object.
#'
#' @param obj A kdensity object.
#' @keywords internal
#' @return S vector of size 1000, used for plotting.

get_range = function(obj) {
  support = obj$support
  minimum = obj$range[1]
  maximum = obj$range[2]
  obj_range = maximum - minimum
  addition = obj_range/10
  xmin = max(minimum - addition, support[1])
  xmax = min(maximum + addition, support[2])
  seq(xmin, xmax, length.out = 1000)
}

#' Plot, Lines and Points Methods for Kernel Density Estimation
#'
#' The `plot` method for `kdensity` objects.
#'
#' @export
#' @param x a `kdensity` object.
#' @param range range of x values.
#' @param plot_start logical; if `TRUE`, plots the parametric start instead of the kernel density estimate.
#' @param zero_line logical; if `TRUE`, add a base line at `y = 0`.
#' @param ... further plotting parameters.
#' @return None.
#' @seealso [kdensity()]
#' @examples
#' ## Using the data set "precip" to eye-ball the similarity between
#' ## a kernel fit, a parametric fit, and a kernel with parametric start fit.
#' kde_gamma = kdensity(precip, kernel = "gaussian", start = "gamma")
#' kde = kdensity(precip, kernel = "gaussian", start = "uniform")
#'
#' plot(kde_gamma, main = "Annual Precipitation in US Cities")
#' lines(kde_gamma, plot_start = TRUE, lty = 2)
#' lines(kde, lty = 3)
#' rug(precip)

plot.kdensity = function(x, range = NULL, plot_start = FALSE, zero_line = TRUE, ...) {
  plot_helper(x, range, plot_start, zero_line, ptype = "plot", ...)
}

#' @rdname plot.kdensity
#' @export
lines.kdensity = function(x, range = NULL, plot_start = FALSE, zero_line = TRUE, ...) {
  plot_helper(x, range, plot_start, zero_line, ptype = "lines", type = "l", ...)
}

#' @rdname plot.kdensity
#' @export
points.kdensity = function(x, range = NULL, plot_start = FALSE, zero_line = TRUE, ...) {
  plot_helper(x, range, plot_start, zero_line, ptype = "points", type = "p", ...)
}

#' Helper function for the plot methods.
#'
#' A helper function for the plot methods that does most of the work under
#' the hood.
#'
#' @param x A `kdensity` object.
#' @param range An optional range vector; like `x` in `plot.default`.
#' @param plot_start Logical; if `TRUE`, plots the parametric start only.
#' @param zero_line Logical; if `TRUE`, adds a line at `y = 0`.
#' @param ptype The kind of plot to make
#' @param ... Passed to plot.default.
#' @keywords internal
#' @return None.

plot_helper = function(x, range = NULL, plot_start = FALSE, zero_line = TRUE, ptype = c("plot", "lines", "points"), ...) {

  ptype = match.arg(ptype)

  if(is.null(range)) range = get_range(x)

  ## Potential arguments in ellipses are handled here. They are modelled
  ## after the structure of the 'density' function.
  supplied = list(...)

  bw_string = NULL
  if(!is.null(x$bw_str)) bw_string = paste0(" ('", x$bw_str, "')")

  defaults = list(type = "l",
                  main = deparse(x$call),
                  ylab = "Density",
                  xlab = paste0("N = ", x$n, "   Bandwidth = ", formatC(x$bw),
                               bw_string),
                  lwd  = 1)

  args = listmerge(x = defaults,
                   y = supplied)
  args$x = range

  if(plot_start) {
    start = x$start_str

    msg = "To use 'plot_start = TRUE', supply a parametric start that is a proper density."
    assertthat::assert_that(!is.null(start), start != "uniform", msg = msg)

    start = get_start(start)

    parameters = stats::coef(x)
    parametric_start = start$density

    args$y = sapply(range, function(y) {
      do.call(parametric_start, as.list(c("x" = y, parameters)))})

  } else {
    args$y = x(range)
  }

  switch(ptype,
         plot   = do.call(graphics::plot, args),
         lines  = do.call(graphics::lines, args),
         points = do.call(graphics::points, args))

  if(zero_line) graphics::abline(h = 0, lwd = 0.1, col = "gray")

  invisible(x)

}

#' @export
print.kdensity <- function(x, ...) {
  digits = list(...)$digits
  cat("\nCall:\n", deparse(x$call), "\n\n",
      "Data:      ", x$data.name, " (",x$n, " obs.)\n",
      "Bandwidth: ", formatC(x$bw, digits = digits), " ('", x$bw_str, "')\n",
      "Support:   (", x$support[1], ", ", x$support[2],   ")\n",
      "Kernel:    ", x$kernel_str, "\n",
      "Start:     ", x$start_str, "\n\n",
      sep = "")
  invisible(x)
}

#' @export
summary.kdensity <- function(object, ...)
{
  digits = list(...)$digits
  parameters = object$estimates
  params = NULL
  if(length(parameters) > 0)
  params = c("Parameter estimates:", "\n",
  sapply(1:length(parameters), function(i) paste0(names(parameters)[i], ": ", formatC(parameters[i], digits), "\n")),
  "\n")
  cat("\nCall: \n", deparse(object$call), "\n\n",
      "Data:        ", object$data.name, " (",object$n, " obs.)\n",
      "Bandwidth:   ", formatC(object$bw, digits = digits), " ('", object$bw_str, "')\n",
      "Support:     (", object$support[1], ", ", object$support[2],   ")\n",
      "Kernel:      ", object$kernel_str, "\n",
      "Start:       ", object$start_str, "\n",
      "Range:       (", formatC(object$range[1], digits), ", ", formatC(object$range[2], digits),   ")\n",
      "NAs in data: ", object$has.na, "\n",
      "Adjustment:  ", object$adjust, "\n\n",
      params,
      sep = "")
  invisible(object)
}