predict method for class "blite". 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 blite
predict(
  object,
  type = c("i", "p", "d", "q", "r"),
  x = NULL,
  x_num = 100,
  n_years = 100,
  ny = NULL,
  level = 95,
  hpd = FALSE,
  lower_tail = TRUE,
  log = FALSE,
  big_q = 1000,
  ...
)

Arguments

object

An object of class "blite" returned from blite.

type

A character vector. Indicates which type of inference is required:

  • "i" for predictive intervals,

  • "p" for the predictive distribution function,

  • "d" for the predictive density function,

  • "q" for the predictive quantile function,

  • "r" for random generation from the predictive distribution.

x

A numeric vector or a matrix with n_years columns. The meaning of x depends on type.

  • type = "p" or type = "d": x contains quantiles at which to evaluate the distribution or density function. No element of x can be less than the threshold attr(object, "inputs")$u.

    If x is not supplied then n_year-specific defaults are set: vectors of length x_num from the 0.1% quantile to the 99% quantile, subject all values being greater than the threshold.

  • type = "q": x contains probabilities in (0,1) at which to evaluate the quantile function. Any values outside (0, 1) will be removed without warning. No element of p can correspond to a predictive quantile that is below the threshold, attr(object, "inputs")$u. That is, no element of p can be less than the value of predict.evpost(object, type = "q", x = attr(object, "inputs")$u).

    If x is not supplied then a default value of c(0.025, 0.25, 0.5, 0.75, 0.975) is used.

  • type = "i" or type = "r": x is not relevant.

x_num

A numeric scalar. If type = "p" or type = "d" and x is not supplied then x_num gives the number of values in x for each value in n_years.

n_years

A numeric vector. Values of \(N\).

ny

A numeric scalar. The (mean) number of observations per year. Setting this appropriately is important. See Details.

level

A numeric vector of values in (0, 100). Only relevant when type = "i". Levels of predictive intervals for the largest value observed in \(N\) years, i.e. level% predictive intervals are returned.

hpd

A logical scalar. Only relevant when type = "i".

If hpd = FALSE then the interval is equi-tailed, with its limits produced by
predict.evpost(object, type ="q", x = p), where p = c((1-level/100)/2, (1+level/100)/2).

If hpd = TRUE then, in addition to the equi-tailed interval, the shortest possible level% interval is calculated. If the predictive distribution is unimodal then this is a highest predictive density (HPD) interval.

lower_tail

A logical scalar. Only relevant when type = "p" or type = "q". If TRUE (default), (output or input) probabilities are \(P[X \leq x]\), otherwise, \(P[X > x]\).

log

A logical scalar. Only relevant when type = "d". If TRUE the log-density is returned.

big_q

A numeric scalar. Only relevant when type = "q". An initial upper bound for the desired quantiles to be passed to uniroot (its argument upper) in the search for the predictive quantiles. If this is not sufficiently large then it is increased until it does provide an upper bound.

...

Additional optional arguments. At present no optional arguments are used.

Value

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

type

The argument type supplied to predict.blite. Which of the following components are present depends type.

x

A matrix containing the argument x supplied to predict.blite, or set within predict.blite if x was not supplied, replicated to have n_years columns if necessary. Only present if type is "p", "d" or "q".

y

The content of y depends on type:

  • type = "p", "d", "q": A matrix with the same dimensions as x. Contains distribution function values (type = "p"), predictive density (type = "d") or quantiles (type = "q").

  • type = "r": A numeric matrix with length(n_years) columns and number of rows equal to the size of the posterior sample.

  • type = "i": y is not present.

long

A length(n_years)*length(level) by 4 numeric matrix containing the equi-tailed limits with columns: lower limit, upper limit, n_years, level. Only present if type = "i". If an interval extends below the threshold then NA is returned.

short

A matrix with the same structure as long containing the HPD limits. Only present if type = "i". Columns 1 and 2 contain NAs if hpd = FALSE or if the corresponding equi-tailed interval extends below the threshold.

The arguments n_years, level, hpd, lower_tail, log supplied to predict.blite are also included, as is the value of ny

and model = "bingp".

Details

The function predict.evpost in the revdbayes package is used to perform the predictive inferences. The effect of adjusting for the values of the extremal index \(\theta\) in the posterior sample in object$sim_vals[, "theta"] is to change the effective time horizon from \(N\) to \(\theta N\).

ny provides information about the (intended) frequency of sampling in time, that is, the number of observations that would be observed in a year if there are no missing values. If the number of observations may vary between years then ny should be set equal to the mean number of observations per year.

Supplying ny. The value of ny may have been set in the call to blite. If ny is supplied by the user in the call to predict.blite then this will be used in preference to the value stored in the fitted model object. If these two values differ then no warning will be given.

Examples

### Cheeseboro wind gusts

cdata <- exdex::cheeseboro
# Each column of the matrix cdata corresponds to data from a different year
# blite() sets cluster automatically to correspond to column (year)
cpost <- blite(cdata, u = 45, k = 3, ny = 31 * 24)

# Interval estimation
predict(cpost)$long
#>         lower    upper n_years level
#> [1,] 76.81578 153.4328     100    95
predict(cpost, hpd = TRUE)$short
#>         lower   upper n_years level
#> [1,] 72.52207 138.269     100    95

# Density function
plot(predict(cpost, type = "d", n_years = c(100, 1000)))


# Distribution function
plot(predict(cpost, type = "p", n_years = c(100, 1000)))


# Quantiles
predict(cpost, type = "q", n_years = c(100, 1000))$y
#>           [,1]      [,2]
#> [1,]  76.81578  85.00125
#> [2,]  86.93571  95.95879
#> [3,]  94.46424 106.13034
#> [4,] 106.00296 124.03650
#> [5,] 153.43282 203.70407

# Random generation
plot(predict(cpost, type = "r"))