R/functionFactory_Neg2LogLik.R

In GeneralizedWendland: Fully Parameterized Generalized Wendland Covariance Function

neg2loglikFactory <- function(y, X = data.frame(), distmat, covariance = NULL, cov.args = list(),
                              chol.args = list(), Rstruct = NULL, covarianceFunction = NULL,
                              choleskyFunction = NULL){

  ######################################################################
  # Checks on input, and generating covarianceFunction and choleskyFunction
  # if they are not already provided as input
  #----------------------------------------------------------------------

  if (is.null(covarianceFunction)) {

    if (is.null(covariance)) {
      stop("Must provide argument covariance or an object ",
           "generated by covarianceFactory.")
    }

    covarianceFunction <- covarianceFactory(covariance = covariance, cov.args = cov.args)
  }

  if (is.null(choleskyFunction)) {
    choleskyFunction <- choleskyFactory(chol.args = chol.args, Rstruct = Rstruct)
  }


  ######################################################################
  # Function which computes the neg2loglikelihood based
  # on input parameters
  #----------------------------------------------------------------------

  neg2loglikFunction <- function(parameters){

    ######################################################################
    # Setting up optional covariates and separating parameters into
    # beta and theta
    #----------------------------------------------------------------------

    if (ncol(X) == 0) {
      theta     <- parameters
      residuals <- y

    } else {
      beta      <- parameters[seq_len(ncol(X))]
      theta     <- parameters[-seq_len(ncol(X))]
      residuals <- y - X %*% beta
    }


    ######################################################################
    # log-likelihood
    #----------------------------------------------------------------------

    n          <- length(residuals)
    Sigma      <- covarianceFunction(h = distmat, theta = theta)
    cholS      <- choleskyFunction(Sigma)

    ######################################################################
    # Handle sparse and dense computations
    # [Note: Considering adding functionFactory or an alternative approach
    #  to eliminate this case separation.]
    #----------------------------------------------------------------------

    if (is.spam(Sigma)) {
      logdet   <- c(determinant.spam.chol.NgPeyton(cholS)$modulus)
      quadform <- t(residuals) %*% solve.spam(cholS, residuals)

    } else {
      piv <- attr(cholS, which = "pivot")
      if (is.null(piv)) piv <- 1:n

      logdet   <- sum(log(diag(cholS)))
      quadform <- t(residuals[piv]) %*%
        backsolve(cholS, forwardsolve(cholS, residuals[piv], upper.tri = TRUE, transpose = TRUE), n)
    }

    ######################################################################
    # Compute neg2loglikelihood
    #----------------------------------------------------------------------

    n2ll <- c(n * log(2*pi) + 2*logdet + quadform)
    return(n2ll)
  }

  return(neg2loglikFunction)
}