rmev {mev}

Exact simulations of multivariate extreme value distributions

Description

Implementation of the random number generators for multivariate extreme-value distributions and max-stable processes based on the two algorithms described in Dombry, Engelke and Oesting (2016).

Usage

rmev(
  n,
  d,
  param,
  asy,
  sigma,
  model = c("log", "alog", "neglog", "aneglog", "bilog", "negbilog", "hr", "br", "xstud",
    "smith", "schlather", "ct", "sdir", "dirmix", "pairbeta", "pairexp", "wdirbs",
    "wexpbs", "maxlin"),
  alg = c("ef", "sm"),
  weights = NULL,
  vario = NULL,
  coord = NULL,
  grid = FALSE,
  dist = NULL,
  ...
)

Arguments

n

number of observations

d

dimension of sample

param

parameter vector for the logistic, bilogistic, negative bilogistic and extremal Dirichlet (Coles and Tawn) model. Parameter matrix for the Dirichlet mixture. Degree of freedoms for extremal student model. See Details.

asy

list of asymmetry parameters, as in function rmvevd from package evd, of 2^d-1 vectors of size corresponding to the power set of d, with sum to one constraints.

sigma

covariance matrix for Brown-Resnick and extremal Student-t distributions. Symmetric matrix of squared coefficients \lambda^2 for the Husler-Reiss model, with zero diagonal elements.

model

for multivariate extreme value distributions, users can choose between 1-parameter logistic and negative logistic, asymmetric logistic and negative logistic, bilogistic, Husler-Reiss, extremal Dirichlet model (Coles and Tawn) or the Dirichlet mixture. Spatial models include the Brown-Resnick, Smith, Schlather and extremal Student max-stable processes. Max linear models are also supported

alg

algorithm, either simulation via extremal function ('ef') or via the spectral measure ('sm'). Default to ef.

weights

vector of length m for the m mixture components that sum to one. For the "maxlin" model, weights should be a matrix with d columns that represent the weight of the components and whose column sum to one (if provided, this argument overrides asy).

vario

semivariogram function whose first argument must be distance. Used only if provided in conjunction with coord and if sigma is missing

coord

d by k matrix of coordinates, used as input in the variogram vario or as parameter for the Smith model. If grid is TRUE, unique entries should be supplied.

grid

Logical. TRUE if the coordinates are two-dimensional grid points (spatial models).

dist

symmetric matrix of pairwise distances. Default to NULL.

...

additional arguments for the vario function

Details

The vector param differs depending on the model

• log: one dimensional parameter greater than 1

• alog: 2^d-d-1 dimensional parameter for dep. Values are recycled if needed.

• neglog: one dimensional positive parameter

• aneglog: 2^d-d-1 dimensional parameter for dep. Values are recycled if needed.

• bilog: d-dimensional vector of parameters in [0,1]

• negbilog: d-dimensional vector of negative parameters

• ct, dir, negdir, sdir: d-dimensional vector of positive (a)symmetry parameters. For dir and negdir, a d+1 vector consisting of the d Dirichlet parameters and the last entry is an index of regular variation in (-\min(\alpha_1, \ldots, \alpha_d), 1] treated as shape parameter

• xstud: one dimensional parameter corresponding to degrees of freedom alpha

• dirmix: d by m-dimensional matrix of positive (a)symmetry parameters

• pairbeta, pairexp: d(d-1)/2+1 vector of parameters, containing the concentration parameter and the coefficients of the pairwise beta, in lexicographical order e.g., \beta_{12}, \beta_{13}, \ldots

• wdirbs, wexpbs: 2d vector of d concentration parameters followed by the d Dirichlet parameters

Stephenson points out that the multivariate asymmetric negative logistic model given in e.g. Coles and Tawn (1991) is not a valid distribution function in dimension d>3 unless additional constraints are imposed on the parameter values. The implementation in mev uses the same construction as the asymmetric logistic distribution (see the vignette). As such it does not match the bivariate implementation of rbvevd.

The dependence parameter of the evd package for the Husler-Reiss distribution can be recovered taking for the Brown–Resnick model 2/r=\sqrt(2\gamma(h)) where h is the lag vector between sites and r=1/\lambda for the Husler–Reiss.

Value

an n by d exact sample from the corresponding multivariate extreme value model

Warning

As of version 1.8 (August 16, 2016), there is a distinction between models hr and br. The latter is meant to be used in conjunction with variograms. The parametrization differs between the two models.

The family of scaled Dirichlet is now parametrized by a parameter in -\min(\alpha) appended to the the d vector param containing the parameter alpha of the Dirichlet model. Arguments model='dir' and model='negdir' are still supported internally, but not listed in the options.

Author(s)

Leo Belzile

References

Dombry, Engelke and Oesting (2016). Exact simulation of max-stable processes, Biometrika, 103(2), 303–317.

See Also

rmevspec, rmvevd, rbvevd

Examples

set.seed(1)
rmev(n=100, d=3, param=2.5, model='log', alg='ef')
rmev(n=100, d=4, param=c(0.2,0.1,0.9,0.5), model='bilog', alg='sm')
## Spatial example using power variogram
#NEW: Semi-variogram must take distance as argument
semivario <- function(x, scale, alpha){ scale*x^alpha }
#grid specification
grid.coord <- as.matrix(expand.grid(runif(4), runif(4)))
rmev(n=100, vario=semivario, coord=grid.coord, model='br', scale = 0.5, alpha = 1)
#using the Brown-Resnick model with a covariance matrix
vario2cov <- function(coord, semivario,...){
 sapply(1:nrow(coord), function(i) sapply(1:nrow(coord), function(j)
  semivario(sqrt(sum((coord[i,])^2)), ...) +
  semivario(sqrt(sum((coord[j,])^2)), ...) -
  semivario(sqrt(sum((coord[i,]-coord[j,])^2)), ...)))
}
rmev(n=100, sigma=vario2cov(grid.coord, semivario = semivario, scale = 0.5, alpha = 1), model='br')
# asymmetric logistic model - see function 'rmvevd' from package 'evd '
asy <- list(0, 0, 0, 0, c(0,0), c(0,0), c(0,0), c(0,0), c(0,0), c(0,0),
  c(.2,.1,.2), c(.1,.1,.2), c(.3,.4,.1), c(.2,.2,.2), c(.4,.6,.2,.5))
rmev(n=1, d=4, param=0.3, asy=asy, model="alog")
#Example with a grid (generating an array)
rmev(n=10, sigma=cbind(c(2,1), c(1,3)), coord=cbind(runif(4), runif(4)), model='smith', grid=TRUE)
## Example with Dirichlet mixture
alpha.mat <- cbind(c(2,1,1),c(1,2,1),c(1,1,2))
rmev(n=100, param=alpha.mat, weights=rep(1/3,3), model='dirmix')
rmev(n=10, param=c(0.1,1,2,3), d=3, model='pairbeta')