predict.evpost | R Documentation |
N
years.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.
## 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,
...
)
object |
An object of class
|
type |
A character vector. Indicates which type of inference is required:
|
x |
A numeric vector or a matrix with
|
x_num |
A numeric scalar. If |
n_years |
A numeric vector. Values of |
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' worth of non-missing data. If Otherwise, a default value will be assumed if
If |
level |
A numeric vector of values in (0, 100).
Only relevant when |
hpd |
A logical scalar.
Only relevant when If If |
lower_tail |
A logical scalar.
Only relevant when |
log |
A logical scalar. Only relevant when |
big_q |
A numeric scalar. Only relevant when |
... |
Additional optional arguments. At present no optional arguments are used. |
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 NA
s 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 NA
s.
An object of class "evpred", a list containing a subset of the following components:
type |
The argument |
x |
A matrix containing the argument |
y |
The content of
|
long |
A |
short |
A matrix with the same structure as |
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
.
Coles, S. G. (2001) An Introduction to Statistical Modeling of Extreme Values, Springer-Verlag, London. Chapter 9: \Sexpr[results=rd]{tools:::Rd_expr_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. \Sexpr[results=rd]{tools:::Rd_expr_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. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1201/b19721")}
plot.evpred
for the S3 plot
method for
objects of class evpred
.
rpost
or rpost_rcpp
for sampling
from an extreme value posterior distribution.
### 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"))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.