R/jsd_continuous_ci.R

In jsdtools: Jensen-Shannon Divergence Estimation, Confidence Intervals, and Distribution Plots

#' Bootstrap confidence interval for continuous JSD
#'
#' @param x Numeric vector for group 1.
#' @param y Numeric vector for group 2.
#' @param B Number of bootstrap replicates.
#' @param conf_level Confidence level. Defaults to 0.95.
#' @param L Optional lower integration bound.
#' @param U Optional upper integration bound.
#' @param base Logarithm base. Defaults to 2. Use `exp(1)` for nats.
#' @param bw Bandwidth passed to [stats::density()].
#' @param kernel Kernel passed to [stats::density()].
#' @param grid_n Number of grid points used for KDE.
#' @param qrange Quantile range used when `L` and `U` are not supplied.
#' @param extend Extension multiplier for the automatically chosen range.
#' @param eps Small constant for numerical stability.
#' @param renormalize Logical; renormalize estimated densities over the grid?
#' @param seed Optional random seed.
#' @param na_rm Logical; remove missing values?
#' @param na_rm_failed Logical; drop failed bootstrap draws when summarizing?
#'
#' @return An object of class `"jsd_ci"`.
#' @export
jsd_continuous_ci <- function(x, y,
                              B = 1000,
                              conf_level = 0.95,
                              L = NULL,
                              U = NULL,
                              base = 2,
                              bw = "nrd0",
                              kernel = "gaussian",
                              grid_n = 4096,
                              qrange = c(0.001, 0.999),
                              extend = 3,
                              eps = 1e-12,
                              renormalize = TRUE,
                              seed = NULL,
                              na_rm = TRUE,
                              na_rm_failed = TRUE) {
  check_base(base)
  check_conf_level(conf_level)

  cleaned <- validate_xy(x, y, min_n = 2, na_rm = na_rm, finite_only = TRUE)
  x <- cleaned$x
  y <- cleaned$y

  if (!is.null(seed)) {
    set.seed(seed)
  }

  if (is.null(L) || is.null(U)) {
    ru <- fixed_range(x, y, qrange = qrange, extend = extend)
    L <- unname(ru["L"])
    U <- unname(ru["U"])
  }

  est_obj <- jsd_continuous(
    x, y,
    L = L, U = U,
    base = base,
    bw = bw,
    kernel = kernel,
    grid_n = grid_n,
    qrange = qrange,
    extend = extend,
    eps = eps,
    renormalize = renormalize,
    na_rm = FALSE
  )

  estimate <- est_obj$estimate
  n <- length(x)
  m <- length(y)

  boot_values <- rep(NA_real_, B)

  for (b in seq_len(B)) {
    xb <- sample(x, size = n, replace = TRUE)
    yb <- sample(y, size = m, replace = TRUE)

    boot_values[b] <- tryCatch(
      jsd_continuous(
        xb, yb,
        L = L, U = U,
        base = base,
        bw = bw,
        kernel = kernel,
        grid_n = grid_n,
        qrange = qrange,
        extend = extend,
        eps = eps,
        renormalize = renormalize,
        na_rm = FALSE
      )$estimate,
      error = function(e) NA_real_
    )
  }

  boot_valid <- if (na_rm_failed) boot_values[is.finite(boot_values)] else boot_values

  if (sum(is.finite(boot_valid)) < 10) {
    warning("Too few valid bootstrap replicates. Check data or bandwidth settings.")
  }

  alpha <- 1 - conf_level
  conf_int <- stats::quantile(
    boot_valid,
    probs = c(alpha / 2, 1 - alpha / 2),
    na.rm = TRUE,
    names = FALSE
  )

  out <- list(
    estimate = estimate,
    conf_int = c(lower = conf_int[1], upper = conf_int[2]),
    conf_level = conf_level,
    type = "continuous",
    method = "KDE bootstrap",
    base = base,
    n_x = length(x),
    n_y = length(y),
    boot_values = boot_values,
    boot_valid_n = sum(is.finite(boot_values)),
    boot_mean = mean(boot_valid, na.rm = TRUE),
    boot_median = stats::median(boot_valid, na.rm = TRUE),
    boot_se = stats::sd(boot_valid, na.rm = TRUE),
    boot_bias = mean(boot_valid, na.rm = TRUE) - estimate,
    range = c(L = L, U = U),
    settings = list(
      B = B,
      bw = bw,
      kernel = kernel,
      grid_n = grid_n,
      qrange = qrange,
      extend = extend,
      eps = eps,
      renormalize = renormalize
    )
  )

  class(out) <- c("jsd_ci", "list")
  out
}