R/getLikelihood.R

In dominikkirchhoff/CCKriging: Kriging for Mixed Continuous and Categorical Input Variables

#' @title Compute Maximum Likelihood Value
#'
#' @description
#'   Compute the Maximum Likelihood value for CCKriging.
#'
#' @param model [\code{\link[=makeCCKriging]{CCKriging}}]\cr
#'   The \code{\link[=makeCCKriging]{CCKriging}} model the likelihood should be computed for.
#' @param par [\code{numeric}]\cr
#'   Optional vector of parameters of the Kriging model. The first values belong to the
#'   continuous variables, the rest to the categorical variables. Defaults to \code{NULL}, which
#'   means the parameters of the \code{model} object are used.
#' @param info [\code{logical(1)}]\cr
#'   Should extra information be printed to the console? Default is \code{TRUE}.
#' @export
getLikelihood = function(model, par = NULL, config, info = TRUE) {
  x = model$x
  y = model$y
  if (is.null(par))
    par = model$par
  cat.type = model$config$cat.type
  cont.type = model$config$cont.type
  nu = model$nu
  n = nrow(x)
  d = ncol(x)
  ones = rep(1, n)

  qual.inds = which(sapply(x, class) == "factor")
  quan.inds = setdiff(1:d, qual.inds)

  corr.matrix = getCorrMatrix(model, par, config = config, info = info)
  chol.status = try(corr.matrix.chol <- chol(corr.matrix), silent = TRUE)
  chol.error = ifelse(class(chol.status) == "try-error", TRUE, FALSE)
  if (chol.error) {
    corr.matrix = repairCorrMat(corr.matrix, info = info)
    corr.matrix.chol = chol(corr.matrix)
  }

  x.dice = backsolve(t(corr.matrix.chol), y, upper.tri = FALSE)
  M.dice = backsolve(t(corr.matrix.chol), ones, upper.tri = FALSE)
  ## compute z.dice in order to compute sigma.hat.sq
  Q = qr.Q(qr(M.dice))
  H = Q %*% t(Q)
  z.dice = as.numeric(x.dice - H %*% x.dice)

  sigma.hat.sq = drop(crossprod(z.dice)/length(z.dice))
  targetFun = n * log(sigma.hat.sq) + 2 * sum(log(diag(corr.matrix.chol)))
  return(targetFun)
}