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
### ===========================================================================

## Fitted objects are read-only: there are no `$<-`, `[[<-` or `update()`
## methods. `exact` is accepted so that getCall() / getElement() work, which
## lets stats::update.default() refit by re-evaluating the stored call.
## Names are matched partially, as `$` does for lists.

#' @export
`[[.kdensity` <- function(x, i, exact = TRUE) {
  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, name) {
  x[[name]]
}

#' @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."
  assert_(object$start_str != "uniform" & object$start_str != "constant", msg = msg)
  val <- object$logLik
  attr(val, "nobs") <- object$n
  attr(val, "df") <- length(stats::coef(object))
  class(val) <- "logLik"
  val
}

#' 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 grid Evaluation backend. `"auto"` uses the fastest verified backend:
#'   FFT for Gaussian and Gaussian-copula kernels, and exact evaluation
#'   otherwise. `"exact"` always retains direct evaluation and `"fft"`
#'   explicitly requests the approximate backend supported by
#'   [kdensity_grid()].
#' @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, grid = "auto", ...) {
  plot_helper(
    x, range, plot_start, zero_line, grid = grid, ptype = "plot", ...
  )
}

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

#' @rdname plot.kdensity
#' @export
points.kdensity <- function(x, range = NULL, plot_start = FALSE,
                            zero_line = TRUE, grid = "auto", ...) {
  plot_helper(
    x, range, plot_start, zero_line, grid = grid,
    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 grid Evaluation backend passed to [kdensity_grid()].
#' @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, grid = "auto",
                        ptype = c("plot", "lines", "points"), ...) {
  ptype <- match.arg(ptype)
  grid <- match.arg(grid, c("exact", "auto", "fft"))

  supplied_range <- !is.null(range)
  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 <- utils::modifyList(
    defaults,
    supplied,
    keep.null = TRUE
  )
  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."
    assert_(!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 {
    if (grid == "exact") {
      args$y <- x(range)
    } else {
      grid_n <- length(range)
      if (supplied_range &&
          identical(x$kernel$compiled_id, 11L)) {
        grid_n <- max(2048, 2 * as.double(grid_n))
      }
      grid_fit <- kdensity_grid(
        x,
        from = min(range),
        to = max(range),
        n = grid_n,
        method = grid
      )
      if (supplied_range) {
        unique_grid <- !duplicated(grid_fit$x)
        args$y <- stats::approx(
          grid_fit$x[unique_grid],
          grid_fit$y[unique_grid],
          xout = range,
          rule = 2
        )$y
      } else {
        args$x <- grid_fit$x
        args$y <- grid_fit$y
      }
    }
  }

  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",
      vapply(seq_along(parameters), function(i) paste0(names(parameters)[i], ": ", formatC(parameters[i], digits), "\n"), character(1)),
      "\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)
}