R/bandwidth.R

In qdfauction: Nonparametric Quantile Density Estimation and Structural Analysis of First-Price Auctions

# Bandwidth selection ----------------------------------------------------------
#
# bcv_criterion() and rlcv_criterion() are revised from R code by Y. P. Chaubey,
# I. Dewan and J. Li (see inst/COPYRIGHTS); wbcv_criterion() is from
# Doosti (2026b).

# internal: evaluate an estimator on the CV grid with and without observation i
cv_fits <- function(x, h, fun, u, loo = TRUE) {
  n  <- length(x)
  q  <- fun(u, x, h)
  q_loo <- if (loo) {
    if (identical(fun, qdf_bernstein)) qdf_bernstein_loo(u, x, h)       # exact, O(n K)
    else if (isFALSE(attr(fun, "bernstein_anchor")))
      qdf_bernstein_loo(u, x, h, zero_at_origin = FALSE)                # exact, O(n K)
    else vapply(seq_len(n), function(i) fun(u[i], x[-i], h), numeric(1))
  } else {
    # Scaling approximation q_{-i} ~ q (n-1)/n. NOT valid for estimators built
    # from order statistics (removing x_(i) shifts the quantile positions);
    # with it, least-squares criteria select the roughest fit on the grid.
    # Kept only for backward compatibility; the default is exact LOO.
    q * (n - 1) / n
  }
  list(q = q, q_loo = q_loo)
}

#' Cross-validation criteria for quantile density estimators
#'
#' @description
#' Three data-driven criteria for the smoothing parameter of a quantile density
#' estimator, written as functions of \eqn{m = 1/h} so that they can be passed
#' to [stats::optimize()]. Each is evaluated on the grid \eqn{u_i = i/n} with
#' spacings \eqn{\Delta_i = X_{(i)} - X_{(i-1)}}, \eqn{X_{(0)} = 0}, so that
#' \eqn{\sum_i \Delta_i\, \phi(\hat q(u_i))} approximates \eqn{\int \phi\{1/\hat g(x)\}\,dx}.
#'
#' * `bcv_criterion()` is the biased (least-squares) cross-validation used by
#'   Chaubey, Dewan and Li (2024) and Doosti, Dewan and Talebian (2025):
#'   \deqn{\mathrm{BCV}(h) = \sum_i \Delta_i\, \hat q(u_i;h)^{-2} - \frac{2}{n}\sum_i \hat q_{-i}(u_i;h)^{-1},}
#'   i.e. the classical \eqn{\int \hat g^2 - 2\,\bar{\hat g}_{-i}} written in
#'   terms of \eqn{\hat q = 1/\hat g}. (Equation (11) of Doosti, Dewan and
#'   Talebian (2025) prints \eqn{\hat q^2} in the first term; the simulation
#'   code, reproduced here, uses \eqn{\hat q^{-2}}.)
#' * `rlcv_criterion()` is the robust likelihood cross-validation of Wu (2019),
#'   \eqn{-\{ \frac1n \sum_i \ln^*(\hat q_{-i}(u_i)^{-1}) - b^*(h) \}}, with
#'   \eqn{\ln^*} linearised below the threshold \eqn{a_n}. It is written so that
#'   it is *minimised*.
#' * `wbcv_criterion()` is the weighted criterion of Doosti (2026b), which
#'   targets \eqn{u\,q(u)} (the quantity entering the first-price auction
#'   inversion) rather than \eqn{q(u)}:
#'   \deqn{\mathrm{WBCV}(h) = \sum_i \Delta_i\{u_i \hat q(u_i)\}^2 - \frac{2}{n}\sum_i u_i / \hat q_{-i}(u_i).}
#'   **Not recommended; kept only for backward compatibility.** The criterion
#'   is not scale-invariant: its two terms have different units (the first scales
#'   like the cube of the data units, the second like their inverse), so the
#'   selected bandwidth depends on the units of the data; on data of order one
#'   or larger the first term dominates and the criterion tends to select the
#'   smoothest fit in `m_range`. The criterion is being revised (Doosti 2026b,
#'   revised version); `wbcv_criterion()` and `bandwidth = "wbcv"` issue a
#'   warning and are kept only so that existing code runs. Use BCV, RLCV or a
#'   numeric bandwidth instead.
#'
#' @param m inverse bandwidth, \eqn{m = 1/h}.
#' @param x numeric sample (sorted internally).
#' @param fun an estimator with signature `fun(u, x, h)`, e.g. [qdf_kernel()].
#'   Extra arguments (such as `H` for [qdf_soni()]) should be bound with a
#'   closure before calling.
#' @param loo logical; if `TRUE` (default) leave-one-out fits are computed
#'   exactly (by refitting without observation \eqn{i}, or by the index shift of
#'   [qdf_bernstein_loo()] for the Bernstein estimator). If `FALSE` the scaling
#'   approximation \eqn{\hat q_{-i} \approx \hat q\,(n-1)/n} is used. The
#'   approximation is not valid for estimators built from order statistics:
#'   removing \eqn{X_{(i)}} shifts every higher quantile position, and with it
#'   least-squares criteria select the roughest fit in `m_range`. It is kept
#'   only for backward compatibility.
#' @param u_scheme grid on which the criterion is evaluated: `"i/n"` (the
#'   published convention) or `"i/(n+1)"`.
#' @return a single numeric value; smaller is better.
#' @references
#' Wu, X. (2019). Robust likelihood cross validation for kernel density
#' estimation. *Journal of Business & Economic Statistics*, 37, 761--770.
#' @seealso [select_bandwidth()], [qdf()]
#' @examples
#' x <- rgamma(60, 5, 1)
#' bcv_criterion(20, x, qdf_kernel)
#' rlcv_criterion(20, x, qdf_kernel)
#' optimize(bcv_criterion, c(10, 50), x = x, fun = qdf_bernstein)$minimum
#' @export
bcv_criterion <- function(m, x, fun, loo = TRUE, u_scheme = c("i/n", "i/(n+1)")) {
  x <- prep_sample(x); n <- length(x); h <- 1 / m
  u_scheme <- match.arg(u_scheme)
  u <- default_u(n, u_scheme)
  d <- diff(c(0, x))
  f <- cv_fits(x, h, fun, u, loo)
  sum(d / f$q^2) - 2 * mean(1 / f$q_loo)
}

#' @rdname bcv_criterion
#' @export
rlcv_criterion <- function(m, x, fun, loo = TRUE, u_scheme = c("i/n", "i/(n+1)")) {
  x <- prep_sample(x); n <- length(x); h <- 1 / m
  u_scheme <- match.arg(u_scheme)
  u <- default_u(n, u_scheme)
  d <- diff(c(0, x))
  f <- cv_fits(x, h, fun, u, loo)
  a <- a_n(x)
  g     <- 1 / f$q                    # implied density along the quantile grid
  g_loo <- 1 / f$q_loo
  ll    <- mean(log_star(g_loo, a))
  bstar <- sum(d * g * (g >= a)) + sum(d * g^2 * (g < a)) / (2 * a)
  -(ll - bstar)
}

#' @rdname bcv_criterion
#' @export
wbcv_criterion <- function(m, x, fun, loo = TRUE, u_scheme = c("i/(n+1)", "i/n")) {
  .wbcv_warning()
  .wbcv_criterion(m, x, fun, loo = loo, u_scheme = u_scheme)
}

.wbcv_warning <- function()
  warning("The WBCV criterion depends on the units of the data and tends to select the ",
          "smoothest fit; it is being revised (Doosti 2026b). Use BCV, RLCV or a numeric ",
          "bandwidth instead.", call. = FALSE)

# WBCV criterion without the warning (used inside select_bandwidth(), which warns once)
.wbcv_criterion <- function(m, x, fun, loo = TRUE, u_scheme = c("i/(n+1)", "i/n")) {
  x <- prep_sample(x); n <- length(x); h <- 1 / m
  u_scheme <- match.arg(u_scheme)
  u <- default_u(n, u_scheme)
  d <- diff(c(0, x))
  f <- cv_fits(x, h, fun, u, loo)
  q <- pmax(f$q, 1e-10); q_loo <- pmax(f$q_loo, 1e-10)
  sum(d * (u * q)^2) - 2 * mean(u / q_loo)
}

#' Select the smoothing parameter of a quantile density estimator
#'
#' Minimises one of the criteria in [bcv_criterion()] over \eqn{m = 1/h} in
#' `m_range` with [stats::optimize()], optionally after a coarse grid search.
#'
#' @inheritParams bcv_criterion
#' @param selector `"bcv"`, `"rlcv"` or `"wbcv"`.
#' @param m_range numeric length-2 range for \eqn{m = 1/h}.
#' @param grid search method: `0` (the default, as in the published
#'   simulations) runs [stats::optimize()] directly over `m_range`; an integer
#'   `grid > 0` first evaluates the criterion on that many equally spaced
#'   points and then runs `optimize()` in the bracket around the grid minimum;
#'   `"exhaustive"` evaluates the criterion at every integer `m` in `m_range`
#'   and returns the global minimiser. When a criterion is multimodal in `m`,
#'   `optimize()` can return an interior local minimum, while the global
#'   minimiser can lie at an end of `m_range`; boundary solutions are flagged
#'   (see below). `NULL` is the same as `0`.
#' @param optimizer optional function `optimizer(f, interval)` used instead of
#'   [stats::optimize()]; it must return a list with elements `minimum` and
#'   `objective`. Useful to plug in a global optimiser for criteria with several
#'   local minima, e.g. `function(f, interval) \{ g <- seq(interval[1],
#'   interval[2], length.out = 200); v <- sapply(g, f); list(minimum = g[which.min(v)],
#'   objective = min(v)) \}` for an exhaustive grid search. When `optimizer`
#'   is supplied it is used, and `grid = "exhaustive"` is treated as `grid = 0`.
#' @return a list with `h`, `m`, `selector`, `criterion` (the minimised value)
#'   and `m_range`; `at_boundary` (logical) and `boundary` (`"lower"`,
#'   `"upper"` or `NA`), which record a minimiser at an end of `m_range`; with
#'   `grid = "exhaustive"` also the evaluated grid `grid_m` and criterion values
#'   `grid_cv`.
#' @section Boundary solutions:
#' If the minimiser lies at either end of `m_range`, whatever the search
#' method, the criterion did not select a bandwidth in that range. A warning is
#' issued and `at_boundary = TRUE` is returned; widen `m_range` or give the
#' bandwidth directly.
#' @examples
#' x <- rgamma(60, 5, 1)
#' select_bandwidth(x, qdf_poisson, "bcv", m_range = c(10, 50))
#' @export
select_bandwidth <- function(x, fun, selector = c("bcv", "rlcv", "wbcv"),
                             m_range = c(10, 50), loo = NULL, grid = NULL,
                             u_scheme = NULL, optimizer = NULL) {
  selector <- match.arg(selector)
  crit <- switch(selector, bcv = bcv_criterion, rlcv = rlcv_criterion,
                 wbcv = wbcv_criterion)
  if (is.null(loo)) loo <- TRUE
  if (is.null(u_scheme)) u_scheme <- if (selector == "wbcv") "i/(n+1)" else "i/n"
  # Default search: optimize() over m_range, as in the published simulations.
  if (is.null(grid)) grid <- 0L
  if (selector == "wbcv") {
    .wbcv_warning()
    crit <- .wbcv_criterion
  }
  if (!(identical(grid, "exhaustive") ||
        (is.numeric(grid) && length(grid) == 1L && is.finite(grid) && grid >= 0 &&
         grid == round(grid))))
    stop("`grid` must be \"exhaustive\" or a single non-negative integer (0 runs ",
         "stats::optimize() over `m_range`).", call. = FALSE)
  obj <- function(m) {
    v <- tryCatch(crit(m, x, fun, loo = loo, u_scheme = u_scheme),
                  error = function(e) NA_real_)
    if (!is.finite(v)) Inf else v
  }
  lo <- m_range[1]; hi <- m_range[2]
  # a user-supplied optimizer takes precedence over the exhaustive search
  if (identical(grid, "exhaustive") && !is.null(optimizer)) grid <- 0L
  if (identical(grid, "exhaustive")) {
    # The criteria can be multimodal in m; stats::optimize() then returns an
    # interior local minimum. Evaluate every integer m.
    mg <- seq(ceiling(lo), floor(hi), by = 1L)
    cv <- vapply(mg, obj, numeric(1))
    if (!any(is.finite(cv))) stop("Criterion is non-finite over the whole range.", call. = FALSE)
    k  <- which.min(cv)
    out <- list(h = 1 / mg[k], m = mg[k], selector = selector, criterion = cv[k],
                m_range = m_range, grid_m = mg, grid_cv = cv)
    return(.flag_boundary(out, at_lo = k == 1L, at_hi = k == length(mg)))
  }
  if (is.numeric(grid) && grid > 0L) {
    mg <- seq(lo, hi, length.out = grid)
    cv <- vapply(mg, obj, numeric(1))
    if (!any(is.finite(cv))) stop("Criterion is non-finite over the whole range.", call. = FALSE)
    k  <- which.min(cv)
    lo <- mg[max(1L, k - 1L)]; hi <- mg[min(grid, k + 1L)]
    if (hi <= lo) { lo <- m_range[1]; hi <- m_range[2] }
  }
  o <- if (is.null(optimizer)) stats::optimize(obj, c(lo, hi)) else {
    oo <- optimizer(obj, c(lo, hi))
    if (!is.list(oo) || is.null(oo$minimum) || is.null(oo$objective))
      stop("`optimizer` must return a list with elements `minimum` and `objective`.", call. = FALSE)
    oo
  }
  if (!is.finite(o$objective))
    stop("Bandwidth selection failed: criterion is non-finite. ",
         "Try a different `m_range` or estimator.", call. = FALSE)
  out <- list(h = 1 / o$minimum, m = o$minimum, selector = selector,
              criterion = o$objective, m_range = m_range)
  # continuous search: optimize() approaches an end of the interval but never
  # returns it exactly, so allow a small tolerance
  tol <- 1e-3 * diff(m_range)
  .flag_boundary(out, at_lo = o$minimum <= m_range[1] + tol,
                 at_hi = o$minimum >= m_range[2] - tol)
}

# A minimiser at an end of m_range means the criterion did not select a
# bandwidth in that range: record it and warn rather than return the endpoint
# silently.
.flag_boundary <- function(out, at_lo, at_hi) {
  out$at_boundary <- at_lo || at_hi
  out$boundary <- if (at_lo) "lower" else if (at_hi) "upper" else NA_character_
  if (out$at_boundary)
    warning(sprintf(paste0("The %s criterion is minimised at the %s end of m_range = [%s, %s] ",
                           "(m = %s): it did not select a bandwidth in this range. ",
                           "Widen `m_range` or give the bandwidth directly."),
                    toupper(out$selector), out$boundary, format(out$m_range[1]),
                    format(out$m_range[2]), format(signif(out$m, 4))), call. = FALSE)
  out
}