hkevp.predict {hkevp}

Predictive distribution of the max-stable process at target positions.

Description

Computes the predictive distribution of Y(\cdot) at a set of ungauged positions (s_1^*, ..., s_k^*), given data at gauged positions (s_1, ..., s_n), by using the output of latent.fit or hkevp.fit.

Two types of prediction are available for the HKEVP, as described in Shaby and Reich (2012). See details.

Usage

hkevp.predict(fit, targets, targets.covariates, predict.type = "kriging")

Arguments

fit

Output from the hkevp.fit procedure.

targets

A matrix of real values giving the spatial coordinates of the ungauged positions. Each row corresponds to an ungauged position.

targets.covariates

A matrix of real values giving the spatial covariates of the ungauged positions. Must match with the covariates used in hkevp.fit or latent.fit.

predict.type

Character string specifying the type of prediction. Must be one of "kriging" (default) or "climat". See details.

Details

The spatial prediction of Y_t(s^*) for a target site s^* and a realisation t of the process is described in Shaby and Reich (2012). This method involves a three-step procedure:

1. Computation of the residual dependence process \theta(\cdot) at the target positions.

2. Computation of the conditional GEV parameters (\mu^*,\sigma^*,\xi^*) at the target sites. See the definition of the HKEVP in Reich and Shaby (2012).

3. Generation of Y_t(s^*) from an independent GEV distribution with parameters (\mu^*,\sigma^*,\xi^*).

As sketched in Shaby and Reich (2012), two types of prediction are possible: the kriging-type and the climatological-type. These two types differ when the residual dependence process \theta is computed (first step of the prediction):

• The kriging-type takes the actual value of A in the MCMC algorithm to compute the residual dependence process. The prediction will be the distribution of the maximum recorded at the specified targets.

• The climatological-type generates A by sampling from the positive stable distribution with characteristic exponent \alpha, where \alpha is the actual value of the MCMC step. The prediction in climatological-type will be the distribution of what could happen in the conditions of the HKEVP dependence structure.

Posterior distribution for each realisation t of the process and each target position s^* is represented with a sample where each element corresponds to a step of the MCMC procedure.

Value

A three-dimensional array where:

• Each row corresponds to a different realisation of the process (a block).

• Each column corresponds to a target position.

• Each slice corresponds to a MCMC step.

Author(s)

Quentin Sebille

References

Reich, B. J., & Shaby, B. A. (2012). A hierarchical max-stable spatial model for extreme precipitation. The annals of applied statistics, 6(4), 1430. <DOI:10.1214/12-AOAS591>

Shaby, B. A., & Reich, B. J. (2012). Bayesian spatial extreme value analysis to assess the changing risk of concurrent high temperatures across large portions of European cropland. Environmetrics, 23(8), 638-648. <DOI:10.1002/env.2178>

Examples

# Simulation of HKEVP:
sites <- as.matrix(expand.grid(1:3,1:3))
targets <- as.matrix(expand.grid(1.5:2.5,1.5:2.5))
all.pos <- rbind(sites, targets)
knots <- sites
loc <- all.pos[,1]*10
scale <- 3
shape <- 0
alpha <- .4
tau <- 1
ysim <- hkevp.rand(10, all.pos, knots, loc, scale, shape, alpha, tau)
yobs <- ysim[,1:9]

# HKEVP fit (omitting first site, used as target):
fit <- hkevp.fit(yobs, sites, niter = 1000)

# Extrapolation:
ypred <- hkevp.predict(fit, targets, predict.type = "kriging")

# Plot of the density and the true value for 4 first realizations:
# par(mfrow = c(2, 2))
# plot(density(ypred[1,1,]), main = "Target 1 / Year 1")
# abline(v = ysim[1,10], col = 2, lwd = 2)
# plot(density(ypred[2,1,]), main = "Target 1 / Year 2")
# abline(v = ysim[2,10], col = 2, lwd = 2)
# plot(density(ypred[1,2,]), main = "Target 2 / Year 1")
# abline(v = ysim[1,11], col = 2, lwd = 2)
# plot(density(ypred[2,2,]), main = "Target 2 / Year 2")
# abline(v = ysim[2,11], col = 2, lwd = 2)