Generalised Pareto (GP) model
The GP distribution is used to model threshold excesses,
that is, the amounts by which threshold exceedances exceed some
threshold. It has two parameters: a scale parameter \(\sigma_u\) and shape parameter \(\xi\). We use the following example to
illustrate the use of the transformation strategies implemented by the
arguments rotate and trans to
rpost and how the user can specify their own prior
distribution.
Gulf of Mexico Wave Height Data
The numeric vector gom contains 315 storm peak
significant wave heights from a location in the Gulf of Mexico, from the
years 1900 to 2005. These data are analysed in Northrop, Attalides, and Jonathan (2017). We set
the threshold at the 65% sample quantile and use set_prior
to set a prior.
We sample first on the \((\sigma_u, \xi)\) scale, with mode relocation only.
t1 <- system.time(
gp1 <- rpost(n = n, model = "gp", prior = fp, thresh = thresh, data = gom,
rotate = FALSE)
)[3]Then we add a rotation of \((\sigma_u, \xi)\) about the estimated posterior mode.
Now we add marginal Box-Cox transformation. We apply Box-Cox transformation to the parameters \(\phi_1 = \sigma_u\) and \(\phi_2 = \xi + \sigma / x_{(m)}\), where \(x_{(m)}\) is the largest threshold excess. The parameters \(\phi_1\) and \(\phi_2\) are positive for all combinations of \((\sigma_u, \xi)\) for which the GP likelihood is positive.
t3 <- system.time(
gp3 <- rpost(n = n, model = "gp", prior = fp, thresh = thresh, data = gom,
rotate = FALSE, trans = "BC")
)[3]
t4 <- system.time(
gp4 <- rpost(n = n, model = "gp", prior = fp, thresh = thresh, data = gom,
trans = "BC")
)[3]We plot the samples obtained with the contours of the corresponding densities superimposed. The plot on the top left is on the original \((\sigma_u, \xi)\) scale. The other plots are on the scale used for the ratio-of-uniforms algorithm, that is, with relocation of the mode to the origin. In the following the \(\rho_i, i = 1, \ldots, d\), \(\rho_1\) and \(\rho_2\) in this example, are the variables to which the ratio-of-uniforms algorithm is applied, i.e. after any transformation (Box-Cox and/or rotation of axes) and relocation of the mode to the origin.
plot(gp1, ru_scale = FALSE, cex.main = 0.75, cex.lab = 0.75,
main = paste("no transformation \n pa = ", round(gp1$pa, 3),
", time = ", round(t1, 2), "s"))
plot(gp2, ru_scale = TRUE, cex.main = 0.75, cex.lab = 0.75,
main = paste("rotation \n pa = ", round(gp2$pa, 3),
", time = ", round(t2, 2), "s"))
plot(gp3, ru_scale = TRUE, cex.main = 0.75, cex.lab = 0.75,
main = paste("Box-Cox \n pa = ", round(gp3$pa, 3),
", time = ", round(t3, 2), "s"))
plot(gp4, ru_scale = TRUE, cex.main = 0.75, cex.lab = 0.75,
main = paste("Box-Cox and rotation \n pa = ", round(gp4$pa, 3),
", time = ", round(t4, 2), "s"))
Comparison of the plots on the right to the plots on the left shows
that rotation of the parameter axes about the mode of the posterior has
reduced dependence between the components. The estimated probabilities
of acceptance for the plots on the right are close to the 0.534 obtained
for a 2-dimensional normal distribution with independent components.
Box-Cox transformation has increased the estimated value of \(p_a\) but not by much. In fact, the extra
computation time required to calculate the Box-Cox transformation each
time that the posterior density is evaluated means that posterior
sampling is slower than the default setting of
rotate = TRUE and trans = "none". That is, in
this example at least, the increase is \(p_a\) resulting from the addition of the
Box-Cox transformation is not sufficient to offset the extra time needed
to compute the posterior density.
We repeat the rotation only and Box-Cox transformation plus rotation analyses for a much higher threshold, set at the 95% sample quantile.
thresh <- quantile(gom, probs = 0.95)
fp <- set_prior(prior = "flat", model = "gp", min_xi = -1)
t2 <- system.time(
gp2 <- rpost(n = n, model = "gp", prior = fp, thresh = thresh, data = gom)
)[3]
t4 <- system.time(
gp4 <- rpost(n = n, model = "gp", prior = fp, thresh = thresh, data = gom,
trans = "BC")
)[3]It is optimistic to use such a high threshold in this example, because it results in only 16 threshold excesses. In the rotation only case, a convergence warning is triggered by one of the optimisations used to created the bounding box of the ratio-of-uniforms method. This appears to be a spurious warning, because the plots of the simulations below seem fine. When a Box-Cox transformation is used there are no warnings, because the optimisations have greater stability and it is now easier for the convergence of the optimisations to be confirmed.
Strong asymmetry in the posterior distribution means that the combination of marginal Box-Cox transformation and rotation produces a larger improvement in \(p_a\) and means that now this strategy is competitive with using rotation alone in terms of computation time.
A User-defined prior
Suppose that we wish to use a prior for \((\sigma_u, \xi)\) with a density that is
proportional to \(\sigma_u^{-1}
(1+\xi)^{\alpha-1} (1-\xi)^{\beta-1}\), for \(\sigma_u > 0, -1 < \xi < 1\) and
for some \(\alpha > 0\) and \(\beta > 0\). This is an improper prior
in which \(\sigma_u\) and \(\xi\) are independent a priori,
\(\log \sigma_u\) is uniform over the
real line and \(\xi\) has beta(\(\alpha, \beta\))-type distribution on the
interval \((-1, 1)\). We can do this by
creating a function that returns the prior log-density and passing this
function to set_prior. The first argument of the prior
log-density function must be the parameter vector (the GP parameters
\((\sigma_u, \xi)\) in this case),
followed by any hyperparameters.
u_prior_fn <- function(x, ab) {
#
# Calculates the the log of the (improper) prior density for GP parameters
# (sigma_u, xi) in which log(sigma_u) is uniform on the real line and xi has
# a beta(alpha, beta)-type prior on the interval (-1, 1).
#
# Args:
# x : A numeric vector. GP parameter vector (sigma, xi).
# ab : A numeric vector. Hyperparameter vector (alpha, beta), where
# alpha and beta must be positive.
#
# Returns : the value of the log-prior at x.
#
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")
gp_u <- rpost(n = n, model = "gp", prior = up, thresh = thresh, data = gom)See also the “Posterior predictive extreme value inference” section of the Posterior Predictive Extreme Value Inference using the revdbayes Package for an analysis that adds inferences about the probability of threshold exceedance to the analysis of threshold excesses using a GP distribution.
