R/4_prediction_SF.R

In funGp: Gaussian Process Models for Scalar and Functional Inputs

#' @author José Betancourt, François Bachoc, Thierry Klein and Jérémy Rohmer
makePreds_SF <- function(sMs.tp, sMs.pp, fMs.tp, fMs.pp, sig2, thetas_s, thetas_f, kerType,
                         L, LInvY, detail, nugget){
  # create empty prediction list
  preds <- list()

  # compute and store conditional mean and standard deviation
  K.tp <- sig2 * setR(thetas_s, sMs.tp, kerType) * setR(thetas_f, fMs.tp, kerType)
  LInvK <- backsolve(L, K.tp, upper.tri = FALSE)
  preds$mean <- t(LInvK) %*% LInvY
  preds$sd <- sqrt(pmax(sig2 - apply(LInvK, 2, crossprod), 0))

  # if user requires details, provide K.tp and K.pp
  if (detail == "full") {
    preds$K.tp <- K.tp
    R <- setR(thetas_s, sMs.pp, kerType) * setR(thetas_f, fMs.pp, kerType)
    preds$K.pp <- sig2 * (R + diag(nugget, nrow = nrow(R), ncol = ncol(R)))
  }

  return(preds)
}

#' @author José Betancourt, François Bachoc, Thierry Klein and Jérémy Rohmer
preMats_SF <- function(sMs, fMs, sOut, sig2, thetas_s, thetas_f, kerType, nugget){
  # precompute L and LInvY matrices
  R <- setR(thetas_s, sMs, kerType) * setR(thetas_f, fMs, kerType)
  K.tt <- sig2 * (R + diag(nugget, nrow = nrow(R), ncol = ncol(R)))
  L <- t(chol(K.tt))
  LInvY <- backsolve(L, sOut, upper.tri = FALSE)

  return(list(L = L, LInvY = LInvY))
}