floodStat: FloodNet Regional Frequency Analysis

###########################################################################
#' Generalized Pareto distribution (GPA)
#'
#' Distribution, density, quantile and random function for the Generalized
#' pareto distriution.
#'
#' @name GPA
#' @param p,q Probabilities or quantiles for the GPA distribution.
#'
#' @param n Number of simulations.
#'
#' @param alpha Scale parameter of the GPA
#'
#' @param kap Shape parameter of the GPA
#'
#' @param lower.tail Should the propability of the lower tail be returned
#'
#' @param log Should the log-density be returned
#'
#' @details
#'
#' These function are reorganisation of some functions in R-package evd
#' using the parametrisation (\eqn{\alpha}, \eqn{\kappa}) that is coherent
#' with the notation frequently used in hydrology.
#'
#' @export
#'
#' @examples
#'
#' kap <- -.2
#' a <- 1
#' xd1 <- rgpa(1e4, a, kap)
#' xd2 <- qgpa(runif(1e4), a, kap)
#'
#' qqplot(xd1, xd2)
#'
#'
#' tt <- seq(0.001,6, len = 100)
#'
#' hist(xd1[xd1<6], main = 'GPA distribution',
#'   freq = FALSE, ylim =c(0,1), xlim = c(0,6))
#'
#' lines(tt, dgpa(tt,a,kap))
#' lines(tt,pgpa(sort(tt), a, kap), col = 2, lty = 2)
#'
#'
pgpa <- function (q,  alpha = 1, kap = 0, lower.tail = TRUE)
{
  if (min(alpha) <= 0)
    stop("invalid alpha")

  if (length(kap) != 1)
    stop("invalid kappa")

  q <- pmax(q, 0)/alpha

  if (kap == 0)
    p <- 1 - exp(-q)
  else {
    p <- pmax(1 - kap * q, 0)
    p <- 1 - p^(1/kap)
  }

  if (!lower.tail)
    p <- 1 - p

  return(p)
}

#' @export
#' @rdname GPA
rgpa <- function (n, alpha = 1, kap = 0)
{
  if (min(alpha) < 0)
    stop("invalid alpha")
  if (length(kap) != 1)
    stop("invalid kappa")
  if (kap == 0)
    return(alpha * rexp(n))
  else return(- alpha * (runif(n)^(kap) - 1)/kap)
}

#' @export
#' @rdname GPA
dgpa <- function (x, alpha = 1, kap = 0, log = FALSE)
{
  if (min(alpha) <= 0)
    stop("invalid alpha")

  if (length(kap) != 1)
    stop("invalid kappa")

  d <- x/alpha
  nn <- length(d)
  alpha <- rep(alpha, length.out = nn)
  index <- (d > 0 & ((1 - kap * d) > 0)) | is.na(d)

  if (kap == 0) {
    d[index] <- log(1/alpha[index]) - d[index]
    d[!index] <- -Inf
  }
  else {
    d[index] <- log(1/alpha[index]) -
                            (1 - 1/kap) * log(1 - kap * d[index])
    d[!index] <- -Inf
  }

  if (!log)
    d <- exp(d)

  return(d)
}

#' @export
#' @rdname GPA
qgpa <- function (p, alpha = 1, kap = 0, lower.tail = TRUE)
{
  if (min(p, na.rm = TRUE) < 0 || max(p, na.rm = TRUE) > 1)
    stop("`p' must contain probabilities in (0,1)")

  if (min(alpha) < 0)
    stop("invalid alpha")

  if (length(kap) != 1)
    stop("invalid kappa")

  if (lower.tail)
    p <- 1 - p

  if (kap == 0)
    return(- alpha * log(p))
  else
    return(-alpha * (p^(kap) - 1)/kap)
}

###############################################################
#' Estimation of the Generalized pareto distribution.
#'
#' Low level functions for the estimation of the generalized pareto
#' distribution(GPA) with two parameters. Can use
#' either maximum likelihoood or the method of L-moments.
#' The algorithm of \code{fgpa2d} is using
#' \code{optim} to directly optimize the log-likelihood (bivariate), while
#' the algorithm of \code{fgpa1d} is using a transformation to use
#' a univariate optimization routine. Moreover, \code{fgpa2d} constraint
#' the shape parameter between -.5 and 1.
#'
#' @name fgpa
#' @param x Sample.
#'
#' @param sol Does solution from \code{optim} be returned.
#'   In case of \code{fgpa1d}, it returns the variance covariance matrix.
#'
#' @param par0 Initial parameter.
#'
#' @param ... aditional arguments to pass to \code{\link{optim}}
#'
#' @section Reference:
#'
#' Davison AC, Smith RL. (1990) Models for Exceedances over High Thresholds.
#'   Journal of the Royal Statistical Society Series B (Methodological).
#'   52(3):393–442.
#'
#' Hosking JRM (1990). L-Moments: Analysis and Estimation of Distributions Using
#'   Linear Combinations of Order Statistics. Journal of the Royal Statistical
#'   Society Series B (Methodological). 52(1):105–24.
#'
#' @export
#'
#' @examples
#'
#' x <- rgpa(1000, 1, -.2)
#' fgpa1d(x)
#' fgpa2d(x)
#' fgpaLmom(x)
#'
fgpa1d <- function(x, sol = FALSE){

  fk <- function(sigma) -mean(log(1-x/sigma))

  fpara <- function(sigma){
    length(x) * (-log(fk(sigma) * sigma) + fk(sigma) - 1)
  }

  mx <- max(x)

  sigma <- NA
  suppressWarnings(try(
    sigma <- optimize(fpara, interval = c(-10000,10000) * mx,
                      maximum = TRUE)$maximum))

  para <- rep(NA,2)
  names(para) <- c('alpha','kappa')
  try(para[2] <- fk(sigma))
  try(para[1] <- para[2]*sigma)

  if(sol){
    ans <- list(par = para, varcov = matrix(NA,2,2))
    ans$varcov[2,2] <- 1-para[2]
    ans$varcov[1,1] <- 2 * para[1]^2
    ans$varcov[1,2] <- ans$varcov[2,1] <- para[1]
    ans$varcov <- ans$varcov * ans$varcov[2,2] / length(x)

  } else {
    ans <- para
  }

  return(ans)

}

#' @export
#' @rdname fgpa
fgpa2d <- function(x, sol = FALSE, par0 = NULL, ...){

  logit0 <- function(z) logit((z+.5)/1.5)
  expit0 <- function(z) expit(z)*1.5 -.5

  ## negative loglikelihood
  nllik <- function(para) -sum(dgpa(x,exp(para[1]),
                                    expit0(para[2]), log = TRUE))

  if(!is.null(par0)){
    par0[1] <- log(par0[1])
    par0[2] <- logit0(par0[2])
  } else
    par0 <- c(log(mean(x)),-.6931472)

  ## estimate the parameter
  para <- optim( par0, nllik, ...)$par
  para[1] <- exp(para[1])
  para[2] <- expit0(para[2])

  names(para) <- c('alpha','kappa')

  if(sol){
    ans <- list(par = para, varcov = matrix(NA,2,2))
    ans$varcov[2,2] <- 1-para[2]
    ans$varcov[1,1] <- 2 * para[1]^2
    ans$varcov[1,2] <- ans$varcov[2,1] <- para[1]
    ans$varcov <- ans$varcov * ans$varcov[2,2] / length(x)

  } else {
    ans <- para
  }


  return(ans)
}

#' @export
#' @rdname fgpa
fgpaLmom <- function(x){
  lmom0 <- lmomco::lmoms(x, nmom = 3)
  return(lmomco::lmom2par(lmom0,'gpa', xi = 0)$para[-1])
}