R/get_crit_approx.R

In smnet: Smoothing for Stream Network Data

#' @importFrom spam update.spam.chol.NgPeyton
#' @importFrom spam backsolve.spam
#' @importFrom spam forwardsolve.spam
#' @importFrom stats rnorm

get_crit_approx<-function(rhoArgs, X, XTX, P, response, cholFactor, n, 
                          np = nrow(XTX), Xw, Xy, n.sm, identifyBit, crit){ 
  
  for(j in 1:n.sm) P[[j]]<-P[[j]]*exp(rhoArgs[j])
  Psum       <- Reduce("+", P)
  info       <- XTX + Psum + identifyBit
  U        <- tryCatch(update.spam.chol.NgPeyton(cholFactor, info), error=function(e) e, warning=function(w) w)
  if(("warning" %in% class(U)) ){ #| any(abs(rhoArgs) > 20)){
    out <- rnorm(n = 1, mean = 10^20)
  } else {
    beta_hat   <- backsolve.spam(U, forwardsolve.spam(U, Xy)) 
    left1      <- forwardsolve.spam(U, Xw)
    ED         <- ifelse(!is.null(dim(left1)), sum(rowMeans(left1*left1)), sum(left1*left1))
    fit        <- X %*% beta_hat  
    if(crit == "AICC") out <- log(sum((response - fit)^2)) + (2*(ED+1)/(n-ED-2))
    if(crit == "AIC")  out <- log(sum((response - fit)^2)) + 2*ED/n
    if(crit == "GCV")  out <- sum((response - fit)^2)/(1-(ED/n))^2
  }
  out
}