revdbayes: Ratio-of-Uniforms Sampling for Bayesian Extreme Value Analysis

Description

Uses the multivariate generalized ratio-of-uniforms method to simulate random samples from the posterior distributions commonly encountered in Bayesian extreme value analyses.

Details

The main functions in the revdbayes package are rpost and rpost_rcpp, which simulate random samples from the posterior distribution of extreme value model parameters using the functions ru and ru_rcpp from the rust package, respectively. The user chooses the extreme value model, the prior density for the parameters and provides the data. There are options to improve the probability of acceptance of the ratio-of-uniforms algorithm by working with transformation of the model parameters.

The functions kgaps_post and dgaps_post simulate from the posterior distribution of the extremal index \theta based on the K-gaps model for threshold interexceedance times of Suveges and Davison (2010) and the similar D-gaps model of Holesovsky and Fusek (2020). See also Attalides (2015).

See vignette("revdbayes-a-vignette", package = "revdbayes") for an overview of the package and vignette("revdbayes-b-using-rcpp-vignette", package = "revdbayes") for an illustration of the improvements in efficiency produced using the Rcpp package. See vignette("revdbayes-c-predictive-vignette", package = "revdbayes") for an outline of how to use revdbayes to perform posterior predictive extreme value inference and vignette("revdbayes-d-kgaps-vignette", package = "revdbayes") considers Bayesian inference for the extremal index \theta using threshold inter-exceedance times.

Author(s)

Maintainer: Paul J. Northrop p.northrop@ucl.ac.uk [copyright holder]

Other contributors:

• Scott D. Grimshaw [contributor]

References

Holesovsky, J. and Fusek, M. Estimation of the extremal index using censored distributions. Extremes 23, 197-213 (2020). doi:10.1007/s10687-020-00374-3

Northrop, P. J. (2016). rust: Ratio-of-Uniforms Simulation with Transformation. R package version 1.2.2. https://cran.r-project.org/package=rust.

Suveges, M. and Davison, A. C. (2010) Model misspecification in peaks over threshold analysis, The Annals of Applied Statistics, 4(1), 203-221. doi:10.1214/09-AOAS292

Attalides, N. (2015) Threshold-based extreme value modelling, PhD thesis, University College London. https://discovery.ucl.ac.uk/1471121/1/Nicolas_Attalides_Thesis.pdf

See Also

set_prior to set a prior density for extreme value parameters.

rpost and rpost_rcpp to perform ratio-of-uniforms sampling from an extreme value posterior distribution.

kgaps_post and dgaps_post to sample from a posterior distribution for the extremal index based on inter-exceedance times.

The ru and ru_rcpp functions in the rust package for details of the arguments that can be passed to ru via rpost and for the form of the object (of class "evpost") returned from rpost, which has the same structure as an object (of class "ru") returned by ru and ru_rcpp.

Random sampling from a binomial posterior distribution

Description

Samples from the posterior distribution of the probability p of a binomial distribution.

Usage

binpost(n, prior, ds_bin, param = c("logit", "p"))

Arguments

Details

If prior$prior == "bin_beta" then the posterior for p is a beta distribution so rbeta is used to sample from the posterior.

If prior$prior == "bin_mdi" then rejection sampling is used to sample from the posterior with an envelope function equal to the density of a beta(ds$m + 1, ds$n_raw - ds$m + 1) density.

If prior$prior == "bin_northrop" then rejection sampling is used to sample from the posterior with an envelope function equal to the posterior density that results from using a Haldane prior.

If prior$prior is a (user-supplied) R function then ru is used to sample from the posterior using the generalised ratio-of-uniforms method.

Value

An object (list) of class "binpost" with components

If prior$prior is a (user-supplied) R function then this list also contains ru_object the object of class "ru" returned by ru.

See Also

set_bin_prior for setting a prior distribution for the binomial probability p.

Examples

u <- quantile(gom, probs = 0.65)
ds_bin <- list()
ds_bin$n_raw <- length(gom)
ds_bin$m <- sum(gom > u)
bp <- set_bin_prior(prior = "jeffreys")
temp <- binpost(n = 1000, prior = bp, ds_bin = ds_bin)
graphics::hist(temp$bin_sim_vals, prob = TRUE)

# Setting a beta prior (Jeffreys in this case) by hand
beta_prior_fn <- function(p, ab) {
  return(stats::dbeta(p, shape1 = ab[1], shape2 = ab[2], log = TRUE))
}
jeffreys <- set_bin_prior(beta_prior_fn, ab = c(1 / 2, 1 / 2))
temp <- binpost(n = 1000, prior = jeffreys, ds_bin = ds_bin)

Create an external pointer to a C++ prior

Description

This function provides an example of a way in which a user can specify their own prior density to rpost_rcpp. More specifically, a function like this (the user will need to create an edited version tailored to their own C++ function(s)) can be used to generate an external pointer to a compiled C++ function that evaluates the log-prior density. Please see the vignette "Faster simulation using revdbayes" for more information.

Usage

create_prior_xptr(fstr)

Arguments

Details

Suppose that the user's C++ functions are in a file called "user_fns.cpp". These functions must be compiled and made available to R before the pointer is created. This can be achieved using the function sourceCpp in the Rcpp package or using RStudio's Source button on the editor toolbar.

For details see the examples in the documentation of the functions rpost_rcpp and set_prior, the vignette "Faster simulation using revdbayes" and the vignette "Rusting Faster: Simulation using Rcpp" in the package rust.

Value

An external pointer.

See Also

set_prior to specify a prior distribution using an external pointer returned by create_prior_xptr and for details of in-built named prior distributions.

The examples in the documentation of rpost_rcpp.

Examples

ptr_gp_flat <- create_prior_xptr("gp_flat")
prior_cfn <- set_prior(prior = ptr_gp_flat, model = "gp", min_xi = -1)

ptr_gev_flat <- create_prior_xptr("gev_flat")
prior_cfn <- set_prior(prior = ptr_gev_flat, model = "gev", min_xi = -1,
                       max_xi = Inf)

mat <- diag(c(10000, 10000, 100))
ptr_gev_norm <- create_prior_xptr("gev_norm")
pn_u <- set_prior(prior = ptr_gev_norm, model = "gev", mean = c(0,0,0),
                  icov = solve(mat))

Random sampling from D-gaps posterior distribution

Description

Uses the rust package to simulate from the posterior distribution of the extremal index \theta based on the D-gaps model for threshold interexceedance times of Holesovsky and Fusek (2020). We refer to this as the D-gaps model, because it uses a tuning parameter D, whereas the related K-gaps model of Suveges and Davison (2010) has a tuning parameter K.

Usage

dgaps_post(
  data,
  thresh,
  D = 1,
  n = 1000,
  inc_cens = TRUE,
  alpha = 1,
  beta = 1,
  param = c("logit", "theta"),
  use_rcpp = TRUE
)

Arguments

Details

A beta(\alpha, \beta) prior distribution is used for \theta so that the posterior from which values are simulated is proportional to

\theta ^ {2 N_1 + \alpha - 1} (1 - \theta e^{-\theta d}) ^ {N_0 + \beta - 1} \exp\{- \theta q (I_0 T_0 + \cdots + I_N T_N)\}.

See dgaps_stat for a description of the variables involved in the contribution of the likelihood to this expression.

The ru function in the rust package simulates from this posterior distribution using the generalised ratio-of-uniforms distribution. To improve the probability of acceptance, and to ensure that the simulation will work even in extreme cases where the posterior density of \theta is unbounded as \theta approaches 0 or 1, we simulate from the posterior distribution of \phi = \log(\theta / (1-\theta)) and then transform back to the \theta-scale.

Value

An object (list) of class "evpost", which has the same structure as an object of class "ru" returned from ru. In addition this list contains

• call: The call to dgaps().

• model: The character scalar "dgaps".

• thresh: The argument thresh.

• ss: The sufficient statistics for the D-gaps likelihood, as calculated by dgaps_stat.

References

Holesovsky, J. and Fusek, M. Estimation of the extremal index using censored distributions. Extremes 23, 197-213 (2020). doi:10.1007/s10687-020-00374-3

Suveges, M. and Davison, A. C. (2010) Model misspecification in peaks over threshold analysis, The Annals of Applied Statistics, 4(1), 203-221. doi:10.1214/09-AOAS292

See Also

ru for the form of the object returned by dgaps_post.

kgaps_post for Bayesian inference about the extremal index \theta using the K-gaps model.

Examples

# Newlyn sea surges

thresh <- quantile(newlyn, probs = 0.90)
d_postsim <- dgaps_post(newlyn, thresh)
plot(d_postsim)

### Cheeseboro wind gusts

d_postsim <- dgaps_post(exdex::cheeseboro, thresh = 45, D = 3)
plot(d_postsim)

The Generalised Extreme Value Distribution

Description

Density function, distribution function, quantile function and random generation for the generalised extreme value (GEV) distribution.

Usage

dgev(x, loc = 0, scale = 1, shape = 0, log = FALSE, m = 1)

pgev(q, loc = 0, scale = 1, shape = 0, lower.tail = TRUE, log.p = FALSE, m = 1)

qgev(p, loc = 0, scale = 1, shape = 0, lower.tail = TRUE, log.p = FALSE, m = 1)

rgev(n, loc = 0, scale = 1, shape = 0, m = 1)

Arguments

Details

The distribution function of a GEV distribution with parameters loc = \mu, scale = \sigma (> 0) and shape = \xi is

F(x) = \exp\{-[1 + \xi (x - \mu) / \sigma] ^ {-1/\xi} \}

for 1 + \xi (x - \mu) / \sigma > 0. If \xi = 0 the distribution function is defined as the limit as \xi tends to zero. The support of the distribution depends on \xi: it is x \leq \mu - \sigma / \xi for \xi < 0; x \geq \mu - \sigma / \xi for \xi > 0; and x is unbounded for \xi = 0. Note that if \xi < -1 the GEV density function becomes infinite as x approaches \mu -\sigma / \xi from below.

If lower.tail = TRUE then if p = 0 (p = 1) then the lower (upper) limit of the distribution is returned, which is -Inf or Inf in some cases. Similarly, but reversed, if lower.tail = FALSE.

See https://en.wikipedia.org/wiki/Generalized_extreme_value_distribution for further information.

The effect of m is to change the location, scale and shape parameters to (\mu + \sigma \log m, \sigma, \xi) if \xi = 0 and (\mu + \sigma (m ^ \xi - 1) / \xi, \sigma m ^ \xi, \xi). For integer m we can think of this as working with the maximum of m independent copies of the original GEV(loc, scale, shape) variable.

Value

dgev gives the density function, pgev gives the distribution function, qgev gives the quantile function, and rgev generates random deviates.

The length of the result is determined by n for rgev, and is the maximum of the lengths of the numerical arguments for the other functions.

The numerical arguments other than n are recycled to the length of the result.

References

Jenkinson, A. F. (1955) The frequency distribution of the annual maximum (or minimum) of meteorological elements. Quart. J. R. Met. Soc., 81, 158-171. doi:10.1002/qj.49708134804

Coles, S. G. (2001) An Introduction to Statistical Modeling of Extreme Values, Springer-Verlag, London. Chapter 3: doi:10.1007/978-1-4471-3675-0_3

Examples

dgev(-1:4, 1, 0.5, 0.8)
dgev(1:6, 1, 0.5, -0.2, log = TRUE)
dgev(1, shape = c(-0.2, 0.4))

pgev(-1:4, 1, 0.5, 0.8)
pgev(1:6, 1, 0.5, -0.2)
pgev(1, c(1, 2), c(1, 2), c(-0.2, 0.4))
pgev(-3, c(1, 2), c(1, 2), c(-0.2, 0.4))
pgev(7, 1, 1, c(-0.2, 0.4))

qgev((1:9)/10, 2, 0.5, 0.8)
qgev(0.5, c(1,2), c(0.5, 1), c(-0.5, 0.5))

p <- (1:9)/10
pgev(qgev(p, 1, 2, 0.8), 1, 2, 0.8)

rgev(6, 1, 0.5, 0.8)

Beta-type prior for GEV shape parameter \xi

Description

For information about this and other priors see set_prior.

Usage

gev_beta(pars, min_xi = -1/2, max_xi = 1/2, pq = c(6, 9), trendsd = 0)

Arguments

Value

The log of the prior density.

Flat prior for GEV parameters (\mu, log \sigma, \xi)

Description

For information about this and other priors see set_prior.

Usage

gev_flat(pars, min_xi = -Inf, max_xi = Inf, trendsd = 0)

Arguments

Value

The log of the prior density.

Flat prior for GEV parameters (\mu, \sigma, \xi)

Description

For information about this and other priors see set_prior.

Usage

gev_flatflat(pars, min_xi = -Inf, max_xi = Inf, trendsd = 0)

Arguments

Value

The log of the prior density.

Trivariate normal prior for GEV parameters (log \mu, log \sigma, \xi)

Description

For information about this and other priors see set_prior.

Usage

gev_loglognorm(pars, mean, icov, min_xi = -Inf, max_xi = Inf, trendsd = 0)

Arguments

Value

The log of the prior density.

Maximal data information (MDI) prior for GEV parameters (\mu, \sigma, \xi)

Description

For information about this and other priors see set_prior.

Usage

gev_mdi(pars, a = 0.577215664901532, min_xi = -1, max_xi = Inf, trendsd = 0)

Arguments

Value

The log of the prior density.

Trivariate normal prior for GEV parameters (\mu, log \sigma, \xi)

Description

For information about this and other priors see set_prior.

Usage

gev_norm(pars, mean, icov, min_xi = -Inf, max_xi = Inf, trendsd = 0)

Arguments

Value

The log of the prior density.

Informative GEV prior on a probability scale

Description

Constructs an informative prior for GEV parameters (\mu, \sigma, \xi), constructed on the probability scale. For information about how to set this prior see set_prior.

Usage

gev_prob(pars, quant, alpha, min_xi = -Inf, max_xi = Inf, trendsd = 0)

Arguments

Details

A prior for GEV parameters (\mu, \sigma, \xi), based on Crowder (1992). This construction is typically used to set an informative prior, based on specified quantiles q_1, q_2, q_3. There are two interpretations of the parameter vector alpha = (\alpha_1, \alpha_2, \alpha_3, \alpha_4): as the parameters of beta distributions for ratio of exceedance probabilities (Stephenson, 2016) and as the parameters of a Dirichlet distribution for differences between non-exceedance probabilities (Northrop et al., 2017). See these publications for details.

Value

The log of the prior density.

References

Crowder, M. (1992) Bayesian priors based on parameter transformation using the distribution function Ann. Inst. Statist. Math., 44, 405-416. https://link.springer.com/article/10.1007/BF00050695.

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. doi:10.1111/rssc.12159

Stephenson, A. (2016) Bayesian inference for extreme value modelling. In Extreme Value Modeling and Risk Analysis: Methods and Applications (eds D. K. Dey and J. Yan), 257-280, Chapman and Hall, London. doi:10.1201/b19721.

See Also

set_prior for setting a prior distribution.

rpost and rpost_rcpp for sampling from an extreme value posterior distribution.

Sets the same prior as the function prior.prob in the evdbayes package.

Informative GEV prior on a quantile scale

Description

Informative GEV prior for GEV parameters (\mu, \sigma, \xi) constructed on the quantile scale. For information about how to set this prior see set_prior.

Usage

gev_quant(pars, prob, shape, scale, min_xi = -Inf, max_xi = Inf, trendsd = 0)

Arguments

Details

See Coles and Tawn (1996) and/or Stephenson (2016) for details.

Note that the lower end point of the distribution of the distribution of the variable in question is assumed to be equal to zero. If this is not the case then the user should shift the data to ensure that this is true.

Value

The log of the prior density.

References

Coles, S. G. and Tawn, J. A. (1996) A Bayesian analysis of extreme rainfall data. Appl. Statist., 45, 463-478.

Stephenson, A. (2016) Bayesian inference for extreme value modelling. In Extreme Value Modeling and Risk Analysis: Methods and Applications (eds D. K. Dey and J. Yan), 257-280, Chapman and Hall, London. doi:10.1201/b19721.

Storm peak significant wave heights from the Gulf of Mexico

Description

A numeric vector containing 315 hindcasts of storm peak significant wave heights, metres, from 1900 to 2005 at an unnamed location in the Gulf of Mexico.

Usage

gom

Format

A vector containing 315 observations.

Source

Oceanweather Inc. (2005) GOMOS – Gulf of Mexico hindcast study.

References

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. doi:10.1111/rssc.12159

The Generalised Pareto Distribution

Description

Density function, distribution function, quantile function and random generation for the generalised Pareto (GP) distribution.

Usage

dgp(x, loc = 0, scale = 1, shape = 0, log = FALSE)

pgp(q, loc = 0, scale = 1, shape = 0, lower.tail = TRUE, log.p = FALSE)

qgp(p, loc = 0, scale = 1, shape = 0, lower.tail = TRUE, log.p = FALSE)

rgp(n, loc = 0, scale = 1, shape = 0)

Arguments

Details

The distribution function of a GP distribution with parameters location = \mu, scale = \sigma (> 0) and shape = \xi is

F(x) = 1 - [1 + \xi (x - \mu) / \sigma] ^ {-1/\xi}

for 1 + \xi (x - \mu) / \sigma > 0. If \xi = 0 the distribution function is defined as the limit as \xi tends to zero. The support of the distribution depends on \xi: it is x \geq \mu for \xi \geq 0; and \mu \leq x \leq \mu - \sigma / \xi for \xi < 0. Note that if \xi < -1 the GP density function becomes infinite as x approaches \mu - \sigma/\xi.

If lower.tail = TRUE then if p = 0 (p = 1) then the lower (upper) limit of the distribution is returned. The upper limit is Inf if shape is non-negative. Similarly, but reversed, if lower.tail = FALSE.

See https://en.wikipedia.org/wiki/Generalized_Pareto_distribution for further information.

Value

dgp gives the density function, pgp gives the distribution function, qgp gives the quantile function, and rgp generates random deviates.

References

Pickands, J. (1975) Statistical inference using extreme order statistics. Annals of Statistics, 3, 119-131. doi:10.1214/aos/1176343003

Coles, S. G. (2001) An Introduction to Statistical Modeling of Extreme Values, Springer-Verlag, London. Chapter 4: doi:10.1007/978-1-4471-3675-0_4

Examples

dgp(0:4, scale = 0.5, shape = 0.8)
dgp(1:6, scale = 0.5, shape = -0.2, log = TRUE)
dgp(1, scale = 1, shape = c(-0.2, 0.4))

pgp(0:4, scale = 0.5, shape = 0.8)
pgp(1:6, scale = 0.5, shape = -0.2)
pgp(1, scale = c(1, 2), shape = c(-0.2, 0.4))
pgp(7, scale = 1, shape = c(-0.2, 0.4))

qgp((0:9)/10, scale = 0.5, shape = 0.8)
qgp(0.5, scale = c(0.5, 1), shape = c(-0.5, 0.5))

p <- (1:9)/10
pgp(qgp(p, scale = 2, shape = 0.8), scale = 2, shape = 0.8)

rgp(6, scale = 0.5, shape = 0.8)

Beta-type prior for GP shape parameter \xi

Description

For information about this and other priors see set_prior.

Usage

gp_beta(pars, min_xi = -1/2, max_xi = 1/2, pq = c(6, 9), trendsd = 0)

Arguments

Value

The log of the prior density.

Flat prior for GP parameters (log \sigma, \xi)

Description

For information about this and other priors see set_prior.

Usage

gp_flat(pars, min_xi = -Inf, max_xi = Inf, trendsd = 0)

Arguments

Value

The log of the prior density.

Flat prior for GP parameters (\sigma, \xi)

Description

For information about this and other priors see set_prior.

Usage

gp_flatflat(pars, min_xi = -Inf, max_xi = Inf, trendsd = 0, upper = NULL)

Arguments

Value

The log of the prior density.

Jeffreys prior for GP parameters (\sigma, \xi)

Description

For information about this and other priors see set_prior.

Usage

gp_jeffreys(pars, min_xi = -1/2, max_xi = Inf, trendsd = 0)

Arguments

Value

The log of the prior density.

Linear Combinations of Ratios of Spacings estimation of generalised Pareto parameters

Description

Uses the Linear Combinations of Ratios of Spacings (LRS) methodology of (Reiss and Thomas, 2007, page 134) to estimate the parameters of the generalised Pareto (GP) distribution, based on a sample of positive values.

Usage

gp_lrs(x)

Arguments

Value

A numeric vector of length 2. The estimates of the scale parameter \sigma and the shape parameter \xi.

References

Reiss, R.-D., Thomas, M. (2007) Statistical Analysis of Extreme Values with Applications to Insurance, Finance, Hydrology and Other Fields.Birkhauser. doi:10.1007/978-3-7643-7399-3.

See Also

gp for details of the parameterisation of the GP distribution.

Examples

u <- quantile(gom, probs = 0.65)
gp_lrs((gom - u)[gom > u])

Maximal data information (MDI) prior for GP parameters (\sigma, \xi)

Description

For information about this and other priors see set_prior.

Usage

gp_mdi(pars, a = 1, min_xi = -1, max_xi = Inf, trendsd = 0)

Arguments

Value

The log of the prior density.

References

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

Bivariate normal prior for GP parameters (log \sigma, \xi)

Description

For information about this and other priors see set_prior.

Usage

gp_norm(pars, mean, icov, min_xi = -Inf, max_xi = Inf, trendsd = 0)

Arguments

Value

The log of the prior density.

Probability-weighted moments estimation of generalised Pareto parameters

Description

Uses the methodology of Hosking and Wallis (1987) to estimate the parameters of the generalised Pareto (GP) distribution.

Usage

gp_pwm(gp_data, u = 0)

Arguments

Value

A list with components

• est: A numeric vector. PWM estimates of GP parameters \sigma (scale) and \xi (shape).

• se: A numeric vector. Estimated standard errors of \sigma and \xi.

• cov: A numeric matrix. Estimate covariance matrix of the the PWM estimators of \sigma and \xi.

References

Hosking, J. R. M. and Wallis, J. R. (1987) Parameter and Quantile Estimation for the Generalized Pareto Distribution. Technometrics, 29(3), 339-349. doi:10.2307/1269343.

See Also

gp for details of the parameterisation of the GP distribution.

Examples

u <- quantile(gom, probs = 0.65)
gp_pwm(gom, u)

Maximum likelihood estimation of generalised Pareto parameters

Description

Uses the methodology of Grimshaw (1993) to find the MLEs of the parameters of the generalised Pareto distribution, based on a sample of positive values. The function is essentially the same as that made available with Grimshaw (1993), with only minor modifications.

Usage

grimshaw_gp_mle(x)

Arguments

Value

A numeric vector of length 2. The estimates of the negated shape parameter k (= -\xi) and the scale parameter a (= \sigma).

References

Grimshaw, S. D. (1993) Computing Maximum Likelihood Estimates for the Generalized Pareto Distribution. Technometrics, 35(2), 185-191. and Computing (1991) 1, 129-133. doi:10.1080/00401706.1993.10485040.

See Also

gp for details of the parameterisation of the GP distribution, in terms of \sigma and \xi.

Examples

u <- quantile(gom, probs = 0.65)
grimshaw_gp_mle((gom - u)[gom > u])

Random sampling from K-gaps posterior distribution

Description

Uses the rust package to simulate from the posterior distribution of the extremal index \theta based on the K-gaps model for threshold interexceedance times of Suveges and Davison (2010).

Usage

kgaps_post(
  data,
  thresh,
  k = 1,
  n = 1000,
  inc_cens = TRUE,
  alpha = 1,
  beta = 1,
  param = c("logit", "theta"),
  use_rcpp = TRUE
)

Arguments

Details

A beta(\alpha, \beta) prior distribution is used for \theta so that the posterior from which values are simulated is proportional to

\theta ^ {2 N_1 + \alpha - 1} (1 - \theta) ^ {N_0 + \beta - 1} \exp\{- \theta q (S_0 + \cdots + S_N)\}.

See kgaps_stat for a description of the variables involved in the contribution of the likelihood to this expression.

The ru function in the rust package simulates from this posterior distribution using the generalised ratio-of-uniforms distribution. To improve the probability of acceptance, and to ensure that the simulation will work even in extreme cases where the posterior density of \theta is unbounded as \theta approaches 0 or 1, we simulate from the posterior distribution of \phi = \log(\theta / (1-\theta)) and then transform back to the \theta-scale.

Value

An object (list) of class "evpost", which has the same structure as an object of class "ru" returned from ru. In addition this list contains

• call: The call to kgaps().

• model: The character scalar "kgaps".

• thresh: The argument thresh.

• ss: The sufficient statistics for the K-gaps likelihood, as calculated by kgaps_stat.

References

Suveges, M. and Davison, A. C. (2010) Model misspecification in peaks over threshold analysis, The Annals of Applied Statistics, 4(1), 203-221. doi:10.1214/09-AOAS292

See Also

ru for the form of the object returned by kgaps_post.

dgaps_post for Bayesian inference about the extremal index \theta using the D-gaps model.

Examples

# Newlyn sea surges

thresh <- quantile(newlyn, probs = 0.90)
k_postsim <- kgaps_post(newlyn, thresh)
plot(k_postsim)

### Cheeseboro wind gusts

k_postsim <- kgaps_post(exdex::cheeseboro, thresh = 45, k = 3)
plot(k_postsim)

Newlyn sea surges

Description

The vector newlyn contains 2894 maximum sea-surges measured at Newlyn, Cornwall, UK over the period 1971-1976. The observations are the maximum hourly sea-surge heights over contiguous 15-hour time periods.

Usage

newlyn

Format

A vector of length 2894.

Source

Coles, S.G. (1991) Modelling extreme multivariate events. PhD thesis, University of Sheffield, U.K.

References

Fawcett, L. and Walshaw, D. (2012) Estimating return levels from serially dependent extremes. Environmetrics, 23(3), 272-283. doi:10.1002/env.2133

Northrop, P. J. (2015) An efficient semiparametric maxima estimator of the extremal index. Extremes, 18, 585-603. doi:10.1007/s10687-015-0221-5

Annual Maximum Temperatures at Oxford

Description

A numeric vector containing annual maximum temperatures, in degrees Fahrenheit, from 1901 to 1980 at Oxford, England.

Usage

oxford

Format

A vector containing 80 observations.

Source

Tabony, R. C. (1983) Extreme value analysis in meteorology. The Meteorological Magazine, 112, 77-98.

Plot diagnostics for an evpost object

Description

plot method for class "evpost". For d = 1 a histogram of the simulated values is plotted with a the density function superimposed. The density is normalized crudely using the trapezium rule. For d = 2 a scatter plot of the simulated values is produced with density contours superimposed. For d > 2 pairwise plots of the simulated values are produced. An interface is also provided to the functions in the bayesplot package that produce plots of Markov chain Monte Carlo (MCMC) simulations. See MCMC-overview for details of these functions.

Usage

## S3 method for class 'evpost'
plot(
  x,
  y,
  ...,
  n = ifelse(x$d == 1, 1001, 101),
  prob = c(0.5, 0.1, 0.25, 0.75, 0.95, 0.99),
  ru_scale = FALSE,
  rows = NULL,
  xlabs = NULL,
  ylabs = NULL,
  points_par = list(col = 8),
  pu_only = FALSE,
  add_pu = FALSE,
  use_bayesplot = FALSE,
  fun_name = c("areas", "intervals", "dens", "hist", "scatter")
)

Arguments

Details

For details of the bayesplot functions available when use_bayesplot = TRUE see MCMC-overview and the bayesplot vignette Plotting MCMC draws.

Value

Nothing is returned unless use_bayesplot = TRUE when a ggplot object, which can be further customized using the ggplot2 package, is returned.

References

Jonah Gabry (2016). bayesplot: Plotting for Bayesian Models. R package version 1.1.0. https://CRAN.R-project.org/package=bayesplot

See Also

summary.evpost for summaries of the simulated values and properties of the ratio-of-uniforms algorithm.

MCMC-overview, MCMC-intervals, MCMC-distributions.

Examples

## GP posterior
u <- stats::quantile(gom, probs = 0.65)
fp <- set_prior(prior = "flat", model = "gp", min_xi = -1)
gpg <- rpost(n = 1000, model = "gp", prior = fp, thresh = u, data = gom)
plot(gpg)


# Using the bayesplot package
plot(gpg, use_bayesplot = TRUE)
plot(gpg, use_bayesplot = TRUE, pars = "xi", prob = 0.95)
plot(gpg, use_bayesplot = TRUE, fun_name = "intervals", pars = "xi")
plot(gpg, use_bayesplot = TRUE, fun_name = "hist")
plot(gpg, use_bayesplot = TRUE, fun_name = "dens")
plot(gpg, use_bayesplot = TRUE, fun_name = "scatter")


## bin-GP posterior
u <- quantile(gom, probs = 0.65)
fp <- set_prior(prior = "flat", model = "gp", min_xi = -1)
bp <- set_bin_prior(prior = "jeffreys")
npy_gom <- length(gom)/105
bgpg <- rpost(n = 1000, model = "bingp", prior = fp, thresh = u,
              data = gom, bin_prior = bp, npy = npy_gom)
plot(bgpg)
plot(bgpg, pu_only = TRUE)
plot(bgpg, add_pu = TRUE)


# Using the bayesplot package
dimnames(bgpg$bin_sim_vals)
plot(bgpg, use_bayesplot = TRUE)
plot(bgpg, use_bayesplot = TRUE, fun_name = "hist")
plot(bgpg, use_bayesplot = TRUE, pars = "p[u]")

Plot diagnostics for an evpred object

Description

plot method for class "evpred". Plots summarising the predictive distribution of the largest value to be observed in N years are produced. The plot produced depends on x$type. If x$type = "d", "p" or "q" then matplot is used to produce a line plot of the predictive density, distribution or quantile function, respectively, with a line for each value of N in x$n_years. If x$type = "r" then estimates of the predictive density (from density) are plotted with a line for each N. If x$type = "i" then lines representing estimated predictive intervals are plotted, with the level of the interval indicated next to the line.

Usage

## S3 method for class 'evpred'
plot(
  x,
  ...,
  leg_pos = NULL,
  leg_text = NULL,
  which_int = c("long", "short", "both")
)

Arguments

Value

Nothing is returned.

See Also

predict.evpost for the S3 predict method for objects of class evpost.

Examples

data(portpirie)
mat <- diag(c(10000, 10000, 100))
pn <- set_prior(prior = "norm", model = "gev", mean = c(0,0,0), cov = mat)
gevp  <- rpost(n = 1000, model = "gev", prior = pn, data = portpirie)

# Predictive density function
d_gevp <- predict(gevp, type = "d", n_years = c(100, 1000))
plot(d_gevp)

# Predictive distribution function
p_gevp <- predict(gevp, type = "p", n_years = c(100, 1000))
plot(p_gevp)

# Predictive quantiles
q_gevp <- predict(gevp, type = "q", n_years = c(100, 1000))
plot(q_gevp)

# Predictive intervals
i_gevp <- predict(gevp, type = "i", n_years = c(100, 1000), hpd = TRUE)
plot(i_gevp, which_int = "both")

# Sample from predictive distribution
r_gevp <- predict(gevp, type = "r", n_years = c(100, 1000))
plot(r_gevp)
plot(r_gevp, xlim = c(4, 10))

Annual Maximum Sea Levels at Port Pirie, South Australia

Description

A numeric vector of length 65 containing annual maximum sea levels, in metres, from 1923 to 1987 at Port Pirie, South Australia.

Usage

portpirie

Format

A numeric vector containing 65 observations.

Source

Coles, S. G. (2001) An Introduction to Statistical Modelling of Extreme Values. London: Springer. doi:10.1007/978-1-4471-3675-0

Posterior predictive checks for an evpost object

Description

pp_check method for class "evpost". This provides an interface to the functions that perform posterior predictive checks in the bayesplot package. See PPC-overview for details of these functions.

Usage

## S3 method for class 'evpost'
pp_check(
  object,
  ...,
  type = c("stat", "overlaid", "multiple", "intervals", "user"),
  subtype = NULL,
  stat = "median",
  nrep = 8,
  fun = NULL
)

Arguments

Details

For details of these functions see PPC-overview. See also the vignette Posterior Predictive Extreme Value Inference and the bayesplot vignette Graphical posterior predictive checks.

The general idea is to compare the observed data object$data with a matrix object$data_rep in which each row is a replication of the observed data simulated from the posterior predictive distribution. For greater detail see Chapter 6 of Gelman et al. (2013).

The format of object$data depends on the model:

• model = "gev". A vector of block maxima.

• model = "gp". Data that lie above the threshold, i.e. threshold exceedances.

• model = "bingp" or model = "pp". The input data are returned but any value lying below the threshold is set to object$thresh.

In all cases any missing values have been removed from the data.

If model = "bingp" or "pp" the rate of threshold exceedance is part of the inference. Therefore, the number of values in object$data_rep that lie above the threshold varies between predictive replications, with values below the threshold being left-censored at the threshold. This limits a little the posterior predictive checks that it is useful to perform. In the examples below we have compared object$data and object$data_rep using only their sample maxima.

Value

A ggplot object that can be further customized using the ggplot2 package.

References

Jonah Gabry (2016). bayesplot: Plotting for Bayesian Models. R package version 1.1.0. https://CRAN.R-project.org/package=bayesplot

Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., and Rubin, D. B. (2013). Bayesian Data Analysis. Chapman & Hall/CRC Press, London, third edition. (Chapter 6)

See Also

rpost and rpost_rcpp for sampling from an extreme value posterior distribution.

bayesplot functions PPC-overview, PPC-distributions, PPC-test-statistics, PPC-intervals, pp_check.

Examples

# GEV model
data(portpirie)
mat <- diag(c(10000, 10000, 100))
pn <- set_prior(prior = "norm", model = "gev", mean = c(0,0,0), cov = mat)
gevp  <- rpost(1000, model = "gev", prior = pn, data = portpirie,
               nrep = 50)

# Posterior predictive test statistics
pp_check(gevp)
pp_check(gevp, stat = "min")
pp_check(gevp, stat = c("min", "max"))
iqr <- function(y) diff(quantile(y, c(0.25, 0.75)))
pp_check(gevp, stat = "iqr")

# Overlaid density and distributions functions
pp_check(gevp, type = "overlaid")
pp_check(gevp, type = "overlaid", subtype = "dens")

# Multiple plots
pp_check(gevp, type = "multiple")
pp_check(gevp, type = "multiple", subtype = "hist")
pp_check(gevp, type = "multiple", subtype = "boxplot")

# Intervals
pp_check(gevp, type = "intervals")
pp_check(gevp, type = "intervals", subtype = "ribbon")

# User-supplied bayesplot function
# Equivalent to p_check(gevp, type = "overlaid")
pp_check(gevp, type = "user", fun = "dens_overlay")

# GP model
u <- quantile(gom, probs = 0.65)
fp <- set_prior(prior = "flat", model = "gp", min_xi = -1)
gpg <- rpost(n = 1000, model = "gp", prior = fp, thresh = u,
             data = gom, nrep = 50)
pp_check(gpg)
pp_check(gpg, type = "overlaid")

# bin-GP model
bp <- set_bin_prior(prior = "jeffreys")
bgpg <- rpost(n = 1000, model = "bingp", prior = fp, thresh = u,
              data = gom, bin_prior = bp, nrep = 50)
pp_check(bgpg, stat = "max")

# PP model
data(rainfall)
rthresh <- 40
pf <- set_prior(prior = "flat", model = "gev", min_xi = -1)
ppr <- rpost(n = 1000, model = "pp", prior = pf, data = rainfall,
             thresh = rthresh, noy = 54, nrep = 50)
pp_check(ppr, stat = "max")

Predictive inference for the largest value observed in N years.

Description

predict method for class "evpost". Performs predictive inference about the largest value to be observed over a future time period of N years. Predictive inferences accounts for uncertainty in model parameters and for uncertainty owing to the variability of future observations.

Usage

## S3 method for class 'evpost'
predict(
  object,
  type = c("i", "p", "d", "q", "r"),
  x = NULL,
  x_num = 100,
  n_years = 100,
  npy = NULL,
  level = 95,
  hpd = FALSE,
  lower_tail = TRUE,
  log = FALSE,
  big_q = 1000,
  ...
)

Arguments

Details

Inferences about future extreme observations are integrated over the posterior distribution of the model parameters, thereby accounting for uncertainty in model parameters and uncertainty owing to the variability of future observations. In practice the integrals involved are estimated using an empirical mean over the posterior sample. See, for example, Coles (2001), Stephenson (2016) or Northrop et al. (2017) for details. See also the vignette Posterior Predictive Extreme Value Inference

GEV / OS / PP. If model = "gev", model = "os" or model = "pp" in the call to rpost or rpost_rcpp we first calculate the number of blocks b in n_years years. To calculate the density function or distribution function of the maximum over n_years we call dgev or pgev with m = b.

• type = "p". We calculate using pgev the GEV distribution function at q for each of the posterior samples of the location, scale and shape parameters. Then we take the mean of these values.

• type = "d". We calculate using dgev the GEV density function at x for each of the posterior samples of the location, scale and shape parameters. Then we take the mean of these values.

• type = "q". We solve numerically predict.evpost(object, type = "p", x = q) = p[i] numerically for q for each element p[i] of p.

• type = "i". If hpd = FALSE then the interval is equi-tailed, equal to predict.evpost() object, type ="q", x = p), where p = c((1-level/100)/2, (1+level/100)/2). If hpd = TRUE then, in addition, we perform a numerical minimisation of the length of level% intervals, after approximating the predictive quantile function using monotonic cubic splines, to reduce computing time.

• type = "r". For each simulated value of the GEV parameters at the n_years level of aggregation we simulate one value from this GEV distribution using rgev. Thus, each sample from the predictive distribution is of a size equal to the size of the posterior sample.

Binomial-GP. If model = "bingp" in the call to rpost or rpost_rcpp then we calculate the mean number of observations in n_years years, i.e. npy * n_years.

Following Northrop et al. (2017), let M_N be the largest value observed in N years, m = npy * n_years and u the threshold object$thresh used in the call to rpost or rpost_rcpp. For fixed values of \theta = (p, \sigma, \xi) the distribution function of M_N is given by F(z, \theta)^m, for z \geq u, where

F(z, \theta) = 1 - p [1 + \xi (x - u) / \sigma] ^ {-1/\xi}.

The distribution function of M_N cannot be evaluated for z < u because no model has been supposed for observations below the threshold.

• type = "p". We calculate F(z, \theta)^m at q for each of the posterior samples \theta. Then we take the mean of these values.

• type = "d". We calculate the density of of M_n, i.e. the derivative of F(z, \theta)^m with respect to z at x for each of the posterior samples \theta. Then we take the mean of these values.

• type = "q" and type = "i". We perform calculations that are analogous to the GEV case above. If n_years is very small and/or level is very close to 100 then a predictive interval may extend below the threshold. In such cases NAs are returned (see Value below).

• type = "r". For each simulated value of the bin-GP parameter we simulate from the distribution of M_N using the inversion method applied to the distribution function of M_N given above. Occasionally a value below the threshold would need to be simulated. If these instances a missing value code NA is returned. Thus, each sample from the predictive distribution is of a size equal to the size of the posterior sample, perhaps with a small number os NAs.

Value

An object of class "evpred", a list containing a subset of the following components:

The arguments n_years, level, hpd, lower_tail, log supplied to predict.evpost are also included, as is the argument npy supplied to, or set within, predict.evpost and the arguments data and model from the original call to rpost or rpost_rcpp.

References

Coles, S. G. (2001) An Introduction to Statistical Modeling of Extreme Values, Springer-Verlag, London. Chapter 9: doi:10.1007/978-1-4471-3675-0_9

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. doi:10.1111/rssc.12159

Stephenson, A. (2016). Bayesian Inference for Extreme Value Modelling. In Extreme Value Modeling and Risk Analysis: Methods and Applications, edited by D. K. Dey and J. Yan, 257-80. London: Chapman and Hall. doi:10.1201/b19721

See Also

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

rpost or rpost_rcpp for sampling from an extreme value posterior distribution.

Examples

### GEV
data(portpirie)
mat <- diag(c(10000, 10000, 100))
pn <- set_prior(prior = "norm", model = "gev", mean = c(0,0,0), cov = mat)
gevp  <- rpost_rcpp(n = 1000, model = "gev", prior = pn, data = portpirie)

# Interval estimation
predict(gevp)$long
predict(gevp, hpd = TRUE)$short
# Density function
x <- 4:7
predict(gevp, type = "d", x = x)$y
plot(predict(gevp, type = "d", n_years = c(100, 1000)))
# Distribution function
predict(gevp, type = "p", x = x)$y
plot(predict(gevp, type = "p", n_years = c(100, 1000)))
# Quantiles
predict(gevp, type = "q", n_years = c(100, 1000))$y
# Random generation
plot(predict(gevp, type = "r"))

### Binomial-GP
u <- quantile(gom, probs = 0.65)
fp <- set_prior(prior = "flat", model = "gp", min_xi = -1)
bp <- set_bin_prior(prior = "jeffreys")
npy_gom <- length(gom)/105
bgpg <- rpost_rcpp(n = 1000, model = "bingp", prior = fp, thresh = u,
                   data = gom, bin_prior = bp)

# Setting npy in call to predict.evpost()
predict(bgpg, npy = npy_gom)$long

# Setting npy in call to rpost() or rpost_rcpp()
bgpg <- rpost_rcpp(n = 1000, model = "bingp", prior = fp, thresh = u,
                   data = gom, bin_prior = bp, npy = npy_gom)

# Interval estimation
predict(bgpg)$long
predict(bgpg, hpd = TRUE)$short
# Density function
plot(predict(bgpg, type = "d", n_years = c(100, 1000)))
# Distribution function
plot(predict(bgpg, type = "p", n_years = c(100, 1000)))
# Quantiles
predict(bgpg, type = "q", n_years = c(100, 1000))$y
# Random generation
plot(predict(bgpg, type = "r"))

Print method for objects of class "evpost"

Description

Print method for objects of class "evpost"

Usage

## S3 method for class 'evpost'
print(x, ...)

Arguments

Details

print.evpost just prints the original function call, to avoid printing a huge list.

Value

The argument x is returned, invisibly.

See Also

plot.evpost for a diagnostic plot.

Examples

# Newlyn sea surges

thresh <- quantile(newlyn, probs = 0.90)
k_postsim <- kgaps_post(newlyn, thresh)
k_postsim

Print method for objects of class "summary.evpost"

Description

print method for an object object of class "summary.evpost".

Usage

## S3 method for class 'summary.evpost'
print(x, ...)

Arguments

Value

Prints

• information about the ratio-of-uniforms bounding box, i.e. object$box

• an estimate of the probability of acceptance, i.e. object$pa

• a summary of the simulated values, via summary(object$sim_vals)

See Also

ru or ru_rcpp for descriptions of object$sim_vals and $box.

plot.evpost for a diagnostic plot.

Examples

# GP posterior
u <- stats::quantile(gom, probs = 0.65)
fp <- set_prior(prior = "flat", model = "gp", min_xi = -1)
gpg <- rpost_rcpp(n = 1000, model = "gp", prior = fp, thresh = u,
                  data = gom)
summary(gpg)

Converts quantiles to GEV parameters

Description

Three quantiles, that is, the value of quantile and their respective exceedance probabilities, are provided. This function attempts to find the location, scale and shape parameters of a GEV distribution that has these quantiles.

Usage

quantile_to_gev(quant, prob)

Arguments

Details

Suppose that G(x) is the distribution function of a GEV(\mu, \sigma, \xi) distribution. This function attempts to solve numerically the set of three non-linear equations

G(q_i) = 1 - p_i, i = 1, 2, 3

where q_i, i=1,2,3 are the quantiles in quant and p_i, i=1,2,3 are the exceedance probabilities in prob. This is reduced to a one-dimensional optimisation over the GEV shape parameter.

Value

A numeric vector of length 3 containing the GEV location, scale and shape parameters.

See Also

rprior_quant for simulation of GEV parameters from a prior constructed on the quantile scale.

Examples

my_q <- c(15, 20, 22.5)
my_p <- 1-c(0.5, 0.9, 0.5^0.01)
x <- quantile_to_gev(quant = my_q, prob = my_p)
# Check
qgev(p = 1 - my_p, loc = x[1], scale = x[2], shape = x[3])

Daily Aggregate Rainfall

Description

A numeric vector of length 20820 containing daily aggregate rainfall observations, in millimetres, recorded at a rain gauge in England over a 57 year period, beginning on a leap year. Three of these years contain only missing values.

Usage

rainfall

Format

A vector containing 20820 observations.

Source

Unknown

Simulation from a Dirichlet distribution

Description

Simulates from a Dirichlet distribution with concentration parameter vector \alpha = (\alpha_1, ..., \alpha_K).

Usage

rDir(n = 1, alpha = c(1, 1))

Arguments

Details

The simulation is based on the property that if Y_1, \ldots, Y_K are independent, Y_i has a gamma(\alpha_i, 1) distribution and S = Y_1 + \cdots + Y_k then (Y_1, \ldots, Y_K) / S has a Dirichlet(\alpha_1, ..., \alpha_K) distribution.

See https://en.wikipedia.org/wiki/Dirichlet_distribution#Gamma_distribution

Value

An n by length(alpha) numeric matrix.

References

Kotz, S., Balakrishnan, N. and Johnson, N. L. (2000) Continuous Multivariate Distributions, vol. 1, Models and Applications, 2nd edn, ch. 49. New York: Wiley. doi:10.1002/0471722065

See Also

rprior_prob for prior simulation of GEV parameters - prior on probability scale.

Examples

rDir(n = 10, alpha = 1:4)

Internal revdbayes functions

Description

Internal revdbayes functions

Usage

process_data(model, data, thresh, noy, use_noy, ros, weights = NULL)

create_ru_list(model, trans, rotate, min_xi, max_xi)

set_which_lam(model)

set_range_phi(model, phi_mid, se_phi, mult)

box_cox(x, lambda = 1, gm = 1, lambda_tol = 1e-06, poly_order = 3)

box_cox_vec(x, lambda = 1, lambda_tol = 1e-06)

box_cox_deriv(x, lambda = 1, lambda_tol = 1e-06, poly_order = 3)

check_sample_size(prior_name, n_check)

check_sample_size_message(prior_name, n_check)

logNegNA(x)

fallback_gp_mle(init, ...)

Details

These functions are not intended to be called by the user.

Random sampling from extreme value posterior distributions

Description

Uses the ru function in the rust package to simulate from the posterior distribution of an extreme value model.

Usage

rpost(
  n,
  model = c("gev", "gp", "bingp", "pp", "os"),
  data,
  prior,
  ...,
  nrep = NULL,
  thresh = NULL,
  noy = NULL,
  use_noy = TRUE,
  npy = NULL,
  ros = NULL,
  bin_prior = structure(list(prior = "bin_beta", ab = c(1/2, 1/2), class = "binprior")),
  bin_param = "logit",
  init_ests = NULL,
  mult = 2,
  use_phi_map = FALSE,
  weights = NULL
)

Arguments

Details

Generalised Pareto (GP): model = "gp". A model for threshold excesses. Required arguments: n, data and prior. If thresh is supplied then only the values in data that exceed thresh are used and the GP distribution is fitted to the amounts by which those values exceed thresh. If thresh is not supplied then the GP distribution is fitted to all values in data, in effect thresh = 0. See also gp.

Binomial-GP: model = "bingp". The GP model for threshold excesses supplemented by a binomial(length(data), p) model for the number of threshold excesses. See Northrop et al. (2017) for details. Currently, the GP and binomial parameters are assumed to be independent a priori.

Generalised extreme value (GEV) model: model = "gev". A model for block maxima. Required arguments: n, data, prior. See also gev.

Point process (PP) model: model = "pp". A model for occurrences of threshold exceedances and threshold excesses. Required arguments: n, data, prior, thresh and noy.

r-largest order statistics (OS) model: model = "os". A model for the largest order statistics within blocks of data. Required arguments: n, data, prior. All the values in data are used unless ros is supplied.

Parameter transformation. The scalar logical arguments (to the function ru) trans and rotate determine, respectively, whether or not Box-Cox transformation is used to reduce asymmetry in the posterior distribution and rotation of parameter axes is used to reduce posterior parameter dependence. The default is trans = "none" and rotate = TRUE.

See the Introducing revdbayes vignette for further details and examples.

Value

An object (list) of class "evpost", which has the same structure as an object of class "ru" returned from ru. In addition this list contains

If nrep is not NULL then this list also contains data_rep, a numerical matrix with nrep rows. Each row contains a replication of the original data data simulated from the posterior predictive distribution. If model = "bingp" or "pp" then the rate of threshold exceedance is part of the inference. Therefore, the number of values in data_rep that lie above the threshold varies between predictive replications (different rows of data_rep). Values below the threshold are left-censored at the threshold, i.e. they are set at the threshold.

If model == "pp" then this list also contains the argument noy to rpost detailed above. If the argument npy was supplied then this list also contains npy.

If model == "gp" or model == "bingp" then this list also contains the argument thresh to rpost detailed above.

If model == "bingp" then this list also contains

References

Coles, S. G. and Powell, E. A. (1996) Bayesian methods in extreme value modelling: a review and new developments. Int. Statist. Rev., 64, 119-136.

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. doi:10.1111/rssc.12159

Stephenson, A. (2016) Bayesian Inference for Extreme Value Modelling. In Extreme Value Modeling and Risk Analysis: Methods and Applications, edited by D. K. Dey and J. Yan, 257-80. London: Chapman and Hall. doi:10.1201/b19721 value posterior using the evdbayes package.

Wadsworth, J. L., Tawn, J. A. and Jonathan, P. (2010) Accounting for choice of measurement scale in extreme value modeling. The Annals of Applied Statistics, 4(3), 1558-1578. doi:10.1214/10-AOAS333

See Also

set_prior for setting a prior distribution.

rpost_rcpp for faster posterior simulation using the Rcpp package.

plot.evpost, summary.evpost and predict.evpost for the S3 plot, summary and predict methods for objects of class evpost.

ru and ru_rcpp in the rust package for details of the arguments that can be passed to ru and the form of the object returned by rpost.

find_lambda and find_lambda_rcpp in the rust package is used inside rpost to set the Box-Cox transformation parameter lambda when the trans = "BC" argument is given.

Examples

# GP model
u <- quantile(gom, probs = 0.65)
fp <- set_prior(prior = "flat", model = "gp", min_xi = -1)
gpg <- rpost(n = 1000, model = "gp", prior = fp, thresh = u, data = gom)
plot(gpg)

# Binomial-GP model
u <- quantile(gom, probs = 0.65)
fp <- set_prior(prior = "flat", model = "gp", min_xi = -1)
bp <- set_bin_prior(prior = "jeffreys")
bgpg <- rpost(n = 1000, model = "bingp", prior = fp, thresh = u, data = gom,
              bin_prior = bp)
plot(bgpg, pu_only = TRUE)
plot(bgpg, add_pu = TRUE)

# Setting the same binomial (Jeffreys) prior by hand
beta_prior_fn <- function(p, ab) {
  return(stats::dbeta(p, shape1 = ab[1], shape2 = ab[2], log = TRUE))
}
jeffreys <- set_bin_prior(beta_prior_fn, ab = c(1 / 2, 1 / 2))
bgpg <- rpost(n = 1000, model = "bingp", prior = fp, thresh = u, data = gom,
              bin_prior = jeffreys)
plot(bgpg, pu_only = TRUE)
plot(bgpg, add_pu = TRUE)

# GEV model
mat <- diag(c(10000, 10000, 100))
pn <- set_prior(prior = "norm", model = "gev", mean = c(0, 0, 0), cov = mat)
gevp  <- rpost(n = 1000, model = "gev", prior = pn, data = portpirie)
plot(gevp)

# GEV model, informative prior constructed on the probability scale
pip  <- set_prior(quant = c(85, 88, 95), alpha = c(4, 2.5, 2.25, 0.25),
                  model = "gev", prior = "prob")
ox_post <- rpost(n = 1000, model = "gev", prior = pip, data = oxford)
plot(ox_post)

# PP model
pf <- set_prior(prior = "flat", model = "gev", min_xi = -1)
ppr <- rpost(n = 1000, model = "pp", prior = pf, data = rainfall,
             thresh = 40, noy = 54)
plot(ppr)

# PP model, informative prior constructed on the quantile scale
piq <- set_prior(prob = 10^-(1:3), shape = c(38.9, 7.1, 47),
                 scale = c(1.5, 6.3, 2.6), model = "gev", prior = "quant")
rn_post <- rpost(n = 1000, model = "pp", prior = piq, data = rainfall,
                 thresh = 40, noy = 54)
plot(rn_post)

# OS model
mat <- diag(c(10000, 10000, 100))
pv <- set_prior(prior = "norm", model = "gev", mean = c(0, 0, 0), cov = mat)
osv <- rpost(n = 1000, model = "os", prior = pv, data = venice)
plot(osv)

Random sampling from extreme value posterior distributions

Description

Uses the ru_rcpp function in the rust package to simulate from the posterior distribution of an extreme value model.

Usage

rpost_rcpp(
  n,
  model = c("gev", "gp", "bingp", "pp", "os"),
  data,
  prior,
  ...,
  nrep = NULL,
  thresh = NULL,
  noy = NULL,
  use_noy = TRUE,
  npy = NULL,
  ros = NULL,
  bin_prior = structure(list(prior = "bin_beta", ab = c(1/2, 1/2), class = "binprior")),
  init_ests = NULL,
  mult = 2,
  use_phi_map = FALSE
)

Arguments

Details

Generalised Pareto (GP): model = "gp". A model for threshold excesses. Required arguments: n, data and prior. If thresh is supplied then only the values in data that exceed thresh are used and the GP distribution is fitted to the amounts by which those values exceed thresh. If thresh is not supplied then the GP distribution is fitted to all values in data, in effect thresh = 0. See also gp.

Binomial-GP: model = "bingp". The GP model for threshold excesses supplemented by a binomial(length(data), p) model for the number of threshold excesses. See Northrop et al. (2017) for details. Currently, the GP and binomial parameters are assumed to be independent a priori.

Generalised extreme value (GEV) model: model = "gev". A model for block maxima. Required arguments: n, data, prior. See also gev.

Point process (PP) model: model = "pp". A model for occurrences of threshold exceedances and threshold excesses. Required arguments: n, data, prior, thresh and noy.

r-largest order statistics (OS) model: model = "os". A model for the largest order statistics within blocks of data. Required arguments: n, data, prior. All the values in data are used unless ros is supplied.

Parameter transformation. The scalar logical arguments (to the function ru) trans and rotate determine, respectively, whether or not Box-Cox transformation is used to reduce asymmetry in the posterior distribution and rotation of parameter axes is used to reduce posterior parameter dependence. The default is trans = "none" and rotate = TRUE.

See the Introducing revdbayes vignette for further details and examples.

Value

An object (list) of class "evpost", which has the same structure as an object of class "ru" returned from ru_rcpp. In addition this list contains

If nrep is not NULL then this list also contains data_rep, a numerical matrix with nrep rows. Each row contains a replication of the original data data simulated from the posterior predictive distribution. If model = "bingp" or "pp" then the rate of threshold exceedance is part of the inference. Therefore, the number of values in data_rep that lie above the threshold varies between predictive replications (different rows of data_rep). Values below the threshold are left-censored at the threshold, i.e. they are set at the threshold.

If model == "pp" then this list also contains the argument noy to rpost detailed above. If the argument npy was supplied then this list also contains npy.

If model == "gp" or model == "bingp" then this list also contains the argument thresh to rpost detailed above.

If model == "bingp" then this list also contains

References

Coles, S. G. and Powell, E. A. (1996) Bayesian methods in extreme value modelling: a review and new developments. Int. Statist. Rev., 64, 119-136.

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. doi:10.1111/rssc.12159

Stephenson, A. (2016) Bayesian Inference for Extreme Value Modelling. In Extreme Value Modeling and Risk Analysis: Methods and Applications, edited by D. K. Dey and J. Yan, 257-80. London: Chapman and Hall. doi:10.1201/b19721 value posterior using the evdbayes package.

Wadsworth, J. L., Tawn, J. A. and Jonathan, P. (2010) Accounting for choice of measurement scale in extreme value modeling. The Annals of Applied Statistics, 4(3), 1558-1578. doi:10.1214/10-AOAS333

See Also

set_prior for setting a prior distribution.

rpost for posterior simulation without using the Rcpp package.

plot.evpost, summary.evpost and predict.evpost for the S3 plot, summary and predict methods for objects of class evpost.

ru_rcpp in the rust package for details of the arguments that can be passed to ru_rcpp and the form of the object returned by rpost_rcpp.

find_lambda in the rust package is used inside rpost to set the Box-Cox transformation parameter lambda when the trans = "BC" argument is given.

Examples

# GP model
u <- quantile(gom, probs = 0.65)
fp <- set_prior(prior = "flat", model = "gp", min_xi = -1)
gpg <- rpost_rcpp(n = 1000, model = "gp", prior = fp, thresh = u,
                  data = gom)
plot(gpg)

# GP model, user-defined prior (same prior as the previous example)
ptr_gp_flat <- create_prior_xptr("gp_flat")
p_user <- set_prior(prior = ptr_gp_flat, model = "gp", min_xi = -1)
gpg <- rpost_rcpp(n = 1000, model = "gp", prior = p_user, thresh = u,
                  data = gom)
plot(gpg)

# Binomial-GP model
u <- quantile(gom, probs = 0.65)
fp <- set_prior(prior = "flat", model = "gp", min_xi = -1)
bp <- set_bin_prior(prior = "jeffreys")
bgpg <- rpost_rcpp(n = 1000, model = "bingp", prior = fp, thresh = u,
                   data = gom, bin_prior = bp)
plot(bgpg, pu_only = TRUE)
plot(bgpg, add_pu = TRUE)

# Setting the same binomial (Jeffreys) prior by hand
beta_prior_fn <- function(p, ab) {
  return(stats::dbeta(p, shape1 = ab[1], shape2 = ab[2], log = TRUE))
}
jeffreys <- set_bin_prior(beta_prior_fn, ab = c(1 / 2, 1 / 2))
bgpg <- rpost_rcpp(n = 1000, model = "bingp", prior = fp, thresh = u,
                   data = gom, bin_prior = jeffreys)
plot(bgpg, pu_only = TRUE)
plot(bgpg, add_pu = TRUE)

# GEV model
mat <- diag(c(10000, 10000, 100))
pn <- set_prior(prior = "norm", model = "gev", mean = c(0, 0, 0), cov = mat)
gevp  <- rpost_rcpp(n = 1000, model = "gev", prior = pn, data = portpirie)
plot(gevp)

# GEV model, user-defined prior (same prior as the previous example)
mat <- diag(c(10000, 10000, 100))
ptr_gev_norm <- create_prior_xptr("gev_norm")
pn_u <- set_prior(prior = ptr_gev_norm, model = "gev", mean = c(0, 0, 0),
                  icov = solve(mat))
gevu <- rpost_rcpp(n = 1000, model = "gev", prior = pn_u, data = portpirie)
plot(gevu)

# GEV model, informative prior constructed on the probability scale
pip  <- set_prior(quant = c(85, 88, 95), alpha = c(4, 2.5, 2.25, 0.25),
  model = "gev", prior = "prob")
ox_post <- rpost_rcpp(n = 1000, model = "gev", prior = pip, data = oxford)
plot(ox_post)

# PP model
pf <- set_prior(prior = "flat", model = "gev", min_xi = -1)
ppr <- rpost_rcpp(n = 1000, model = "pp", prior = pf, data = rainfall,
                  thresh = 40, noy = 54)
plot(ppr)

# PP model, user-defined prior (same prior as the previous example)
ptr_gev_flat <- create_prior_xptr("gev_flat")
pf_u <- set_prior(prior = ptr_gev_flat, model = "gev", min_xi = -1,
                  max_xi = Inf)
ppru <- rpost_rcpp(n = 1000, model = "pp", prior = pf_u, data = rainfall,
                   thresh = 40, noy = 54)
plot(ppru)

# PP model, informative prior constructed on the quantile scale
piq <- set_prior(prob = 10^-(1:3), shape = c(38.9, 7.1, 47),
                 scale = c(1.5, 6.3, 2.6), model = "gev", prior = "quant")
rn_post <- rpost_rcpp(n = 1000, model = "pp", prior = piq, data = rainfall,
                      thresh = 40, noy = 54)
plot(rn_post)

# OS model
mat <- diag(c(10000, 10000, 100))
pv <- set_prior(prior = "norm", model = "gev", mean = c(0, 0, 0), cov = mat)
osv <- rpost_rcpp(n = 1000, model = "os", prior = pv, data = venice)
plot(osv)

Prior simulation of GEV parameters - prior on probability scale

Description

Simulates from the prior distribution for GEV parameters based on Crowder (1992), in which independent beta priors are specified for ratios of probabilities (which is equivalent to a Dirichlet prior on differences between these probabilities).

Usage

rprior_prob(n, quant, alpha, exc = FALSE, lb = NULL, lb_prob = 0.001)

Arguments

Details

The simulation is based on the way that the prior is constructed. See Stephenson (1996) the evdbayes user guide or Northrop et al. (2017) Northrop et al. (2017) for details of the construction of the prior. First, differences between probabilities are simulated from a Dirichlet distribution. Then the GEV location, scale and shape parameters that correspond to these quantile values are found, by solving numerically a set of three non-linear equations in which the GEV quantile function evaluated at the simulated probabilities is equated to the quantiles in quant. This is reduced to a one-dimensional optimisation over the GEV shape parameter.

Value

An n by 3 numeric matrix.

References

Crowder, M. (1992) Bayesian priors based on parameter transformation using the distribution function. Ann. Inst. Statist. Math., 44(3), 405-416. https://link.springer.com/article/10.1007/BF00050695

Stephenson, A. 2016. Bayesian Inference for Extreme Value Modelling. In Extreme Value Modeling and Risk Analysis: Methods and Applications, edited by D. K. Dey and J. Yan, 257-80. London: Chapman and Hall. doi:10.1201/b19721

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. doi:10.1111/rssc.12159

See Also

rpost and rpost_rcpp for sampling from an extreme value posterior distribution.

Examples

quant <- c(85, 88, 95)
alpha <- c(4, 2.5, 2.25, 0.25)
x <- rprior_prob(n = 1000, quant = quant, alpha = alpha, exc = TRUE)
x <- rprior_prob(n = 1000, quant = quant, alpha = alpha, exc = TRUE, lb = 0)

Prior simulation of GEV parameters - prior on quantile scale

Description

Simulates from the prior distribution for GEV parameters proposed in Coles and Tawn (1996), based on independent gamma priors for differences between quantiles.

Usage

rprior_quant(n, prob, shape, scale, lb = NULL, lb_prob = 0.001)

Arguments

Details

The simulation is based on the way that the prior is constructed. See Coles and Tawn (1996), Stephenson (2016) or the evdbayes user guide for details of the construction of the prior. First, the quantile differences are simulated from the specified gamma distributions. Then the simulated quantiles are calculated. Then the GEV location, scale and shape parameters that give these quantile values are found, by solving numerically a set of three non-linear equations in which the GEV quantile function evaluated at the values in prob is equated to the simulated quantiles. This is reduced to a one-dimensional optimisation over the GEV shape parameter.

Value

An n by 3 numeric matrix.

References

Coles, S. G. and Tawn, J. A. (1996) A Bayesian analysis of extreme rainfall data. Appl. Statist., 45, 463-478.

Stephenson, A. 2016. Bayesian Inference for Extreme Value Modelling. In Extreme Value Modeling and Risk Analysis: Methods and Applications, edited by D. K. Dey and J. Yan, 257-80. London: Chapman and Hall. doi:10.1201/b19721

See Also

rpost and rpost_rcpp for sampling from an extreme value posterior distribution.

Examples

pr <- 10 ^ -(1:3)
sh <- c(38.9, 7.1, 47)
sc <- c(1.5, 6.3, 2.6)
x <- rprior_quant(n = 1000, prob = pr, shape = sh, scale = sc)
x <- rprior_quant(n = 1000, prob = pr, shape = sh, scale = sc, lb = 0)

Construction of a prior distribution for a binomial probability p

Description

Constructs a prior distribution for use as the argument bin_prior in rpost or in binpost. The user can choose from a list of in-built priors or specify their own prior function, returning the log of the prior density, using an R function and arguments for hyperparameters.

Usage

set_bin_prior(
  prior = c("jeffreys", "laplace", "haldane", "beta", "mdi", "northrop"),
  ...
)

Arguments

Details

Binomial priors. The names of the binomial priors set using bin_prior are:

• "jeffreys": the Jeffreys beta(1/2, 1/2) prior.

• "laplace": the Bayes-Laplace beta(1, 1) prior.

• "haldane": the Haldane beta(0, 0) prior.

• "beta": a beta(\alpha, \beta) prior. The argument ab is a vector containing c(\alpha, \beta). The default is ab = c(1, 1).

• "mdi": the MDI prior \pi(p) = 1.6186 p^p (1-p)^{1-p}, for 0 < p < 1.

• "northrop": the improper prior \pi(p)=\{-\ln(1-p)\}^{-1}(1-p)^{-1}, for 0 < p < 1.

Apart from the last two priors these are all beta distributions.

Value

A list of class "binprior". The first component is the name of the input prior. Apart from the MDI prior this will be "beta", in which case the other component of the list is a vector of length two giving the corresponding values of the beta parameters.

See Also

binpost for sampling from a binomial posterior distribution.

Examples

bp <- set_bin_prior(prior = "jeffreys")

# Setting the Jeffreys prior by hand
beta_prior_fn <- function(p, ab) {
  return(stats::dbeta(p, shape1 = ab[1], shape2 = ab[2], log = TRUE))
}
jeffreys <- set_bin_prior(beta_prior_fn, ab = c(1 / 2, 1 / 2))

Construction of prior distributions for extreme value model parameters

Description

Constructs a prior distribution for use as the argument prior in rpost and rpost_rcpp. The user can either specify their own prior function, returning the log of the prior density, (using an R function or an external pointer to a compiled C++ function) and arguments for hyperparameters or choose from a list of in-built model-specific priors. Note that the arguments model = "gev", model = "pp" and model =="os" are equivalent because a prior is specified is the GEV parameterisation in each of these cases. Note also that for model = "pp" the prior GEV parameterisation relates to the value of noy subsequently supplied to rpost or rpost_rcpp. The argument model is used for consistency with rpost.

Usage

set_prior(
  prior = c("norm", "loglognorm", "mdi", "flat", "flatflat", "jeffreys", "beta", "prob",
    "quant"),
  model = c("gev", "gp", "pp", "os"),
  ...
)

Arguments

Details

Of the in-built named priors available in revdbayes only those specified using prior = "prob" (gev_prob), prior = "quant" (gev_quant) prior = "norm" (gev_norm) or prior = "loglognorm" (gev_loglognorm) are proper. If model = "gev" these priors are equivalent to priors available in the evdbayes package, namely prior.prob, prior.quant, prior.norm and prior.loglognorm.

The other in-built priors are improper, that is, the integral of the prior function over its support is not finite. Such priors do not necessarily result in a proper posterior distribution. Northrop and Attalides (2016) consider the issue of posterior propriety in Bayesian extreme value analyses. In most of improper priors below the prior for the scale parameter \sigma is taken to be 1/\sigma, i.e. a flat prior for \log \sigma. Here we denote the scale parameter of the GP distribution by \sigma, whereas we use \sigma_u in the revdbayes vignette.

For all in-built priors the arguments min_xi and max_xi may be supplied by the user. The prior density is set to zero for any value of the shape parameter \xi that is outside (min_xi, max_xi). This will override the default values of min_xi and max_xi in the named priors detailed above.

Extreme value priors. It is typical to use either prior = "prob" (gev_prob) or prior = "quant" (gev_quant) to set an informative prior and one of the other prior (or a user-supplied function) otherwise. The names of the in-built extreme value priors set using prior and details of hyperparameters are:

• "prob". A prior for GEV parameters (\mu, \sigma, \xi) based on Crowder (1992). See gev_prob for details. See also Northrop et al. (2017) and Stephenson (2016).

• "quant". A prior for GEV parameters (\mu, \sigma, \xi) based on Coles and Tawn (1996). See gev_quant for details.

• "norm".

  For model = "gp": (\log \sigma, \xi), is bivariate normal with mean mean (a numeric vector of length 2) and covariance matrix cov (a symmetric positive definite 2 by 2 matrix).

  For model = "gev": (\mu, \log \sigma, \xi), is trivariate normal with mean mean (a numeric vector of length 3) and covariance matrix cov (a symmetric positive definite 3 by 3 matrix).

• "loglognorm". For model = "gev" only: (\log \mu, \log \sigma, \xi), is trivariate normal with mean mean (a numeric vector of length 3) and covariance matrix cov (a symmetric positive definite 3 by 3 matrix).

• "mdi".

  For model = "gp": (an extended version of) the maximal data information (MDI) prior, that is,

  \pi(\sigma, \xi) = \sigma^{-1} \exp[-a(\xi + 1)], {\rm ~for~} \sigma > 0, \xi \geq -1, a \geq 0.

  The value of a is set using the argument a. The default value is a = 1, which gives the MDI prior.

  For model = "gev": (an extended version of) the maximal data information (MDI) prior, that is,

  \pi(\mu, \sigma, \xi) = \sigma^{-1} \exp[-a(\xi + 1)], {\rm ~for~} \sigma > 0, \xi \geq -1, a \geq 0.

  The value of a is set using the argument a. The default value is a = \gamma, where \gamma = 0.57721 is Euler's constant, which gives the MDI prior.

  For each of these cases \xi must be is bounded below a priori for the posterior to be proper (Northrop and Attalides, 2016). An argument for the bound \xi \geq -1 is that for \xi < -1 the GP (GEV) likelihood is unbounded above as -\sigma/\xi (\mu - \sigma/\xi)) approaches the sample maximum. In maximum likelihood estimation of GP parameters (Grimshaw, 1993) and GEV parameters a local maximum of the likelihood is sought on the region \sigma > 0, \xi \geq -1.

• "flat".

  For model = "gp": a flat prior for \xi (and for \log \sigma):

  \pi(\sigma, \xi) = \sigma^{-1}, {\rm ~for~} \sigma > 0.

  For model = "gev": a flat prior for \xi (and for \mu and \log \sigma):

  \pi(\mu, \sigma, \xi) = \sigma^{-1}, {\rm ~for~} \sigma > 0.

• "flatflat".

  For model = "gp": flat priors for \sigma and \xi:

  \pi(\sigma, \xi) = 1, {\rm ~for~} \sigma > 0.

  For model = "gev": flat priors for \mu, \sigma and \xi:

  \pi(\mu, \sigma, \xi) = 1, {\rm ~for~} \sigma > 0.

  Therefore, the posterior is proportional to the likelihood.

• "jeffreys". For model = "gp" only: the Jeffreys prior (Castellanos and Cabras, 2007):

  \pi(\sigma, \xi) = \sigma^{-1}(1+\xi)^{-1}(1+2\xi)^{-1/2}, {\rm ~for~} \sigma > 0, \xi > -1 / 2.

  In the GEV case the Jeffreys prior doesn't yield a proper posterior for any sample size. See Northrop and Attalides (2016) for details.

• "beta". For model = "gp": a beta-type(p, q) prior is used for xi on the interval (min_xi, max_xi):

  \pi(\sigma, \xi) = \sigma^{-1} (\xi - {\min}_{\xi}) ^ {p-1} ({\max}_{\xi} - \xi) ^ {q-1}, {\rm ~for~} {\min}_{\xi} < \xi < {\max}_{\xi}.

  For model = "gev": similarly ...

  \pi(\mu, \sigma, \xi) = \sigma^{-1} (\xi - {\min}_{\xi}) ^ {p-1} ({\max}_{\xi} - \xi) ^ {q-1}, {\rm ~for~} {\min}_{\xi} < \xi < {\max}_{\xi}.

  The argument pq is a vector containing c(p,q). The default settings for this prior are p = 6, q = 9 and min_xi = -1/2, max_xi = 1/2, which corresponds to the prior for \xi proposed in Martins and Stedinger (2000, 2001).

Value

A list with class "evprior". The first component is the input prior, i.e. either the name of the prior or a user-supplied function. The remaining components contain the numeric values of any hyperparameters in the prior.

References

Castellanos, E. M. and Cabras, S. (2007) A default Bayesian procedure for the generalized Pareto distribution. Journal of Statistical Planning and Inference 137(2), 473-483. doi:10.1016/j.jspi.2006.01.006.

Coles, S. G. and Tawn, J. A. (1996) A Bayesian analysis of extreme rainfall data. Appl. Statist., 45, 463-478.

Crowder, M. (1992) Bayesian priors based on parameter transformation using the distribution function Ann. Inst. Statist. Math., 44, 405-416. https://link.springer.com/article/10.1007/BF00050695.

Grimshaw, S. D. (1993) Computing Maximum Likelihood Estimates for the Generalized Pareto Distribution. Technometrics, 35(2), 185-191. doi:10.1080/00401706.1993.10485040.

Hosking, J. R. M. and Wallis, J. R. (1987) Parameter and Quantile Estimation for the Generalized Pareto Distribution. Technometrics, 29(3), 339-349. doi:10.2307/1269343.

Martins, E. S. and Stedinger, J. R. (2000) Generalized maximum likelihood generalized extreme value quantile estimators for hydrologic data, Water Resources Research, 36(3), 737-744. doi:10.1029/1999WR900330.

Martins, E. S. and Stedinger, J. R. (2001) Generalized maximum likelihood Pareto-Poisson estimators for partial duration series, Water Resources Research, 37(10), 2551-2557. doi:10.1029/2001WR000367.

Northrop, P.J. and Attalides, N. (2016) Posterior propriety in Bayesian extreme value analyses using reference priors Statistica Sinica, 26(2), 721–743 doi: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. doi:10.1111/rssc.12159

Stephenson, A. (2016) Bayesian inference for extreme value modelling. In Extreme Value Modeling and Risk Analysis: Methods and Applications (eds D. K. Dey and J. Yan), 257-280, Chapman and Hall, London. doi:10.1201/b19721.

See Also

rpost and rpost_rcpp for sampling from an extreme value posterior distribution.

create_prior_xptr for creating an external pointer to a C++ function to evaluate the log-prior density.

rprior_prob and rprior_quant for sampling from informative prior distributions for GEV parameters.

gp_norm, gp_mdi, gp_flat, gp_flatflat, gp_jeffreys, gp_beta to see the arguments for priors for GP parameters.

gev_norm, gev_loglognorm, gev_mdi, gev_flat, gev_flatflat, gev_beta, gev_prob, gev_quant to see the arguments for priors for GEV parameters.

Examples

# Normal prior for GEV parameters (mu, log(sigma), xi).
mat <- diag(c(10000, 10000, 100))
pn <- set_prior(prior = "norm", model = "gev", mean = c(0,0,0), cov = mat)
pn

# Prior for GP parameters with flat prior for xi on (-1, infinity).
fp <- set_prior(prior = "flat", model = "gp", min_xi = -1)
fp

# A user-defined prior (see the vignette for details).
u_prior_fn <- function(x, ab) {
  if (x[1] <= 0 | x[2] <= -1 | x[2] >= 1) {
    return(-Inf)
  }
  return(-log(x[1]) + (ab[1] - 1) * log(1 + x[2]) +
         (ab[2] - 1) * log(1 - x[2]))
}
up <- set_prior(prior = u_prior_fn, ab = c(2, 2), model = "gp")

# A user-defined prior using a pointer to a C++ function
ptr_gp_flat <- create_prior_xptr("gp_flat")
u_prior_ptr <- set_prior(prior = ptr_gp_flat, model = "gp")

Summarizing an evpost object

Description

summary method for class "evpost"

Usage

## S3 method for class 'evpost'
summary(object, add_pu = FALSE, ...)

Arguments

Value

Prints

• information about the ratio-of-uniforms bounding box, i.e. object$box

• an estimate of the probability of acceptance, i.e. object$pa

• a summary of the simulated values, via summary(object$sim_vals)

See Also

ru or ru_rcpp for descriptions of object$sim_vals and object$box.

plot.evpost for a diagnostic plot.

Examples

# GP posterior
u <- stats::quantile(gom, probs = 0.65)
fp <- set_prior(prior = "flat", model = "gp", min_xi = -1)
gpg <- rpost_rcpp(n = 1000, model = "gp", prior = fp, thresh = u,
                  data = gom)
summary(gpg)

Largest Sea Levels in Venice

Description

The venice data frame has 51 rows and 10 columns. The jth column contains the jth largest sea levels in Venice, for the years 1931-1981. Only the largest six measurements are available for the year 1935; the corresponding row contains four missing values. The years for each set of measurements are given as row names.

Usage

venice

Format

A data frame with 51 rows and 10 columns.

Source

Smith, R. L. (1986) Extreme value theory based on the r largest annual events. Journal of Hydrology, 86, 27-43. doi:10.1016/0022-1694(86)90004-1

References

Coles, S. G. (2001) An Introduction to Statistical Modelling of Extreme Values. London: Springer. doi:10.1007/978-1-4471-3675-0

Random sampling from a binomial posterior distribution, using weights

Description

Samples from the posterior distribution of the probability p of a binomial distribution. User-supplied weights are applied to each observation when constructing the log-likelihood.

Usage

wbinpost(n, prior, ds_bin)

Arguments

Details

For prior$prior == "bin_beta" the posterior for p is a beta distribution so rbeta is used to sample from the posterior.

Value

An object (list) of class "binpost" with components

See Also

set_bin_prior for setting a prior distribution for the binomial probability p.

Examples

u <- quantile(gom, probs = 0.65)
ds_bin <- list(sf = gom > u, w = rep(1, length(gom)))
bp <- set_bin_prior(prior = "jeffreys")
temp <- wbinpost(n = 1000, prior = bp, ds_bin = ds_bin)
graphics::hist(temp$bin_sim_vals, prob = TRUE)