R/pelt.R

In PCpluS: Piecewise Constant Plus Smooth Regression

# restimate cps by PELT
.pelt <- function(obs, sd = NULL) {
  # little hack since PELT does not scale
  if (is.null(sd)) {
    sd <- as.numeric(diff(quantile(diff(obs), c(0.25, 0.75))) / diff(qnorm(c(0.25, 0.75))) / sqrt(2))
  }
  
  pen.value <- 2 * log(length(obs[6:(length(obs) - 5)])) 
  # pelt <- fpop::Fpop(obs[6:(length(obs) - 5)] / sd, lambda = pen.value)
  
  pelt <- changepoint::cpt.mean(data = obs[6:(length(obs) - 5)] / sd, penalty = "Manual",
                                pen.value = pen.value, method = "PELT")
  
  left <- c(1, pelt@cpts[-length(pelt@cpts)] + 6)
  right <- c(left[-1] - 1, length(obs))
  values <- pelt@param.est$mean
  est <- numeric(length(obs))
  for (segment in seq_along(left)) {
    est[left[segment]:right[segment]] <- values[segment]
  }
  est <- est * sd
  
  list(est = est, cps = pelt@cpts[-length(pelt@cpts)] + 6L)
}