ithresh: Threshold selection in the i.i.d. case (peaks over threshold)

Description Usage Arguments Details Value References See Also Examples

View source: R/ithresh.R

Description

Produces a diagnostic plot to assist in the selection of an extreme value threshold in the case where the data can be treated as independent and identically distributed (i.i.d.) observations. For example, it could be that these observations are the cluster maxima resulting from the declustering of time series data. The predictive ability of models fitted using each of a user-supplied set of thresholds is assessed using leave-one-out cross-validation in a Bayesian setup. These models are based on a Generalized Pareto (GP) distribution for threshold excesses and a binomial model for the probability of threshold exceedance. See Northrop et al. (2017) for details.

Usage

1
ithresh(data, u_vec, ..., n_v = 1, npy = NULL, use_rcpp = TRUE)

Arguments

data

A numeric vector of observations. Any missing values will be removed. The argument npy (see below) may be supplied as an attribute of data using attr(data, "npy") <- value, where value is the value of npy (see attr). If npy is supplied twice, as both attr(data, "npy")) and using the npy argument, then the former is used.

u_vec

A numeric vector. A vector of training thresholds at which inferences are made from a binomial-GP model. These could be set at sample quantiles of data using quantile. Any duplicated values will be removed.

...

Further (optional) arguments to be passed to the revdbayes function rpost_rcpp (or rpost), which use the generalized ratio-of-uniforms method to simulate from extreme value posterior distributions. In particular:

  • n The size of the posterior sample used to perform predictive inference. Default: n = 1000.

  • prior A prior for GP parameters to be passed to the revdbayes function set_prior. Can be either a character scalar that chooses an in-built prior, or a user-supplied R function or pointer to a compiled C++ function. See the set_prior documentation for details of the in-built priors. See the revdbayes vignette Faster simulation using revdbayes for information about creating a pointer to a C++ function. See also the Examples section.

    If the user supplies and R function then rpost will be used for posterior simulation, rather than (the faster) rpost_rcpp, regardless of the input value of use_rcpp.

    Default: prior = "mdi" with a = 0.6 and min_xi = -1. This particular prior is studied in Northrop et al. (2017).

  • h_prior A list of further arguments (hyperparameters) for the GP prior specified in prior.

  • bin_prior A character scalar that chooses an in-built prior for the threshold exceedance probability p, to be passed to the revdbayes function set_bin_prior. only relevant if prior == "beta".

    Default: prior = "jeffreys", i.e. Beta(1/2, 1/2).

  • h_bin_prior A list of further arguments (hyperparameters) for the binomial prior specified in bin_prior. See the set_bin_prior documentation for details of the in-built priors.

  • trans A character scalar: either "none" or "BC". See rpost_rcpp for details. The default is "none", which is usually faster than "BC". However, if there are very few threshold excesses then using trans = "BC" can make the optimizations involved in the generalized ratio-of-uniforms algorithm more stable. If using trans = "none" produces an error for a particular posterior simulation then trans = "BC" is used instead.

n_v

A numeric scalar. Each of the n_v largest values in u_vec will be used (separately) as a validation threshold for the training thresholds in u_vec that lie at or below that validation threshold. If n_v = 1 then all the training thresholds are used with validation performed using the threshold u_vec[length(u_vec)]. If n_v = 2 then, in addition, the assessment is performed using u_vec[1], ..., u_vec[length(u_vec) - 1] with validation threshold u_vec[length(u_vec) - 1], and so on.

npy

A numeric scalar. The mean number of observations per year of data, after excluding any missing values, i.e. the number of non-missing observations divided by total number of years of non-missing data. May be supplied using as an attribute attr(data, "npy") of data instead.

The value of npy does not affect any calculation in ithresh, it only affects subsequent extreme value inferences using predict.ithresh. However, setting npy in the call to rpost, or as an attribute of data avoids the need to supply npy when calling predict.ithresh.

use_rcpp

A logical scalar. If TRUE (the default) the revdbayes function rpost_rcpp is used for posterior simulation. If FALSE, or if the user supplies an R function to set the prior for GP parameters, the (slower) function rpost is used.

Details

For a given threshold in u_vec:

See Northrop et al. (2017) and the introductory threshr vignette for further details and examples.

Value

An object (list) of class "ithresh", containing the components

References

Northrop, P.J. and Attalides, N. (2016) Posterior propriety in Bayesian extreme value analyses using reference priors Statistica Sinica, 26(2), 721–743 http://dx.doi.org/10.5705/ss.2014.034.

Northrop, P. J., Attalides, N. and Jonathan, P. (2017) Cross-validatory extreme value threshold selection and uncertainty with application to ocean storm severity. Journal of the Royal Statistical Society Series C: Applied Statistics, 66(1), 93-120. http://dx.doi.org/10.1111/rssc.12159

Jonathan, P. and Ewans, K. (2013) Statistical modelling of extreme ocean environments for marine design : a review. Ocean Engineering, 62, 91-109. http://dx.doi.org/10.1016/j.oceaneng.2013.01.004

See Also

plot.ithresh for the S3 plot method for objects of class ithresh.

summary.ithresh Summarizing measures of threshold predictive performance.

predict.ithresh for predictive inference for the largest value observed in N years.

rpost in the revdbayes package for details of the arguments that can be passed to rpost_rcpp/rpost.

set_prior and set_bin_prior in the revdbayes package for details of how to set a prior distributions for GP parameters and for the exceedance probability p.

quantile.

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
# Note:
# 1. Smoother plots result from making n larger than the default n = 1000.
# 2. In the examples below validation thresholds rather higher than is
#    advisable have been used, with far fewer excesses than the minimum of
#    50 suggested by Jonathan and Ewans (2013).

## North Sea significant wave heights, default prior -----------------------
#' # A plot akin to the top left of Figure 7 in Northrop et al. (2017)
#' # ... but with fewer training thresholds

u_vec_ns <- quantile(ns, probs = seq(0.05, 0.95, by = 0.1))
ns_cv <- ithresh(data = ns, u_vec = u_vec_ns, n_v = 3)
plot(ns_cv, lwd = 2, add_legend = TRUE, legend_pos = "topright")
mtext("significant wave height / m", side = 3, line = 2.5)

## Gulf of Mexico significant wave heights, default prior ------------------

u_vec_gom <- quantile(gom, probs = seq(0.25, 0.95, by = 0.1))
# Setting a prior using its name and parameter value(s) --------------------
# This example gives the same prior as the default
gom_cv <- ithresh(data = gom, u_vec = u_vec_gom, n_v = 4, prior = "mdi",
                  h_prior = list(a = 0.6))

## Setting a user-defined (log-)prior R function ---------------------------
# This example also gives the same prior as the default
# (It will take longer to run than the example above because ithresh detects
#  that the prior is an R function and sets use_rcpp to FALSE.)

user_prior <- function(pars, a, min_xi = -1) {
  if (pars[1] <= 0 | pars[2] < min_xi) {
    return(-Inf)
  }
  return(-log(pars[1]) - a * pars[2])
}
gom_cv <- ithresh(data = gom, u_vec = u_vec_gom, n_v = 4, prior = user_prior,
                  h_prior = list(a = 0.6))

## Setting a user-defined (log-)prior (pointer to a) C++ function ----------
# We make use of a C++ function and function create_prior_xptr() to create
# the required pointer from the revdbayes package

prior_ptr <- revdbayes:::create_prior_xptr("gp_flat")
gom_cv <- ithresh(data = gom, u_vec = u_vec_gom, n_v = 4, prior = prior_ptr,
                  h_prior = list(min_xi = -1))

threshr documentation built on Sept. 4, 2017, 9:03 a.m.