R/maq.R

In martindurocher/floodRFA: A R-package for regional frequency of floods

################################################################################
#' Flood quantile estimation by model averaging
#'
#' Estimate the flood quantile of a given period by averaging the predictions of
#' severals distributions. The weights are computed using either the BIC
#' or the AIC criterion.
#'
#' @param x Sample
#' @param distr List of distributions. See \link{vec2par}.
#' @param rt Return period in years
#' @param bic Logical. Use the BIC criterion for weight (TRUE)
#'            or the AIC (FALSE).
#'
#' @param method String determining the method of estimation.
#'               Either L-moments ('lmom') or
#'               Maximum likelihood ('mle').
#'
#' @export
#'
#' @examples
#'
#' xd <- exp(rnorm(50)+2)
#' maq(xd, method = 'lmom')
#'
#' xd <- exp(rlogis(50)+1)
#' maq(xd, distr = c('gev','glo'), bic = TRUE)

maq <- function(x, distr = c('gev','glo','pe3','ln3'),
         rt = c(2,5,10,20,50,100),
         bic = FALSE, method = 'mle'){

  ans <- list()
  ans$bic <- bic
  ans$distr <- distr
  ans$method <- method

  if(bic) pn <- log(length(x))/2
  else pn <- 1

  ## compute parameters
  if(method == 'lmom'){
    ll <- lmomco::lmoms(x)
    f <- function(z) lmomco::lmom2par(ll,z)
    pp <- lapply(distr, f)

  } else if(method == 'mle'){
    f <- function(z) try(ans <- lmomco::mle2par(x,z) [c('type','para','source')])
    pp <- lapply(distr, f)
  }

  ## compute weights
  wf <- function(z){
    ans <- NA
    try(ans <- sum(log(lmomco::dlmomco(x,z))) - pn * length(z$para),TRUE)
    return(ans)
  }
  ww <- sapply(pp, wf)
  fid <- is.finite(ww)
  ww[!fid]  <- 0 # set to zero non finite
  ww[fid] <- exp(ww[fid]-min(ww[fid]))
  ww[fid] <- ww[fid] / sum(ww[fid])
  ans$weights <- ww

  ## compute quantile
  f <- function(z){
    if(is.null(z)) return(rep(0,length(rt)))
    else return(lmomco::qlmomco(1-1/rt,z))
  }
  qq <- sapply(pp, f)

  ## Model averaging quantile
  ans$quantiles <- apply(qq, 1,function(z) sum(z*ww))
  names(ans$quantiles) <- as.character(rt)

  class(ans) <- 'maq'

  return(ans)
}

#' @export
print.maq <- function(obj){
  cat('\nModel averaging for flood quantile estimate\n\n')

  cat('Criteria : ')
  if(obj$bic) cat('BIC\n')
  else cat('AIC\n')
  cat(paste('Method : ', obj$method,'\n'))

  print(data.frame( Distr = obj$distr,
                    Weights = round(obj$weights,3)))

  cat('\nQuantiles\n')
  print(obj$quantiles)
}