Nothing
#' Information Criteria
#'
#' Compute AIC and (approximate) DIC for \code{evmOpt} objects
#'
#' @aliases AIC.evmSim
#' @param object fit model object
#' @param penalized whether to use the penalized log-likelihood
#' @param nsamp Number of approximate Gaussian sample to use in computing DIC.
#' Defaults to \code{nsamp=1e3}. Only used when the object has class 'evmOpt'.
#' @param ... other arguments currently ignored
#' @param k numeric, the penalty per parameter to be used; the
#' default \code{k = 2} is the classical AIC.
#' @param DIC Logical. Whether to compute DIC. Defaults to \code{DIC = TRUE}.
#' Only applicable to objects of class 'evmSim'.
#' @param WAIC Logical. Whether to compute WAIC. Defaults to \code{WAIC = TRUE}.
#' Only applicable to objects of class 'evmSim'.
#' @details If the object has class 'evmOpt', \code{nsamp} random draws are
#' made from the Gaussian distribution with mean and covariance inferred from
#' the model object. The result will be an approximate DIC. Note that AIC should
#' not be trusted if priors are not flat. For example, if you use a regularizing
#' prior on xi, say xi ~ N(0, 0.25), AIC can be misleading and DIC should be
#' preferred. If the object has class 'evmSim', the actual posterior draws are
#' used in the computation. Also note that sometimes the optimizer returns
#' an approximatae covariance that is not postive-semidefinite, in which case
#' the DIC will be reported as NA.
#' @return The AIC and DIC
#' @seealso \code{\link[stats]{AIC}}
#' @importFrom stats AIC logLik
#' @export
AIC.evmOpt <- function(object, penalized=FALSE, nsamp=1e3, DIC, WAIC, ..., k=2){
ll <- unclass(logLik(object, penalized=penalized))
AIC <- -2*ll + k * attr(ll, 'df')
attr(AIC, "df") <- NULL
AIC
}
DIC.evmSim <- function(object){
ll <- object$map$family$log.lik(object$map$data, object$map$th)
samp <- object$param
dev <- -2 * apply(samp, 1, ll)
#dev <- dev[dev < Inf]
mean(dev) + (mean(dev) + 2 * ll(colMeans(samp)))
}
WAIC.evmSim <- function(object){
# This follows Richard McElreath's Statistical Rethinking, Section 6.4.2
# Get density function
dfun <- object$map$family$density
samp <- object$param
# Need matrix with one column for each parameter
npar <- length(coef(object))
D <- texmexMakeNewdataD(object$map, newdata=NULL)
n <- length(object$map$data$y)
# Get loglik for every row in samp, every value of (phi, xi)
param <- texmexGetParam(data = D, co = samp)
lps <- lapply(1:nrow(D[[1]]), # For each observation get matrix of parameters
function(i, x, p){
wh <- lapply(1:length(D),
function(j, x, p, i){
rowSums(t(t(p[[j]]) * c(x[[j]][i, ])))
}, x=x, p=p, i=i)
wh <- do.call("cbind", wh)
colnames(wh) <- names(x)
wh
}, x=D, p=param)
# lps is a list of matrices, one for each observation in the original dataset,
# each matrix having columns for the parameters (phi and xi for GPD), and
# with number of rows equal to the number of retained samples from the
# Markov chains
ll <- sapply(1:length(lps), function(X){
dfun(object$map$data$y[X], lps[[X]], object$map, log.d=TRUE)
})
lppd <- apply(ll, 2, function(X){
xm <- max(X)
s <- sum(exp(X - xm))
xm + log(s)
}) - log(n)
# Get effective number of parameters
ep <- apply(ll, 2, var)
-2 * (sum(lppd) - sum(ep))
}
#' @export
AIC.evmSim <- function(object, DIC = TRUE, WAIC = TRUE, ..., k=2){
aic <- AIC.evmOpt(object$map, ..., k=2)
aic <- c(AIC = aic)
if (DIC){
dic <- DIC.evmSim(object)
aic <- c(aic, DIC = dic)
}
if (WAIC){
waic <- WAIC.evmSim(object)
aic <- c(aic, WAIC = waic)
}
aic
}
#' Log-likelihood for evmOpt objects
#'
#' Return the log-likelihood or penalized log-likelihood for
#' \code{evmOpt} objects.
#'
#' @param object fit model object
#' @param penalized whether to return the penalized log-likelihood
#' @param ... some methods need more arguments
#' @return an object of class \code{logLik}
#' @seealso \code{\link[stats]{logLik}}
#' @importFrom stats logLik
#' @export
logLik.evmOpt <- function(object, penalized=FALSE, ...){
ll <- if (penalized) {object$ploglik} else {object$loglik}
attr(ll, 'df') <- length(coef(object))
class(ll) <- 'logLik'
ll
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.