#' Spectral log-likelihood function
#'
#' Compute the negative spectral log-likelihood function for Brown--Resnick model with peaks-over-threshold.
#'
#' The function compute the negative log-likelihood function based on the spectral representation developed
#' by Engelke et al. (2015). This simplified expression is obtained by conditioning on the event
#' `\code{sum(x)} exceeds a high threshold \code{u > 1}'. Margins must have been standardized.
#'
#' @param obs List of observations vectors for which \code{sum(x)} exceeds a high threshold.
#' @param loc Matrix of coordinates as given by \code{expand.grid()}.
#' @param vario Semi-variogram function taking a vector of coordinates as input.
#' @param nCores Number of cores used for the computation
#' @param cl Cluster instance as created by \code{makeCluster} of the \code{parallel} package.
#' @return Negative spectral log-likelihood function evaluated at the set of observations \code{obs} with semi-variogram \code{vario}.
#' @examples
#' #Define semi-variogram function
#' vario <- function(h){
#' 1 / 2 * norm(h,type = "2")^1.5
#' }
#'
#' #Define locations
#' loc <- expand.grid(1:4, 1:4)
#'
#' #Simulate data
#' obs <- simulPareto(1000, loc, vario)
#'
#' #Evaluate risk functional
#' sums <- sapply(obs, sum)
#'
#' #Select exceedances
#' exceedances <- obs[sums > quantile(sums, 0.9)]
#'
#' #Evaluate the spectral function
#' spectralLikelihood(exceedances, loc, vario)
#' @export
#' @references Engelke, S. et al. (2015). Estimation of Huesler-Reiss distributions and Brown-Resnick processes. Journal of the Royal Statistical Society: Series B, 77(1), 239-265
spectralLikelihood <- function(obs, loc, vario, nCores = 1L, cl = NULL){
if(is.matrix(obs)){ #Not converted to list
obs <- split(obs, row(obs)) #create list
}
if(!inherits(obs, "list") || length(obs) < 1 || !inherits(obs[[1]], c("numeric","integer"))){
stop('obs must be a list of vectors')
}
if(!inherits(loc, c("matrix", "data.frame"))) {
stop('`loc` must be a data frame of coordinates as generated by `expand.grid()` or a matrix of locations (one site per row)')
}
if(!is.numeric(nCores) || nCores < 1) {
stop('`nCores` must a positive number of cores to use for parallel computing.')
}
if(nCores > 1 && !inherits(cl, "cluster")) {
stop('For parallel computation, `cl` must an cluster created by `makeCluster` of the package `parallel`.')
}
n <- length(obs)
dim <- nrow(loc)
if(dim != length(obs[[1]])){
stop('The size of the vectors of observations does not match the grid size.')
}
gamma <- tryCatch({
dists <- lapply(1:ncol(loc), function(i) {
outer(loc[,i],loc[,i], "-")
})
computeVarMat <- sapply(1:length(dists[[1]]), function(i){
h <- rep(0,ncol(loc))
for(j in 1:ncol(loc)){
h[j] = dists[[j]][i]
}
vario(h)
})
matrix(computeVarMat, dim, dim)
}, warning = function(war) {
war
}, error = function(err) {
stop('The semi-variogram is not valid for the locations provided.')
})
psi <- (outer(gamma[-1,1], gamma[-1,1], "+") - (gamma[-1,-1]))
invPsi <- MASS::ginv(psi)
productCovMatrix = function(i){
omega <- log(obs[[i]][-1]/obs[[i]][1]) + gamma[-1,1]
res <- t(omega) %*% invPsi %*% omega
}
if(nCores > 1){
likelihood <- parallel::parLapply(cl,1:n, productCovMatrix)
} else {
likelihood <- lapply(1:n, productCovMatrix)
}
logdetA = determinant(psi, logarithm = TRUE)$modulus
(n/2 * logdetA + 1/2 * sum(unlist(likelihood)))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.