R/evd_example.R

In paulnorthrop/oola: Objected-oriented Loglikelihood Adjustment

# ================================ evd::fgev ================================ #

# We create S3 method functions for objects inheriting class "gev":

# loglikVec.gev()
# estfun.gev()
# coef.gev()
# nobs.gev()

#' @export
logLikVec.gev <- function(object, pars = NULL, ...) {
  if (!missing(...)) {
    warning("extra arguments discarded")
  }
  # Extract from object all the parameter estimates: free and fixed
  all_pars <- coef(object, complete = TRUE)
  free_pars <- coef(object, complete = FALSE)
  # If pars is supplied then overwrite the values of the free parameters
  if (!is.null(pars)) {
    names_all_pars <- names(all_pars)
    names_free_pars <- names(free_pars)
    which_free <- which(all_pars %in% free_pars)
    all_pars[which_free] <- pars
  }
  mu <- all_pars["loc"]
  sigma <- all_pars["scale"]
  xi <- all_pars["shape"]
  if (!is.null(object$nsloc)) {
    mu_reg <- all_pars[paste0("loc", names(object$nsloc))]
    mu <- mu + as.matrix(object$nsloc) %*% mu_reg
  }
  # Calculate the weighted loglikelihood contributions
  if (sigma <= 0) {
    val <- -Inf
  } else {
    val <- revdbayes::dgev(object$data, loc = as.vector(mu), scale = sigma,
                           shape = xi, log = TRUE)
  }
  # Return the usual attributes for a "logLik" object
  attr(val, "nobs") <- nobs(object)
  attr(val, "df") <- length(free_pars)
  class(val) <- "logLikVec"
  return(val)
}

#' @export
nobs.gev <- function(object, ...) {
  return(object$n)
}

#' @export
coef.gev <- function(object, complete = FALSE, ...) {
  if (complete) {
    val <- object$param
  } else {
    val <- object$estimate
  }
  return(val)
}

#' @export
vcov.gev <- function(object, complete = FALSE, ...) {
  # Variance-covariance matrix for the free parameters
  vc <- object$var.cov
  free_pars <- names(coef(object, complete = FALSE))
  dimnames(vc) <- list(free_pars, free_pars)
  if (complete) {
    all_pars <- names(coef(object, complete = TRUE))
    np <- length(all_pars)
    dummy <- matrix(NA, np, np)
    which_free <- which(all_pars %in% free_pars)
    dummy[which_free, which_free] <- vc
    vc <- dummy
    dimnames(vc) <- list(all_pars, all_pars)
  }
  return(vc)
}