R/pcplus.R

In PCpluS: Piecewise Constant Plus Smooth Regression

pcplus <- function(y, bandwidth, lambda, sd = NULL, nlambda = 30L, thresh = 1e-7, maxit = 1e5L) {
  if (!is.numeric(y)) {
    stop("'y' must be a numeric vector")
  }
  
  if (!all(is.finite(y))) {
    stop("'y' must contain only finite values")
  }
  
  n <- length(y)
  
  if (!is.numeric(bandwidth) || length(bandwidth) != 1) {
    stop("'bandwidth' must be a single numeric")
  }
  
  if (!is.finite(bandwidth)) {
    test <- bandwidth == Inf
    if (is.na(test) || !test) {
      stop("'bandwidth' must be finite or Inf")
    }
  }
  
  if (bandwidth <= 0) {
    stop("'bandwidth' must be a positive value")
  }
  
  if (bandwidth > 0.25 && bandwidth != Inf) {
    bandwidth <- Inf
    warning("'bandwidth' larger than 0.25 has been replaced by Inf")
  }
  
  if (is.null(sd)) {
    sd <- as.numeric(diff(quantile(diff(y), c(0.25, 0.75))) / diff(qnorm(c(0.25, 0.75))) / sqrt(2))
  } else {
    if (!is.numeric(sd) || length(sd) != 1 || !is.finite(sd)) {
      stop("'sd' must be a single finite numeric value or NULL")
    }
    
    if (sd <= 0) {
      stop("'sd' must be a positive value or NULL")
    }
  }
  
  if (bandwidth == Inf) {
    if (length(lambda) == 1 && !is.na(lambda) && lambda == Inf) {
      return(list(est = rep(mean(y), n), smooth = rep(mean(y), n), cpfun = rep(0, n), cps = integer(0)))
    }
    
    pelt <- .pelt(obs = y - mean(y), sd = sd)
    return(list(est = pelt$est + mean(y), smooth = rep(mean(y), n), cpfun = pelt$est, cps = pelt$cps))
  }
  
  if (!is.numeric(lambda)) {
    stop("'lambda' must be a numeric")
  }
  
  if (!all(is.finite(lambda))) {
    test <- lambda[!is.finite(lambda)] == Inf
    if (any(is.na(test)) | any(!test)) {
      stop("'lambda' must be finite or Inf")
    }
  }
  
  if (any(lambda <= 0)) {
    stop("'lambda' must contain only positive values")
  }
  
  L <- as.integer(bandwidth * n + 1e-14)
  if (L < 1L) {
    L <- 1L
    bandwidth <- 1 / n
    warning("'bandwidth' smaller than 1 / n has been replaced by 1/ n")
  }
  
  if (!is.numeric(thresh) || length(thresh) != 1 || !is.finite(thresh)) {
    stop("'thresh' must be a single finite numeric value")
  }
  
  if (thresh <= 0) {
    stop("'thresh' must be a positive value")
  }
  
  if (!is.numeric(maxit) || length(maxit) != 1 || !is.finite(maxit)) {
    stop("'maxit' should be a single finite integer.")
  }
  
  if (maxit <= 0) {
    stop("'maxit' must be a positive integer")
  }
  
  if (!is.integer(maxit)) {
    maxit <- as.integer(maxit + 1e-6)
  }
  
  ImSy <- y - .kernelSmoothingEpanechnikov(y, bandwidth)
  null_dev <- mean(ImSy^2)
  
  listKernel <- list(cusumKernel = numeric(2 * L + 1),
                     Xty = numeric(n - 1),
                     XtX = matrix(0, 2 * L - 1, 4 * L - 2),
                     isComputedXtX = matrix(FALSE, 2 * L - 1, 4 * L - 2), 
                     XtXgap = numeric(2 * L),
                     ImSX = matrix(0, 3 * L, 2 * L),
                     isComputedImSX = logical(2 * L))
  
  .initializeKernel(ImSy, bandwidth, listKernel$cusumKernel, listKernel$Xty, listKernel$XtX, listKernel$isComputedXtX,
                    listKernel$XtXgap, listKernel$ImSX, listKernel$isComputedImSX)
  
  if (length(lambda) == 1) {
    if (lambda == Inf) {
      smooth <- .kernelSmoothingEpanechnikov(y, bandwidth)
      return(list(est = smooth, smooth = smooth, cpfun = rep(0, n), cps = integer(0)))
    }
    
    
    if (!is.numeric(nlambda) || length(nlambda) != 1 || !is.finite(nlambda) || nlambda <= 0) {
      stop("nlambda must be a single positive integer")
    }

    if (nlambda > 1) {
      Xty <- .getXty(ImSy, bandwidth, listKernel$cusumKernel, listKernel$Xty, listKernel$XtX, listKernel$isComputedXtX,
                     listKernel$XtXgap, listKernel$ImSX, listKernel$isComputedImSX)
      
      lambda_max <- max(abs(Xty))
      lambda <- exp(seq(from = log(lambda_max), to = log(lambda), length.out = nlambda))
    }
  } else {
    if (any(diff(lambda) > 0)) {
      warning("lambda is not a decreasing sequence; we continue with the sorted sequence")
      lambda <- sort(lambda, decreasing = TRUE)
    }
  }

  beta <- .lassoImSX(ImSy, bandwidth, listKernel$cusumKernel, listKernel$Xty, listKernel$XtX, listKernel$isComputedXtX,
                     listKernel$XtXgap, listKernel$ImSX, listKernel$isComputedImSX,
                     lambda, thresh * null_dev, maxit)

  fit <- c(0, cumsum(beta[, length(lambda)]))
  smooth <- .kernelSmoothingEpanechnikov(y - fit, bandwidth)
  
  # restimate cps by PELT
  pelt <- .pelt(obs = as.numeric(y - smooth), sd = sd)
  
  # reestimate without penalty but condition on cps of PELT's estimate
  cps <- pelt$cps
  
  if (length(cps) > 0) {
    fitCps <- .postProcessing(cps - 2, ImSy, bandwidth, # - 1 for C++ style, - 1 for dealing with extended ImSX matrix, ImSX <- cbind(rep(0, n), ImSX)
                              listKernel$cusumKernel, listKernel$Xty, listKernel$XtX, listKernel$isComputedXtX,
                              listKernel$XtXgap, listKernel$ImSX, listKernel$isComputedImSX) 
    fit <- numeric(n)
    fit[cps] <- fitCps
    fit <- cumsum(fit)
    smooth <- .kernelSmoothingEpanechnikov(y - fit, bandwidth)
  } else {
    fit <- rep(0, n)
    smooth <- .kernelSmoothingEpanechnikov(y, bandwidth)
  }
  
  list(est = fit + smooth, smooth = smooth, cpfun = fit, cps = cps)
}