R/gaussian_kde.R

In BDIB: Bayesian Dynamic Information Borrowing for Longitudinal Data with Informative Dropout

#' Gaussian Kernel Density Estimation
#' @param samples Numeric vector of samples.
#' @param eval_points Numeric vector of locations to evaluate.
#' @param bandwidth Optional bandwidth; Silverman's rule if \code{NULL}.
#' @return Numeric vector of density estimates.
#' @export
#' @examples
#' set.seed(1000)
#' xs <- rnorm(1000, mean = -3, sd = 0.5)
#' pts <- seq(-4.5, -1.5, length.out = 5)
#' dens <- gaussian_kde(xs, pts)      # KDE values at 'pts'
#' cbind(pts, round(dens, 4))
gaussian_kde <- function(samples, eval_points, bandwidth = NULL) {
  stopifnot(is.numeric(samples), length(samples) > 1)
  stopifnot(is.numeric(eval_points))
  if (is.null(bandwidth)) {
    n <- length(samples)
    sigma <- stats::sd(samples)
    sigma_robust <- stats::IQR(samples) / 1.34
    bandwidth <- 0.9 * min(sigma, sigma_robust) * n^(-1/5)
    bandwidth <- max(bandwidth, .Machine$double.eps)
  }
  z <- outer(eval_points, samples, FUN = "-") / bandwidth
  rowMeans(stats::dnorm(z)) / bandwidth
}