Add maxstable

Author: lbelzile

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: mev
 Type: Package
 Title: Modelling of Extreme Values
-Version: 1.17.10006
+Version: 1.17.10007
 Authors@R: c(person(given="Leo", family="Belzile", role = c("aut", "cre"), email = "belzilel@gmail.com", comment = c(ORCID = "0000-0002-9135-014X")), person(given="Jennifer L.", family="Wadsworth", role=c("aut")), person(given="Paul J.", family="Northrop", role=c("aut"), comment = c(ORCID = "0000-0002-1992-4882")), person(given="Scott D.", family="Grimshaw", role=c("aut"), comment = c(ORCID = "0000-0002-6326-9360")), person(given="Jin", family="Zhang", role=c("ctb")), person(given="Michael A.", family="Stephens", role=c("ctb")), person(given="Art B.", family="Owen", role=c("ctb")), person(given="Raphael", family="Huser", role=c("aut"), comment = c(ORCID = "0000-0002-1228-2071")))
 Description: Various tools for the analysis of univariate, multivariate and functional extremes. Exact simulation from max-stable processes [Dombry, Engelke and Oesting (2016) <doi:10.1093/biomet/asw008>, R-Pareto processes for various parametric models, including Brown-Resnick (Wadsworth and Tawn, 2014, <doi:10.1093/biomet/ast042>) and Extremal Student (Thibaud and Opitz, 2015, <doi:10.1093/biomet/asv045>). Threshold selection methods, including Wadsworth (2016) <doi:10.1080/00401706.2014.998345>, and Northrop and Coleman (2014) <doi:10.1007/s10687-014-0183-z>. Multivariate extreme diagnostics. Estimation and likelihoods for univariate extremes, e.g., Coles (2001) <doi:10.1007/978-1-4471-3675-0>.
 License: GPL-3

NAMESPACE

@@ -198,6 +198,7 @@ export(likmgp)
 export(maxstabtest)
 export(mvrnorm)
 export(pegp)
+export(penultimate)
 export(pextgp)
 export(pextgp.G)
 export(pgev)

R/asymcoef.R

@@ -5,12 +5,12 @@
 #'  Semadeni's (2021) PhD thesis.
 #'  Two estimators are implemented: one based on empirical distributions, the second using empirical likelihood.
 #' @details Let \code{U}, \code{V} be uniform random variables and define the partial extremal dependence coefficients
-#' \eqn{\varphi_{+}(u) = \Pr(V > U | U > u, V > u)},
-#' \eqn{\varphi_{-}(u) = \Pr(V < U | U > u, V > u)} and
-#' \eqn{\varphi_0(u) = \Pr(V = U | U > u, V > u)}
+#' \deqn{\varphi_{+}(u) = \Pr(V > U | U > u, V > u),},
+#' \deqn{\varphi_{-}(u) = \Pr(V < U | U > u, V > u),}
+#' \deqn{\varphi_0(u) = \Pr(V = U | U > u, V > u).}
 #' Define
 #' \deqn{ \varphi(u) = \frac{\varphi_{+} - \varphi_{-}}{\varphi_{+} + \varphi_{-}}}
-# and the coefficient of extremal asymmetry as \eqn{\varphi = \lim_{u \to 1} \varphi(u)}.
+#' and the coefficient of extremal asymmetry as \eqn{\varphi = \lim_{u \to 1} \varphi(u)}.
 #'
 #' The empirical likelihood estimator, derived for max-stable vectors with unit Frechet margins, is
 #' \deqn{\frac{\sum_i p_i I(w_i \leq 0.5) - 0.5}{0.5 - 2\sum_i p_i(0.5-w_i) I(w_i \leq 0.5)}}

R/extgp.R

@@ -755,7 +755,24 @@ tstab.egp <- function(
       stop("Threshold vector not provided")
     }
   }
-  if (model %in% c("egp1", "egp2", "egp3") & length(model) == 1) {
+  model <- match.arg(
+    model,
+    choices = c(
+      "pt-beta",
+      "pt-gamma",
+      "pt-power",
+      "gj-tnorm",
+      "gj-beta",
+      "exptilt",
+      "logist",
+      "egp1",
+      "egp2",
+      "egp3"
+    ),
+    several.ok = FALSE
+  )
+
+  if (model %in% c("egp1", "egp2", "egp3")) {
     model <- switch(
       model,
       egp1 = "pt-beta",

R/kernel_exponential_thselect.R

@@ -10,7 +10,7 @@
 #' \item{\code{k0}:} number of exceedances
 #' \item{\code{shape}:} Hill's shape estimate
 #' \item{\code{rho}:} second-order regular variation parameter estimate
-#' \item{\code{gof}:} goodness-of-fit statistic for the \"best\" threshold.
+#' \item{\code{gof}:} goodness-of-fit statistic for the chosen threshold.
 #' }
 #' @references Goegebeur , Y., Beirlant , J., and de Wet , T. (2008). Linking Pareto-Tail Kernel Goodness-of-fit Statistics with Tail Index at Optimal Threshold and Second Order Estimation. REVSTAT-Statistical Journal, 6(\bold{1}), 51–69. <doi:10.57805/revstat.v6i1.57>
 #' @export

R/lthill.R

@@ -243,7 +243,7 @@ tstab.lthill <- function(
 
 #' Lower truncated Hill threshold selection
 #'
-#' Given a sample of positive data with Pareto tail, the algorithm computes the optimal number of order statistics that minimizes the variance of the average left truncated tail index estimator, and uses the relationship to the Hill estimator for the Hall class of distributions to derive the optimal number (minimizing the asymptotic mean squared error) of the Hill estimator. The default value for the second order regular variation index is taken to be \eqn{rho=-1}.
+#' Given a sample of positive data with Pareto tail, the algorithm computes the optimal number of order statistics that minimizes the variance of the average left truncated tail index estimator, and uses the relationship to the Hill estimator for the Hall class of distributions to derive the optimal number (minimizing the asymptotic mean squared error) of the Hill estimator. The default value for the second order regular variation index is taken to be \eqn{\rho=-1}.
 #' @param xdat [vector] positive vector of exceedances
 #' @param kmin [int] minimum number of exceedances
 #' @param kmax [int] maximum number of exceedances for the estimation of the shape parameter.

R/maxstabtest.R

@@ -17,7 +17,7 @@
 #' @param ties.method string indicating the method for \code{\link[base]{rank}}. Default to \code{"random"}.
 #' @param ... additional arguments for backward compatibility
 #' @param plot logical indicating whether a graph should be produced (default to \code{TRUE}).
-#' @return a Tukey probability-probability plot with 95% confidence intervals obtained using a nonparametric bootstrap
+#' @return a Tukey probability-probability plot with 95\% confidence intervals obtained using a nonparametric bootstrap
 #' @export
 #' @examples
 #' \dontrun{

R/mle.R

@@ -1003,7 +1003,7 @@ fit.rlarg <- function(
 #' @export
 plot.mev_gpd <- function(
   x,
-  which = 1:2,
+  which = c("pp", "qq"),
   main,
   xlab = "Theoretical quantiles",
   ylab = "Sample quantiles",
@@ -1015,12 +1015,17 @@ plot.mev_gpd <- function(
       "Object \"x\" does not contain \"exceedances\", or else the latter is not a vector"
     )
   }
-  if (!is.numeric(which) || any(which < 1) || any(which > 2)) {
-    stop("\"which\" must be in 1:2")
-  }
   show <- rep(FALSE, 2)
-  show[which] <- TRUE
-  if (!add) {
+  if (!is.numeric(which)) {
+    which <- match.arg(which)
+    show <- c("pp", "qq") %in% which
+  } else {
+    if (!isTRUE(all(which %in% 1:2))) {
+      stop("\"which\" must be in 1:2 if not a string")
+    }
+    show[which] <- TRUE
+  }
+  if (!isTRUE(add)) {
     old.par <- par(no.readonly = TRUE)
     if (sum(show) == 2) {
       par(mfrow = c(1, 2), mar = c(5, 5, 4, 1))
@@ -1113,7 +1118,7 @@ plot.mev_gpd <- function(
 #' @export
 plot.mev_gpdbayes <- function(
   x,
-  which = 1:2,
+  which = c("pp", "qq"),
   main,
   xlab = "Theoretical quantiles",
   ylab = "Sample quantiles",
@@ -1125,12 +1130,17 @@ plot.mev_gpdbayes <- function(
       "Object \"x\" does not contain \"exceedances\", or else the latter is not a vector"
     )
   }
-  if (!is.numeric(which) || any(which < 1) || any(which > 2)) {
-    stop("\"which\" must be in 1:2")
-  }
   show <- rep(FALSE, 2)
-  show[which] <- TRUE
-  if (!add) {
+  if (!is.numeric(which)) {
+    which <- match.arg(which)
+    show <- c("pp", "qq") %in% which
+  } else {
+    if (!isTRUE(all(which %in% 1:2))) {
+      stop("\"which\" must be in 1:2 if not a string")
+    }
+    show[which] <- TRUE
+  }
+  if (!isTRUE(add)) {
     old.par <- par(no.readonly = TRUE)
     if (sum(show) == 2) {
       par(mfrow = c(1, 2), mar = c(5, 5, 4, 1))
@@ -1229,27 +1239,35 @@ plot.mev_gpdbayes <- function(
 #' @export
 plot.mev_gev <- function(
   x,
-  which = 1:2,
+  which = c("pp", "qq"),
   main,
   xlab = "Theoretical quantiles",
   ylab = "Sample quantiles",
+  add = TRUE,
   ...
 ) {
   if (!is.vector(x$xdat)) {
     stop(
       "Object \"x\" does not contain \"exceedances\", or else the latter is not a vector"
     )
   }
-  if (!is.numeric(which) || any(which < 1) || any(which > 2)) {
-    stop("\"which\" must be in 1:2")
-  }
   show <- rep(FALSE, 2)
-  show[which] <- TRUE
-  old.par <- par(no.readonly = TRUE)
-  if (sum(show) == 2) {
-    par(mfrow = c(1, 2), mar = c(5, 5, 4, 1))
+  if (!is.numeric(which)) {
+    which <- match.arg(which)
+    show <- c("pp", "qq") %in% which
+  } else {
+    if (!isTRUE(all(which %in% 1:2))) {
+      stop("\"which\" must be in 1:2 if not a string")
+    }
+    show[which] <- TRUE
+  }
+  if (!isTRUE(add)) {
+    old.par <- par(no.readonly = TRUE)
+    if (sum(show) == 2) {
+      par(mfrow = c(1, 2), mar = c(5, 5, 4, 1))
+    }
+    on.exit(par(old.par))
   }
-  on.exit(par(old.par))
   if (missing(main)) {
     main <- c("Probability-probability plot", "Quantile-quantile plot")
   } else {
@@ -1329,10 +1347,11 @@ plot.mev_gev <- function(
 #' @export
 plot.mev_rlarg <- function(
   x,
-  which = 1:2,
+  which = c("pp", "qq"),
   main,
   xlab = "Theoretical quantiles",
   ylab = "Sample quantiles",
+  add = TRUE,
   ...
 ) {
   if (!isTRUE(all.equal(x$param[3], 0, check.attributes = FALSE))) {
@@ -1348,16 +1367,24 @@ plot.mev_rlarg <- function(
     dat <- sort(as.vector(c(ppdat[, 1], apply(ppdat, 1, diff))))
   }
   n <- length(dat)
-  if (!is.numeric(which) || any(which < 1) || any(which > 2)) {
-    stop("\"which\" must be in 1:2")
-  }
   show <- rep(FALSE, 2)
+  if (!is.numeric(which)) {
+    which <- match.arg(which)
+    show <- c("pp", "qq") %in% which
+  } else {
+    if (!isTRUE(all(which %in% 1:2))) {
+      stop("\"which\" must be in 1:2 if not a string")
+    }
+    show[which] <- TRUE
+  }
   show[which] <- TRUE
-  old.par <- par(no.readonly = TRUE)
-  if (sum(show) == 2) {
-    par(mfrow = c(1, 2), mar = c(5, 5, 4, 1))
+  if (!isTRUE(add)) {
+    old.par <- par(no.readonly = TRUE)
+    if (sum(show) == 2) {
+      par(mfrow = c(1, 2), mar = c(5, 5, 4, 1))
+    }
+    on.exit(par(old.par))
   }
-  on.exit(par(old.par))
   if (missing(main)) {
     main <- c("Probability-probability plot", "Quantile-quantile plot")
   } else {
@@ -1432,27 +1459,35 @@ plot.mev_rlarg <- function(
 #' @export
 plot.mev_pp <- function(
   x,
-  which = 1:2,
+  which = c("pp", "qq"),
   main = "Quantile-quantile plot",
   xlab = "Theoretical quantiles",
   ylab = "Sample quantiles",
+  add = TRUE,
   ...
 ) {
   if (!is.vector(x$exceedances)) {
     stop(
       "Object \"x\" does not contain \"exceedances\", or else the latter is not a vector"
     )
   }
-  if (!is.numeric(which) || any(which < 1) || any(which > 2)) {
-    stop("\"which\" must be in 1:2")
-  }
   show <- rep(FALSE, 2)
-  show[which] <- TRUE
-  old.par <- par(no.readonly = TRUE)
-  if (sum(show) == 2) {
-    par(mfrow = c(1, 2), mar = c(5, 5, 4, 1))
+  if (!is.numeric(which)) {
+    which <- match.arg(which)
+    show <- c("pp", "qq") %in% which
+  } else {
+    if (!isTRUE(all(which %in% 1:2))) {
+      stop("\"which\" must be in 1:2 if not a string")
+    }
+    show[which] <- TRUE
+  }
+  if (!isTRUE(add)) {
+    old.par <- par(no.readonly = TRUE)
+    if (sum(show) == 2) {
+      par(mfrow = c(1, 2), mar = c(5, 5, 4, 1))
+    }
+    on.exit(par(old.par))
   }
-  on.exit(par(old.par))
   if (missing(main)) {
     main <- c("Probability-probability plot", "Quantile-quantile plot")
   } else {

R/mrl.R

@@ -204,7 +204,7 @@ automrl <- function(
 
 #' Mean residual life plot
 #'
-#' Computes mean of sample exceedances over a range of thresholds or for a prespecified number of largest order statistics, and returns a plot with 95% Wald-based confidence intervals as a function of either the threshold or the number of exceedances. The main purpose is the plotting method, which generates the so-called mean residual life plot. The latter should be approximately linear over the threshold for a generalized Pareto distribution
+#' Computes mean of sample exceedances over a range of thresholds or for a pre-specified number of largest order statistics, and returns a plot with 95\% Wald-based confidence intervals as a function of either the threshold or the number of exceedances. The main purpose is the plotting method, which generates the so-called mean residual life plot. The latter should be approximately linear over the threshold for a generalized Pareto distribution
 #'
 #' @references Davison, A.C. and R.L. Smith (1990). Models for Exceedances over High Thresholds (with discussion), \emph{Journal of the Royal Statistical Society. Series B (Methodological)}, \bold{52}(3), 393--442.
 #' @export