R/kernelSmoothing.R

In PCpluS: Piecewise Constant Plus Smooth Regression

kernelSmoothing <- function(y, bandwidth, nbandwidth = 30, kernel = c("epanechnikov", "gaussian", "rectangular",
                                                                      "triangular", "biweight", "silverman")) {
  # TODO test y
  n <- length(y)
  
  if (missing(bandwidth)) {
    # TODO test nbandwidth
    
    bandwidth <- exp(seq(log(1 / n * 1.01), log(0.5), length.out = nbandwidth))
  }
  
  # TODO test bandwidth
  
  if (!is.function(kernel)) {
    kernel <- match.arg(kernel)
    
    kernel <- switch(kernel,
                     rectangular = function(x) 1 / 2,
                     triangular = function(x) 1 - abs(x),
                     epanechnikov = function(x) 3 / 4 * (1 - x^2),
                     biweight = function(x) 5 / 16 * (1 - x^2)^2,
                     gaussian = function(x) dnorm(x, 0, 1),
                     silverman = function(x) exp(-abs(x) / sqrt(2)) * sin(abs(x) / sqrt(2) + pi / 4) / 2,
                     stop("unknown kernel")
    )
  }
  
  if (length(bandwidth) > 1) {
    cv <- numeric(length(bandwidth))
    for (i in seq_along(bandwidth)) {
      b <- as.integer(n * bandwidth[i] + 1e-12)
      cv[i] <- .CVkernelSmoothing(Y = y, K = kernel((1:b) / (n * bandwidth[i])))
    }
    bandwidth <- bandwidth[which.min(cv)]
  }
  
  if (bandwidth == Inf) {
    return(rep(mean(y), length(y)))
  }

  b <- as.integer(n * bandwidth + 1e-12)
  .kernelSmoothing(y, kernel((-b:b) / (n * bandwidth)))
}