graphicalExtremes: Statistical methodology for graphical extreme value models.

Description

An implementation of the statistical methodology paper Engelke and Hitz (2020) for sparse multivariate extreme value models. Includes exact simulation algorithms and statistical inference methods for multivariate Pareto distributions on graphical structures. Also contains implementations of statistical methods from Engelke and Volgushev (2022), Röttger et al. (2021), and Hentschel et al. (2022).

Details

The following global options are used by functions in the package. Their values can be changed using base::options().

"graphicalExtremes.mc.cores"

The (maximal) number of cores to use in parallel tasks. Will always be overwritten by 1 on Windows.

"graphicalExtremes.tol.small"

The "small" tolerance is used in internal computations for values that should mathematically be exactly equal to zero, but deviate due to inherent limitations of numerical computations. This value is used e.g. when checking matrices for symmetry and definiteness. In general, this value is used only as a "permissive" tolerance, in the sense that if a value has to be positive, it is compared to actual zero, but if it has to be zero, its absolute value is compared to this tolerance.

"graphicalExtremes.tol.large"

The "large" tolerance is used for values that converge to zero, but are mathematically not supposed to be equal to zero. This value is used e.g. when converting a precision matrix \Theta to an adjacency matrix of a graph.

⁠"graphicalExtremes.default.alert⁠

The default alert function to be used in validity checks of Huesler-Reiss parameter matrix transformations. Can be a function that takes an arbitrary number of strings as arguments (e.g. cat(), stop()), FALSE to ignore the alerts, or TRUE/NULL to use the default function warning().

Author(s)

Maintainer: Sebastian Engelke sebastian.engelke@unige.ch

Authors:

• Adrien S. Hitz

• Nicola Gnecco

• Manuel Hentschel

References

Engelke S, Hitz AS (2020). “Graphical models for extremes (with discussion).” J. R. Stat. Soc. Ser. B Stat. Methodol., 82, 871–932.

Engelke S, Volgushev S (2022). “Structure learning for extremal tree models.” J. R. Stat. Soc. Ser. B Stat. Methodol.. doi:10.1111/rssb.12556, Forthcoming, https://rss.onlinelibrary.wiley.com/doi/pdf/10.1111/rssb.12556.

Hentschel M, Engelke S, Segers J (2022). “Statistical Inference for Hüsler-Reiss Graphical Models Through Matrix Completions.” doi:10.48550/ARXIV.2210.14292, https://arxiv.org/abs/2210.14292.

Röttger F, Engelke S, Zwiernik P (2021). “Total positivity in multivariate extremes.” doi:10.48550/ARXIV.2112.14727, https://arxiv.org/abs/2112.14727.

See Also

Useful links:

• https://github.com/sebastian-engelke/graphicalExtremes

• Report bugs at https://github.com/sebastian-engelke/graphicalExtremes/issues

Conversion between Huesler-Reiss parameter matrices

Description

Converts between different matrices that parametrize the same Huesler-Reiss distribution: \Gamma, \Sigma, \Theta, \Sigma^k, \Theta^k. The (d-1) \times (d-1) matrices \Sigma^k and \Theta^k can also be given/returned as d \times d matrices with the kth row and column filled with zeros.

Usage

Gamma2Sigma(Gamma, k = NULL, full = FALSE, check = TRUE)

Gamma2Theta(Gamma, k = NULL, full = FALSE, check = TRUE)

Sigma2Gamma(Sigma, k = NULL, full = FALSE, check = TRUE)

Theta2Gamma(Theta, k = NULL, full = FALSE, check = TRUE)

Sigma2Theta(
  Sigma,
  k1 = NULL,
  k2 = NULL,
  full1 = FALSE,
  full2 = FALSE,
  check = TRUE
)

Theta2Sigma(
  Theta,
  k1 = NULL,
  k2 = NULL,
  full1 = FALSE,
  full2 = FALSE,
  check = TRUE
)

Theta2Theta(
  Theta,
  k1 = NULL,
  k2 = NULL,
  full1 = FALSE,
  full2 = FALSE,
  check = TRUE
)

Sigma2Sigma(
  Sigma,
  k1 = NULL,
  k2 = NULL,
  full1 = FALSE,
  full2 = FALSE,
  check = TRUE
)

Gamma2Gamma(Gamma, check = TRUE)

matrix2matrix(
  M,
  name1 = c("Gamma", "Sigma", "Theta")[1],
  name2 = c("Gamma", "Sigma", "Theta")[1],
  k1 = NULL,
  k2 = NULL,
  full1 = FALSE,
  full2 = FALSE,
  check = TRUE
)

Arguments

Details

If k, k1, or k2 is NULL, the corresponding ⁠full*⁠ argument is ignored.

Gamma2Gamma only checks and returns the input.

matrix2matrix is a wrapper function that calls the corresponding conversion function implied by name1, name2.

Value

The desired parameter matrix corresponding to the specified inputs.

See Also

Other parameter matrix transformations: Gamma2graph(), chi2Gamma(), par2Matrix()

Compute theoretical \chi in 3D

Description

Computes the theoretical \chi coefficient in 3 dimensions.

Usage

Gamma2chi_3D(Gamma)

Arguments

Value

The 3-dimensional \chi coefficient, i.e., the extremal correlation coefficient for the HR distribution. Note that 0 \leq \chi \leq 1.

Convert matrix to graph

Description

Creates a graph object representing the graph structure implied by a parameter matrix.

Usage

Gamma2graph(Gamma, tol = get_large_tol(), check = TRUE)

Sigma2graph(Sigma, tol = get_large_tol(), k = NULL, full = FALSE, check = TRUE)

Theta2graph(Theta, tol = get_large_tol(), k = NULL, full = FALSE, check = TRUE)

partialMatrixToGraph(M)

Arguments

Value

An igraph::graph object.

See Also

Other parameter matrix transformations: Gamma2Sigma(), chi2Gamma(), par2Matrix()

Compute exponent measure

Description

Computes the exponent measure of HR distribution.

Usage

V_HR(x, Gamma = NULL, Theta = NULL)

Arguments

Details

Only Gamma is needed for the computation. Theta is only used to compute Gamma if necessary.

Value

Numeric. The exponent measure of the HR distribution.

Computes the Z-matrix

Description

Copied from the R package "golazo" with kind permission by Piotr Zwiernik piotr.zwiernik@utoronto.ca. This function outputs the Z matrix, that is, the unique ultrametric matrix dominating S. This matrix is used to construct a starting point in the GOLAZO algorithm when L=0 but U has strictly positive (off-diagonal entries).

Usage

Zmatrix(S)

Arguments

Censor dataset

Description

Censors each row of matrix x with vector p.

Usage

censor(x, p)

Arguments

Value

Numeric n \times d matrix.

HR parameter matrix checks

Description

Checks whether the matrix given is a valid Huesler-Reiss parameter matrix in the role of \Gamma, \Theta, or \Sigma, respectively.

Usage

checkGamma(Gamma, alert = NULL, tol = get_small_tol(), returnBoolean = FALSE)

checkSigmaTheta(
  M,
  k,
  full,
  matrixName = "Sigma",
  tol = get_small_tol(),
  alert = NULL,
  returnBoolean = FALSE
)

checkTheta(
  Theta,
  k = NULL,
  full = FALSE,
  tol = get_small_tol(),
  alert = NULL,
  returnBoolean = FALSE
)

checkSigma(
  Sigma,
  k = NULL,
  full = FALSE,
  tol = get_small_tol(),
  alert = NULL,
  returnBoolean = FALSE
)

checkMatrix(
  M,
  name = c("Gamma", "Sigma", "Theta")[1],
  k = NULL,
  full = FALSE,
  tol = get_small_tol(),
  alert = NULL,
  returnBoolean = FALSE
)

is_valid_Gamma(M, tol = get_small_tol())

is_valid_Theta(Theta, k = NULL, full = FALSE, tol = get_small_tol())

is_valid_Sigma(Sigma, k = NULL, full = FALSE, tol = get_small_tol())

Arguments

Details

The function ⁠is_valid_*⁠ are a wrapper around ⁠check*⁠, with arguments alert=FALSE and returnBoolean=TRUE.

Value

For ⁠check*⁠, the input matrix, passed through ensure_matrix_symmetry_and_truncate_zeros.

For ⁠is_valid_*⁠, a boolean indicating whether the input is a valid parameter matrix.

See Also

Other input validation functions: check_graph(), check_partial_matrix_and_graph(), ensure_matrix_symmetry()

Check input graph

Description

Checks that the input graph is a valid graph for an extremal graphical model. If necessary, converts the graph to an undirected graph. Removes vertex labels if present.

Usage

check_graph(
  graph,
  graph_type = c("general", "decomposable", "block", "tree"),
  check_connected = TRUE,
  nVertices = NULL
)

Arguments

Value

The given graph, if necessary converted to undirected. If the graph is not valid an error is thrown.

See Also

Other input validation functions: checkGamma(), check_partial_matrix_and_graph(), ensure_matrix_symmetry()

Check input graph and partial matrix

Description

Checks and converts the partial matrix and graph given for a HR graphical model.

Usage

check_partial_matrix_and_graph(M, graph = NULL, graph_type = "general")

Arguments

Value

A list consisting of

Throws an error if the input is not valid.

See Also

Other input validation functions: checkGamma(), check_graph(), ensure_matrix_symmetry()

Identify pairs of vertices that are split by a separator

Description

Identify pairs of vertices that are split by a separator

Usage

check_split_by_sep(graph, sep, edgeMat)

Arguments

Value

A logical vector, indicating for each edge whether it is split by sep

Transformation between \chi and \Gamma

Description

Transforms between the extremal correlation \chi and the variogram \Gamma. Only valid for Huesler-Reiss distributions. Done element-wise, no checks of the entire matrix structure are performed.

Usage

chi2Gamma(chi)

Gamma2chi(Gamma)

Arguments

Details

The formula for transformation from \chi to \Gamma is element-wise

\Gamma = (2 \Phi^{-1}(1 - 0.5 \chi))^2,

where \Phi^{-1} is the inverse of the standard normal distribution function.

The formula for transformation from \Gamma to \chi is element-wise

\chi = 2 - 2 \Phi(\sqrt{\Gamma} / 2),

where \Phi is the standard normal distribution function.

Value

Numeric vector or matrix containing the implied \Gamma.

Numeric vector or matrix containing the implied \chi.

See Also

Other parameter matrix transformations: Gamma2Sigma(), Gamma2graph(), par2Matrix()

Completion of Gamma matrices

Description

Given a graph and a (partial) variogram matrix Gamma, returns a full variogram matrix that agrees with Gamma in entries corresponding to edges of graph and whose corresponding precision matrix, obtained by Gamma2Theta(), has zeros in entries corresponding to non-edges of graph. For results on the existence and uniqueness of this completion, see Hentschel et al. (2022).

Usage

complete_Gamma(Gamma, graph = NULL, ...)

Arguments

Details

If graph is decomposable, Gamma only needs to be specified on the edges of the graph, other entries are ignored. If graph is not decomposable, the graphical completion algorithm requires a fully specified (but non-graphical) variogram matrix Gamma to begin with.

If not initial completion is provided, the function edmcr::npf() can be used to compute one. The package edmcr might need to be installed manually from GitHub.

Value

Completed d \times d variogram matrix.

References

Hentschel M, Engelke S, Segers J (2022). “Statistical Inference for Hüsler-Reiss Graphical Models Through Matrix Completions.” doi:10.48550/ARXIV.2210.14292, https://arxiv.org/abs/2210.14292.

See Also

Gamma2Theta()

Other matrix completion related topics: complete_Gamma_decomposable(), complete_Gamma_general(), complete_Gamma_general_demo(), complete_Gamma_general_split()

Examples

## Block graph:
Gamma <- rbind(
  c(0, .5, NA, NA),
  c(.5, 0, 1, 1.5),
  c(NA, 1, 0, .8),
  c(NA, 1.5, .8, 0)
)

complete_Gamma(Gamma)

## Alternative representation of the same completion problem:
my_graph <- igraph::graph_from_adjacency_matrix(rbind(
  c(0, 1, 0, 0),
  c(1, 0, 1, 1),
  c(0, 1, 0, 1),
  c(0, 1, 1, 0)
), mode = "undirected")
Gamma_vec <- c(.5, 1, 1.5, .8)
complete_Gamma(Gamma_vec, my_graph)

## Decomposable graph:
G <- rbind(
c(0, 5, 7, 6, NA),
c(5, 0, 14, 15, NA),
c(7, 14, 0, 5, 5),
c(6, 15, 5, 0, 6),
c(NA, NA, 5, 6, 0)
)

complete_Gamma(G)

## Non-decomposable graph:
G <- rbind(
c(0, 5, 7, 6, 6),
c(5, 0, 14, 15, 13),
c(7, 14, 0, 5, 5),
c(6, 15, 5, 0, 6),
c(6, 13, 5, 6, 0)
)
g <- igraph::make_ring(5)

complete_Gamma(G, g)

Completion of decomposable Gamma matrices

Description

Given a decomposable graph and incomplete variogram matrix Gamma, returns the full Gamma matrix implied by the conditional independencies.

Usage

complete_Gamma_decomposable(Gamma, graph = NULL)

Arguments

Value

A complete variogram matrix that agrees with Gamma on the entries corresponding to edges in graph and the diagonals. The corresponding \Theta matrix produced by Gamma2Theta() has zeros in the remaining entries.

See Also

Other matrix completion related topics: complete_Gamma(), complete_Gamma_general(), complete_Gamma_general_demo(), complete_Gamma_general_split()

Non-decomposable completion of variogram matrices

Description

Given a non-decomposable graph, and (non-graphical) variogram matrix Gamma, modifies Gamma in non-edge entries, such that the resulting matrix is a variogram matrix with graphical structure described by graph.

Usage

complete_Gamma_general(
  Gamma,
  graph,
  N = 10000,
  tol = get_large_tol(),
  check_tol = 100
)

Arguments

Value

A completed d \times d variogram matrix.

See Also

Other matrix completion related topics: complete_Gamma(), complete_Gamma_decomposable(), complete_Gamma_general_demo(), complete_Gamma_general_split()

DEMO-VERSION: Completion of non-decomposable Gamma matrices

Description

Given a graph and variogram matrix Gamma, returns the full Gamma matrix implied by the conditional independencies. DEMO VERSION: Returns a lot of details and allows specifying the graph list that is used. Is way slower than other functions.

Usage

complete_Gamma_general_demo(Gamma, graph, N = 1000, tol = 0, gList = NULL)

Arguments

Value

A nested list, containing the following details. The "error term" is the maximal absolute value of Theta in a non-edge entry.

See Also

Other matrix completion related topics: complete_Gamma(), complete_Gamma_decomposable(), complete_Gamma_general(), complete_Gamma_general_split()

Non-decomposable completion of variogram matrices

Description

Given a non-decomposable graph, and (non-graphical) variogram matrix Gamma, modifies Gamma in non-edge entries, such that the resulting matrix is a variogram matrix with graphical structure described by graph. Does so by splitting graph at complete separators into smaller subgraphs, and calling complete_Gamma_general for each subgraph/submatrix, using multiple cores if available.

Usage

complete_Gamma_general_split(
  Gamma,
  graph,
  N = 10000,
  sub_tol = get_large_tol() * 0.001,
  check_tol = 100,
  mc_cores_overwrite = NULL,
  final_tol = get_large_tol()
)

Arguments

Value

A completed d \times d variogram matrix.

See Also

Other matrix completion related topics: complete_Gamma(), complete_Gamma_decomposable(), complete_Gamma_general(), complete_Gamma_general_demo()

Completion of two-clique decomposable Gamma matrices

Description

Given a decomposable graph consisting of two cliques and incomplete variogram matrix Gamma, returns the full Gamma matrix implied by the conditional independencies. The rows/columns of Gamma must be ordered such that the clique of size nA (excluding separator) comes first, then the separator of size nC, and then the remaining nB vertices.

Usage

complete_Gamma_one_step(Gamma, nA, nC, nB)

Compute plot limits

Description

Helper function to compute the axis limits of a plot with given x, y data and optionally a fixed x-y-ratio and correcting the latitude/longitude scale at different latitudes

Usage

computeLimits(xData, yData, xyRatio = 1, convertLatLong = TRUE, stretch = 1)

Upper Danube basin dataset

Description

A dataset containing river discharge data for tributaries of the Danube.

Usage

danube

Format

A named list with four entries

data_clustered

A numeric matrix, containing pre-processed discharge data for each gauging station

data_raw

A numeric matrix, containing daily (raw) discharge data for each gauging station

info

A data frame, containing information about each gauging station

flow_edges

A two-column numeric matrix. Each row contains the indices (in info) of a pair of gauging stations that are directly connected by a river.

Details

To obtain the matrix data_clustered, the daily discharge data from the summer months of 1960 to 2010, given in data_raw, was declustered, yielding between seven and ten observations per year. Each row corresponds to one observation from this declustered time series, the non-unique rownames indicate which year an observation is from. Each column corresponds to one of the gauging stations, with column indices in data_raw/data_clustered corresponding to row indices in info. See (Asadi et al. 2015) for details on the preprocessing and declustering.

info is a data frame containing the following information for each of the gauging stations or its corresponding catchment area:

RivNames

Name of the river at the gauging station

Lat, Long

Coordinates of the gauging station

Lat_Center, Long_Center

Coordinates of the center of the catchment corresponding to the gauging station

Alt

Mean altitude of the catchment

Area

Area of the catchment corresponding to the gauging station

Slope

Mean slope of the catchment

PlotCoordX, PlotCoordY

X-Y-coordinates which can be used to arrange the gauging stations when plotting a flow graph.

Source

Bavarian Environmental Agency https://www.gkd.bayern.de.

References

Asadi P, Davison AC, Engelke S (2015). “Extremes on river networks.” Ann. Appl. Stat., 9(4), 2023 – 2050. doi:10.1214/15-AOAS863.

See Also

Other danubeData: getDanubeFlowGraph(), plotDanube()

Other datasets: flights

Examples

dim(danube$data_clustered)
colnames(danube$info)

Data standardization to multivariate Pareto scale

Description

Transforms the data matrix empirically to the multivariate Pareto scale.

Usage

data2mpareto(data, p, na.rm = FALSE)

Arguments

Details

The columns of the data matrix are first transformed empirically to standard Pareto distributions. Then, only the observations where at least one component exceeds the p-quantile of the standard Pareto distribution are kept. Those observations are finally divided by the p-quantile of the standard Pareto distribution to standardize them to the multivariate Pareto scale.

If na.rm is FALSE, missing entries are left as such during the transformation of univariate marginals. In the thresholding step, missing values are considered as -Inf.

Value

Numeric m \times d matrix, where m is the number of rows in the original data matrix that are above the threshold.

See Also

Other parameter estimation methods: emp_chi(), emp_chi_multdim(), emp_vario(), emtp2(), fmpareto_HR_MLE(), fmpareto_graph_HR(), loglik_HR()

Other structure estimation methods: eglatent(), eglearn(), emst(), fit_graph_to_Theta()

Examples

n <- 20
d <- 4
p <- .8
G <- cbind(
  c(0, 1.5, 1.5, 2),
  c(1.5, 0, 2, 1.5),
  c(1.5, 2, 0, 1.5),
  c(2, 1.5, 1.5, 0)
)

set.seed(123)
my_data <- rmstable(n, "HR", d = d, par = G)
data2mpareto(my_data, p)

Learning extremal graph structure with latent variables

Description

Following the methodology from Engelke and Taeb (2024), fits an extremal graph structure with latent variables.

Usage

eglatent(
  Gamma,
  lam1_list = c(0.1, 0.15, 0.19, 0.205),
  lam2_list = c(2),
  refit = TRUE,
  verbose = FALSE
)

Arguments

Value

The function fits one model for each combination of values in lam1_list and lam2_list. All returned objects have one entry per model. List consisting of:

References

Engelke S, Taeb A (2024). “Extremal graphical modeling with latent variables.” 2403.09604.

See Also

Other structure estimation methods: data2mpareto(), eglearn(), emst(), fit_graph_to_Theta()

Learning extremal graph structure

Description

Following the methodology from Engelke et al. (2022), fits an extremal graph structure using the neighborhood selection approach (see Meinshausen and Bühlmann (2006)) or graphical lasso (see Friedman et al. (2008)).

Usage

eglearn(
  data,
  p = NULL,
  rholist = c(0.1, 0.15, 0.19, 0.205),
  reg_method = c("ns", "glasso"),
  complete_Gamma = FALSE
)

Arguments

Value

List made of:

References

Engelke S, Lalancette M, Volgushev S (2022). “Learning extremal graphical structures in high dimensions.” doi:10.48550/ARXIV.2111.00840, Available from https://arxiv.org/abs/2111.00840., 2111.00840, https://arxiv.org/abs/2111.00840.

Friedman J, Hastie T, Tibshirani R (2008). “Sparse inverse covariance estimation with the graphical lasso.” Biostatistics, 9(3), 432–441.

Meinshausen N, Bühlmann P (2006). “High-dimensional graphs and variable selection with the Lasso.” Ann. Statist., 34(3), 1436 – 1462. doi:10.1214/009053606000000281.

See Also

Other structure estimation methods: data2mpareto(), eglatent(), emst(), fit_graph_to_Theta()

Examples

set.seed(2)
m <- generate_random_model(d=6)
y <- rmpareto(n=500, par=m$Gamma)
ret <- eglearn(y)

Empirical estimation of extremal correlation matrix \chi

Description

Estimates empirically the matrix of bivariate extremal correlation coefficients \chi.

Usage

emp_chi(data, p = NULL)

emp_chi_pairwise(data, p = NULL, verbose = FALSE)

Arguments

Details

emp_chi_pairwise calls emp_chi for each pair of observations. This is more robust if the data contains many NAs, but can take rather long.

Value

Numeric matrix d \times d. The matrix contains the bivariate extremal coefficients \chi_{ij}, for i, j = 1, ..., d.

See Also

Other parameter estimation methods: data2mpareto(), emp_chi_multdim(), emp_vario(), emtp2(), fmpareto_HR_MLE(), fmpareto_graph_HR(), loglik_HR()

Examples

n <- 100
d <- 4
p <- .8
Gamma <- cbind(
  c(0, 1.5, 1.5, 2),
  c(1.5, 0, 2, 1.5),
  c(1.5, 2, 0, 1.5),
  c(2, 1.5, 1.5, 0)
)

set.seed(123)
my_data <- rmstable(n, "HR", d = d, par = Gamma)
emp_chi(my_data, p)

Empirical estimation of extremal correlation \chi

Description

Estimates the d-dimensional extremal correlation coefficient \chi empirically.

Usage

emp_chi_multdim(data, p = NULL)

Arguments

Value

Numeric scalar. The empirical d-dimensional extremal correlation coefficient \chi for the data.

See Also

Other parameter estimation methods: data2mpareto(), emp_chi(), emp_vario(), emtp2(), fmpareto_HR_MLE(), fmpareto_graph_HR(), loglik_HR()

Examples

n <- 100
d <- 2
p <- .8
G <- cbind(
  c(0, 1.5),
  c(1.5, 0)
)

set.seed(123)
my_data <- rmstable(n, "HR", d = d, par = G)
emp_chi_multdim(my_data, p)

Estimation of the variogram matrix \Gamma of a Huesler-Reiss distribution

Description

Estimates the variogram of the Huesler-Reiss distribution empirically.

Usage

emp_vario(data, k = NULL, p = NULL)

emp_vario_pairwise(data, k = NULL, p = NULL, verbose = FALSE)

Arguments

Details

emp_vario_pairwise calls emp_vario for each pair of observations. This is more robust if the data contains many NAs, but can take rather long.

Value

Numeric d \times d matrix. The estimated variogram of the Huesler-Reiss distribution.

See Also

Other parameter estimation methods: data2mpareto(), emp_chi(), emp_chi_multdim(), emtp2(), fmpareto_HR_MLE(), fmpareto_graph_HR(), loglik_HR()

Examples

G <- generate_random_Gamma(d=5)
y <- rmpareto(n=100, par=G)
Ghat <- emp_vario(y)

Fitting extremal minimum spanning tree

Description

Fits an extremal minimum spanning tree, where the edge weights are:

• negative maximized log-likelihoods of the bivariate Huesler-Reiss distributions, if method = "ML". See Engelke and Hitz (2020) for details.

• empirical extremal variogram, if method = "vario". See Engelke and Volgushev (2022) for details.

• empirical extremal correlation, if method = "chi". See Engelke and Volgushev (2022) for details.

Usage

emst(data, p = NULL, method = c("vario", "ML", "chi"), cens = FALSE)

Arguments

Value

List consisting of:

References

Engelke S, Hitz AS (2020). “Graphical models for extremes (with discussion).” J. R. Stat. Soc. Ser. B Stat. Methodol., 82, 871–932.

Engelke S, Volgushev S (2022). “Structure learning for extremal tree models.” J. R. Stat. Soc. Ser. B Stat. Methodol.. doi:10.1111/rssb.12556, Forthcoming, https://rss.onlinelibrary.wiley.com/doi/pdf/10.1111/rssb.12556.

See Also

Other structure estimation methods: data2mpareto(), eglatent(), eglearn(), fit_graph_to_Theta()

Examples

## Fitting a 4-dimensional HR minimum spanning tree
my_graph <- igraph::graph_from_adjacency_matrix(
  rbind(
    c(0, 1, 0, 0),
    c(1, 0, 1, 1),
    c(0, 1, 0, 0),
    c(0, 1, 0, 0)
  ),
  mode = "undirected"
)
n <- 100
Gamma_vec <- c(.5, 1.4, .8)
complete_Gamma(Gamma = Gamma_vec, graph = my_graph) ## full Gamma matrix

set.seed(123)
my_data <- rmpareto_tree(n, "HR", tree = my_graph, par = Gamma_vec)
my_fit <- emst(my_data, p = NULL, method = "ML", cens = FALSE)

Performs Gaussian likelihood optimization under Laplacian matrix constraints.

Description

This function implements a block descent algorithm to find the maximum of the Gaussian log-likelihood under the constraint that the concentration matrix is a Laplacian matrix. See Röttger et al. (2021) for details.

Usage

emtp2(Gamma, tol = 1e-06, verbose = TRUE, initial_point = TRUE)

Arguments

Value

A list consisting of:

References

Röttger F, Engelke S, Zwiernik P (2021). “Total positivity in multivariate extremes.” doi:10.48550/ARXIV.2112.14727, https://arxiv.org/abs/2112.14727.

See Also

Other parameter estimation methods: data2mpareto(), emp_chi(), emp_chi_multdim(), emp_vario(), fmpareto_HR_MLE(), fmpareto_graph_HR(), loglik_HR()

Ensure numerical matrix symmetry/zero values

Description

Ensures the symmetry of a square matrix by averaging it with its transpose. Optionally verifies that the matrix was close to symmetric before.

Makes sure zeros are "numerically zero", by truncating all small values.

Usage

ensure_matrix_symmetry(M, checkTol = Inf, alert = NULL)

truncate_zeros(M, tol = get_small_tol())

ensure_matrix_symmetry_and_truncate_zeros(
  M,
  tol = get_small_tol(),
  checkTol = Inf
)

Arguments

Value

The adjusted value of M.

See Also

Other input validation functions: checkGamma(), check_graph(), check_partial_matrix_and_graph()

Fast computation of diag(y %% M %% t(y))

Description

Fast computation of diag(y %% M %% t(y))

Usage

fast_diag(y, M)

Arguments

Value

Numeric vector

Helper function to combine par with fixed params (in init)

Description

Helper function to combine par with fixed params (in init)

Usage

fillFixedParams(par, init, fixParams)

Arguments

Find a separator set for two vertices

Description

Finds a reasonably small set of vertices that separate v0 and v1 in graph.

Usage

findVsep(graph, v0, v1)

Fit value(s) in interval

Description

Fit value(s) in interval, all arguments are recycled where necessary.

Usage

fitInInterval(x, xMin = -Inf, xMax = Inf)

Arguments

Value

Numeric vector

Experimental: Fit graph using empirical Theta matrix

Description

Fits a graph to an empirical Gamma matrix by computing the corresponding Theta matrix using Gamma2Theta() and greedily chooses m edges that correspond to high absolute values in Theta.

Usage

fit_graph_to_Theta(data, m = NULL, Gamma_emp = NULL)

Arguments

Value

A list containing an [igraph::graph] object and a fitted Gamma matrix

See Also

Other structure estimation methods: data2mpareto(), eglatent(), eglearn(), emst()

Convert flight counts to connection list

Description

Convert a numeric matrix containing flight counts between airports to a data frame containing a list of connections.

Usage

flightCountMatrixToConnectionList(nFlightsPerConnection, directed = TRUE)

Arguments

Value

A data frame with columns departureAirport, arrivalAirport, nFlights. Each row represents one connection with >=1 flights in the input matrix.

See Also

Other flight data related topics: flights, getFlightDelayData(), getFlightGraph(), plotFlights()

Examples

flightCountMatrixToConnectionList(flights$flightCounts[1:100, 1:100, 1])

Flights delay data

Description

A dataset containing daily total delays of major U.S. airlines. The raw data was obtained from the U.S. Bureau of Transportation Statistics, and pre-processed as described in Hentschel et al. (2022). Note: The CRAN version of this package contains only data from 2010-2013. The full dataset is available in the GitHub version of this package.

Usage

flights

Format

A named list with three entries:

airports

A data.frame, containing information about US airports

delays

A numeric matrix, containing daily aggregated delays at the airports in the dataset

flightCounts

A numeric array, containing yearly flight numbers between airports in the dataset

Details

flightCounts is a three-dimensional array, containing the number of flights in the dataset between each pair of airports, aggregated on a yearly basis. Each entry is the total number of flights between the departure airport (row) and destination airport (column) in a given year (dimension 3). This array does not contain any NAs, even if an airport did not operate at all in a given year, which is simply indicated by zeros.

delays is a three-dimensional array containing daily total positive delays, in minutes, of incoming and outgoing flights respectively. Each column corresponds to an airport in the dataset and each row corresponds to a day. The third dimension has length two, 'arrivals' containing delays of incoming flights and 'departures' containing delays of outgoing flights. Zeros indicate that there were flights arriving/departing at that airport on a given day, but none of them had delays. NAs indicate that there were no flights arriving/departing at that airport on that day at all.

airports is a data frame containing the following information about a number of US airports. Some entries are missing, which is indicated by NAs.

IATA

3-letter IATA code

Name

name of the airport

City

main city served by the airport

Country

country or territory where the airport is located (mostly "United States")

ICAO

4-letter ICAO code

Latitude

latitude of the airport, in decimal degrees

Longitude

longitude of the airport, in decimal degrees

Altitude

altitude of the airport, in feet

Timezone

timezone of the airport, in hours offset from UTC

DST

Daylight savings time used at the airport. 'A'=US/Canada, 'N'=None.

Timezone2

name of the timezone of the airport

Source

Raw delays data:

• https://www.bts.gov/browse-statistical-products-and-data/bts-publications/airline-service-quality-performance-234-time

Fields/Forms used in the raw data:

• https://esubmit.rita.dot.gov/ViewReports.aspx

• https://esubmit.rita.dot.gov/On-Time-Form1.aspx

• https://esubmit.rita.dot.gov/On-Time-Form3A.aspx

Airports (includes license information):

• https://openflights.org/data

References

Hentschel M, Engelke S, Segers J (2022). “Statistical Inference for Hüsler-Reiss Graphical Models Through Matrix Completions.” doi:10.48550/ARXIV.2210.14292, https://arxiv.org/abs/2210.14292.

See Also

Other flight data related topics: flightCountMatrixToConnectionList(), getFlightDelayData(), getFlightGraph(), plotFlights()

Other datasets: danube

Examples

# Get total number of flights in the dataset:
totalFlightCounts <- apply(flights$flightCounts, c(1,2), sum)

# Get number of flights for specific years in the dataset:
flightCounts_10_11 <- apply(flights$flightCounts[,,c('2010', '2011')], c(1,2), sum)

# Get list of connections from 2008:
connections_10 <- flightCountMatrixToConnectionList(flights$flightCounts[,,'2010'])

Parameter fitting for multivariate Huesler-Reiss Pareto distribution

Description

Fits the parameters of a multivariate Huesler-Reiss Pareto distribution using (censored) maximum likelihood estimation.

Usage

fmpareto_HR_MLE(
  data,
  p = NULL,
  cens = FALSE,
  init = NULL,
  fixParams = integer(0),
  useTheta = TRUE,
  maxit = 100,
  graph = NULL,
  optMethod = "BFGS",
  nAttemptsFixInit = 3
)

Arguments

Details

Only the parameters corresponding to edges in graph are optimized, the remaining entries are implied by the graphical structure. If graph is NULL, the complete graph is used. The optimization is done either in terms of the variogram (Gamma) or precision matrix (Theta), depending on the value of useTheta. If graph is non-decomposable, useTheta=TRUE is significantly faster, otherwise they are similar in performance.

Value

List consisting of:

See Also

Other parameter estimation methods: data2mpareto(), emp_chi(), emp_chi_multdim(), emp_vario(), emtp2(), fmpareto_graph_HR(), loglik_HR()

Parameter fitting for Huesler-Reiss graphical models

Description

Fits the parameter matrix (variogram) of a multivariate Huesler-Reiss Pareto distribution with a given graphical structure, using maximum-likelihood estimation or the empirical variogram.

Usage

fmpareto_graph_HR(
  data,
  graph,
  p = NULL,
  method = c("vario", "ML"),
  handleCliques = c("average", "full", "sequential"),
  ...
)

Arguments

Details

If handleCliques='average', the marginal parameter matrix is estimated for each maximal clique of the graph and then combined into a partial parameter matrix by taking the average of entries from overlapping cliques. Lastly, the full parameter matrix is computed using complete_Gamma().

If handleCliques='full', first the full parameter matrix is estimated using the specified method and then the non-edge entries are adjusted such that the final parameter matrix has the graphical structure indicated by graph.

If handleCliques='sequential', graph must be decomposable, and method='ML' must be specified. The parameter matrix is first estimated on the (recursive) separators and then on the rest of the cliques, keeping previously estimated entries fixed.

If method='ML', the computational cost is mostly influenced by the total size of the graph (if handleCliques='full') or the size of the cliques, and can already take a significant amount of time for modest dimensions (e.g. d=3).

Value

The estimated parameter matrix.

See Also

Other parameter estimation methods: data2mpareto(), emp_chi(), emp_chi_multdim(), emp_vario(), emtp2(), fmpareto_HR_MLE(), loglik_HR()

HR Parameter fitting - Helper functions

Description

Helper functions called by fmpareto_HR_MLE.

Usage

fmpareto_graph_HR_clique_average(data, graph, method = c("ML", "vario"), ...)

fmpareto_graph_HR_clique_sequential(data, graph, ...)

combine_clique_estimates_by_averaging(cliques, subGammas)

Generate a random Gamma matrix

Description

Generates a valid Gamma matrix with a given dimension

Usage

generate_random_Gamma(d, ...)

Arguments

See Also

Other example generation functions: generate_random_chordal_graph(), generate_random_graphical_Gamma(), generate_random_integer_Gamma(), generate_random_model(), generate_random_spd_matrix()

Generate random graphs

Description

Generate random graphs with different structures. These do not follow well-defined distributions and are mostly meant for quickly generating test models.

Usage

generate_random_chordal_graph(
  d,
  cMin = 2,
  cMax = 6,
  sMin = 1,
  sMax = 4,
  block_graph = FALSE,
  ...
)

generate_random_connected_graph(
  d,
  m = NULL,
  p = 2/(d + 1),
  maxTries = 1000,
  ...
)

generate_random_tree(d)

generate_random_cactus(d, cMin = 2, cMax = 6)

Arguments

Details

generate_random_chordal_graph generates a random chordal graph by starting with a (small) complete graph and then adding new cliques until the specified size is reached. The sizes of cliques and separators can be specified.

generate_random_connected_graph first tries to generate an Erdoes-Renyi graph, if that fails, falls back to producing a tree and adding random edges to that tree.

generate_random_cactus generates a random cactus graph (mostly useful for benchmarking).

Value

An [igraph::graph] object

See Also

Other example generation functions: generate_random_Gamma(), generate_random_graphical_Gamma(), generate_random_integer_Gamma(), generate_random_model(), generate_random_spd_matrix()

Generate a random Gamma matrix for a given graph

Description

Generates a valid Gamma matrix with conditional independence structure specified by a graph

Usage

generate_random_graphical_Gamma(graph, ...)

Arguments

See Also

Other example generation functions: generate_random_Gamma(), generate_random_chordal_graph(), generate_random_integer_Gamma(), generate_random_model(), generate_random_spd_matrix()

Generate a random Gamma matrix containing only integers

Description

Generates a random variogram Matrix by producing a (d-1) \times (d-1) matrix B with random integer entries between -b and b, computing S = B %*% t(B), and passing this S to Sigma2Gamma(). This process is repeated with an increasing b until a valid Gamma matrix is produced.

Usage

generate_random_integer_Gamma(d, b = 2, b_step = 1)

Arguments

Value

A numeric d \times d variogram matrix with integer entries

See Also

Other example generation functions: generate_random_Gamma(), generate_random_chordal_graph(), generate_random_graphical_Gamma(), generate_random_model(), generate_random_spd_matrix()

Examples

generate_random_integer_Gamma(5, 2, 0.1)

Generate random Huesler-Reiss Models

Description

Generates a random connected graph and Gamma matrix with conditional independence structure corresponding to that graph.

Usage

generate_random_model(d, graph_type = "general", ...)

Arguments

See Also

Other example generation functions: generate_random_Gamma(), generate_random_chordal_graph(), generate_random_graphical_Gamma(), generate_random_integer_Gamma(), generate_random_spd_matrix()

Examples

set.seed(1)
d <- 12

generate_random_model(d, 'tree')
generate_random_model(d, 'block')
generate_random_model(d, 'decomposable')
generate_random_model(d, 'general')
generate_random_model(d, 'complete')

Generate a random symmetric positive definite matrix

Description

Generates a random d \times d symmetric positive definite matrix. This is done by generating a random d \times d matrix B, then computing B %*% t(B), and then normalizing the matrix to approximately single digit entries.

Usage

generate_random_spd_matrix(d, bMin = -10, bMax = 10, ...)

Arguments

See Also

Other example generation functions: generate_random_Gamma(), generate_random_chordal_graph(), generate_random_graphical_Gamma(), generate_random_integer_Gamma(), generate_random_model()

Get Danube flow graph

Description

Returns an igraph::graph object representing the flow graph of the danube dataset.

Usage

getDanubeFlowGraph(stationIndices = NULL, directed = FALSE)

Arguments

Value

An igraph::graph object.

See Also

Other danubeData: danube, plotDanube()

Get filtered flight delays

Description

Get filtered flight delay data, containing only a selection of dates and airports. Currently, all possible selections correspond to the case study in Hentschel et al. (2022).

Usage

getFlightDelayData(
  what = c("delays", "IATAs", "dates"),
  airportFilter = c("all", "tcCluster", "tcAll"),
  dateFilter = c("all", "tcTrain", "tcTest", "tcAll"),
  delayFilter = c("totals", "arrivals", "departures")[1]
)

Arguments

Details

The provided lists of airports and dates correspond to the ones used in the case study of Hentschel et al. (2022). The argument airportFilter="tcCluster" corresponds to the airports in the analyzed "Texas Cluster", airportFilter="tcAll" corresponds to all airports used in the previous clustering step, airportFilter="all" corresponds to all airports in the dataset.

Similarly, dateFilter="tcTrain" selects the dates from the training set, dateFilter="tcTest" the ones from the test/validation set. To get the union of these sets, specify dateFilter="tcAll". To get all dates in the dataset (possibly more than for "tcAll"), specify dateFilter="all".

Value

If what="IATAs" or what="dates", a character vector. If required, it can be converted to Date objects using as.Date().

If what="delays", a three-dimensional array or two-dimensional matrix, with dimensions corresponding to dates, airports, and delay types.

References

Hentschel M, Engelke S, Segers J (2022). “Statistical Inference for Hüsler-Reiss Graphical Models Through Matrix Completions.” doi:10.48550/ARXIV.2210.14292, https://arxiv.org/abs/2210.14292.

See Also

Other flight data related topics: flightCountMatrixToConnectionList(), flights, getFlightGraph(), plotFlights()

Get flight graph

Description

Convert the info from flights$flightCounts to an igraph::graph object.

Usage

getFlightGraph(IATAs = NULL, years = NULL, minNFlights = 1, directed = FALSE)

Arguments

Value

An igraph::graph object containing a vertex for each airport and an edge whenever there are at least minNFlights between two airports.

See Also

Other flight data related topics: flightCountMatrixToConnectionList(), flights, getFlightDelayData(), plotFlights()

Examples

g <- getFlightGraph()

Get package data

Description

Get private data sets from inst/extdata.

Usage

getPackageData(filename, isRDS = TRUE)

Arguments

Value

If isRDS=TRUE an R object. If isRDS=FALSE an environment, containing the R objects from the file.

Get the submatrix corresponding to a subgraph

Description

Both the graph and subgraph need to have persistent IDs If graph==NULL it is assumed to have pIDs 1, 2, ...

Usage

getSubMatrixForSubgraph(fullMatrix, subgraph, graph = NULL)

Get alert function

Description

Get a function that can be used to alert the user of invalid inputs. Returns the value implied by the overwrite argument, or the option "graphicalExtremes.default.alert", falling back to warning() if neither is specified.

Usage

get_alert_function(overwrite = NULL)

Arguments

Value

A function that takes an arbitrary number of strings as arguments.

See Also

graphicalExtremes-package

Other default parameters: get_mc_cores(), get_small_tol()

Get Cliques and Separators of a graph

Description

Finds all cliques, separators, and (recursively) separators of separators in a graph.

Usage

get_cliques_and_separators(
  graph,
  sortIntoLayers = FALSE,
  includeSingletons = FALSE
)

Arguments

Value

A list of vertex sets that represent the cliques and (recursive) separators of graph, ordered such that separators come before cliques they separate.

Get Cliques and Separators of a graph

Description

Finds all cliques, separators, and (recursively) separators of separators in a graph.

Usage

get_cliques_and_separators_OLD(graph)

Arguments

Value

A list of vertex sets that represent the cliques and (recursive) separators of graph

Number of cores to be used in parallel computations

Description

Helper function that returns the number of cores to be used in parallel computations. Will always be 1 on Windows. On other systems, this value can be set using setOption('graphicalExtremes.mc.cores', ...).

Usage

get_mc_cores(overwrite = NULL)

Arguments

Value

An integer to be used as number of cores

See Also

graphicalExtremes-package

Other default parameters: get_alert_function(), get_small_tol()

Tolerances to be used in computations

Description

Helper function that returns the tolerance to be used in internal computations.

Usage

get_small_tol(overwrite = NULL)

get_large_tol(overwrite = NULL)

Arguments

Details

There are two different tolerances used in the package, for details see graphicalExtremes-package. The default values for these tolerances can be set using the options "graphicalExtremes.tol.small" and "graphicalExtremes.tol.large".

Value

A non-negative numerical scalar

See Also

graphicalExtremes-package

Other default parameters: get_alert_function(), get_mc_cores()

Graph equality

Description

Produce true if two graphs have same vertices and edges (ordered)

Usage

graphs_equal(g1, g2)

Arguments

Value

logical indicating if the graphs are equal

Full censored log-likelihood of HR model

Description

Computes the full (censored) log-likelihood of HR model.

Usage

logLH_HR(data, Gamma, cens = FALSE)

Arguments

Value

Numeric. The full censored log-likelihood of HR model.

Compute censored exponent measure

Description

Computes the (censored) exponent measure density of HR distribution.

Usage

logdVK_HR(x, K, Gamma)

Arguments

Value

Numeric. The censored exponent measure of the HR distribution. If no entries are censored, the result of ⁠logdV_HR(x, par⁠ is returned.

Compute the exponent measure density of HR distribution

Description

Computes the exponent measure density of HR distribution.

Usage

logdV_HR(x, Gamma = NULL, Theta = NULL)

Arguments

Details

Both Gamma and Theta are needed internally, but if one is missing it is computed from the other one.

Value

Numeric. The censored exponent measure of the HR distribution.

Compute Huesler-Reiss log-likelihood, AIC, and BIC

Description

Computes (censored) Huesler-Reiss log-likelihood, AIC, and BIC values.

Usage

loglik_HR(data, p = NULL, graph = NULL, Gamma, cens = FALSE)

Arguments

Value

Numeric vector c("loglik"=..., "aic"=..., "bic"=...) with the evaluated log-likelihood, AIC, and BIC values.

See Also

Other parameter estimation methods: data2mpareto(), emp_chi(), emp_chi_multdim(), emp_vario(), emtp2(), fmpareto_HR_MLE(), fmpareto_graph_HR()

Convert indices to numerical indices

Description

Converts (possibly) logical indices to numerical ones. Also ensures unique indices and sorts them if specified.

Usage

make_numeric_indices(ind, n = NULL, unique = TRUE, sort = TRUE)

Arguments

Value

A numerical index vector

Create a list of separators

Description

Creates a list of separator set, such that every pair of non-adjacent vertices in graph is completely disconnected by the removal of (at least) one of the separator sets from the graph.

Usage

make_sep_list(graph, details = TRUE)

Arguments

Value

A list of numeric vectors

Marginalize multivariate Pareto dataset

Description

Marginalize a multivariate Pareto dataset data with respect to the variables in set_indices.

Usage

mparetomargins(data, set_indices)

Arguments

Value

Numeric n \times m matrix, where m is the length of set_indices. Marginalized multivariate Pareto data.

Order Cliques

Description

Orders the cliques in a connected decomposable graph so that they fulfill the running intersection property.

Usage

order_cliques(cliques)

Create Gamma or Theta from vector

Description

Convert parameter vector par (upper triangular part of Gamma/Theta matrix) to full Gamma/Theta, or vice versa.

Usage

par2Matrix(par, allowMatrix = FALSE, allowNull = FALSE, zeroRowSums = FALSE)

par2Gamma(par, allowMatrix = FALSE, allowNull = FALSE)

par2Theta(par, allowMatrix = FALSE, allowNull = FALSE)

matrix2par(M, allowVector = FALSE, allowNull = FALSE)

Arguments

Value

Numeric matrix d \times d. Full Gamma/Theta matrix corresponding to par.

Upper triangular part of M (or M itself/NULL if allowed)

See Also

Other parameter matrix transformations: Gamma2Sigma(), Gamma2graph(), chi2Gamma()

Factory: parToMatrices

Description

Creates a helper function to convert a parameter vector par to a Gamma and/or Theta matrix.

Usage

parToMatricesFactory(
  graph,
  init = NULL,
  fixParams = integer(0),
  parIsTheta = FALSE,
  checkValidity = TRUE
)

Arguments

Value

A function parToMatrices(par, forceGamma=FALSE, forceTheta=FALSE), which takes a parameter vector and returns either NULL or a list with entries Gamma, Theta. The function returns NULL if checkValidity==TRUE and par implies an invalid matrix. Otherwise, depending on parIsTheta, forceTheta, and forceGamma, one or both of Gamma and Theta are matrices implied by par.

Plot Danube River Flow Data

Description

Plotting function to visualize the river flow data from the danube dataset. Requires ggplot2 to be installed.

Usage

plotDanube(
  stationIndices = NULL,
  graph = NULL,
  directed = NULL,
  plotStations = TRUE,
  plotConnections = TRUE,
  labelStations = FALSE,
  returnGGPlot = FALSE,
  useStationVolume = FALSE,
  useConnectionVolume = FALSE,
  mapCountries = c("Germany"),
  vertexColors = NULL,
  vertexShapes = NULL,
  edgeColors = NULL,
  xyRatio = NULL,
  clipMap = 1.2,
  useLatex = FALSE,
  edgeAlpha = 0.2
)

plotDanubeIGraph(
  stationIndices = NULL,
  graph = NULL,
  directed = NULL,
  labelStations = TRUE,
  vertexColors = NULL,
  vertexShapes = NULL,
  edgeColors = NULL,
  ...
)

Arguments

Details

The values of vertexColors, vertexShapes, and edgeColors are interpreted differently by ggplot2::geom_point/ggplot2::geom_segment and igraph::plot.igraph().

plotDanube uses a combination of ggplot2 functions to plot the graph.

plotDanubeIGraph uses igraph::plot.igraph to plot the graph.

See Also

plotFlights

Other danubeData: danube, getDanubeFlowGraph()

Examples

# Basic plot
graphicalExtremes::plotDanube()

# Plot flow volumes
graphicalExtremes::plotDanube(
    clipMap = 1.2,
    useConnectionVolume = TRUE,
    useStationVolume = TRUE
)

# Plot other graph structures
nStations <- nrow(graphicalExtremes::danube$info)
g <- igraph::erdos.renyi.game(nStations, nStations, 'gnm')
graphicalExtremes::plotDanube(
    clipMap = 1.2,
    graph = g
)

Plot flight data

Description

Plotting function to visualize the flight connections from the flights dataset. This function requires the package ggplot2 to be installed.

Usage

plotFlights(
  airportIndices = NULL,
  airports_sel = NULL,
  connections_sel = NULL,
  graph = NULL,
  plotAirports = TRUE,
  plotConnections = TRUE,
  labelAirports = FALSE,
  returnGGPlot = FALSE,
  useAirportNFlights = FALSE,
  useConnectionNFlights = FALSE,
  minNFlights = 0,
  map = "state",
  vertexColors = NULL,
  vertexShapes = NULL,
  edgeColors = NULL,
  xyRatio = NULL,
  clipMap = FALSE,
  useLatex = FALSE,
  edgeAlpha = 0.2
)

Arguments

Value

If returnGGPlot is TRUE, a ggplot2::ggplot object, otherwise NULL.

See Also

plotDanube

Other flight data related topics: flightCountMatrixToConnectionList(), flights, getFlightDelayData(), getFlightGraph()

Examples

# Plot all airports in the dataset
plotFlights(plotConnections = FALSE, map = 'world')

# Plot a selection of airports
plotFlights(c('JFK', 'SFO', 'LAX'), useConnectionNFlights = TRUE, useAirportNFlights = TRUE)

# Plot airports with a custom connections graph
IATAs <- c('ACV', 'BFL', 'EUG', 'SFO', 'MRY')
graph <- igraph::make_full_graph(length(IATAs))
plotFlights(IATAs, graph=graph, clipMap = 1.5)

Sampling of a multivariate Pareto distribution

Description

Simulates exact samples of a multivariate Pareto distribution.

Usage

rmpareto(
  n,
  model = c("HR", "logistic", "neglogistic", "dirichlet"),
  d = NULL,
  par
)

Arguments

Details

The simulation follows the algorithm in Engelke and Hitz (2020). For details on the parameters of the Huesler-Reiss, logistic and negative logistic distributions see Dombry et al. (2016), and for the Dirichlet distribution see Coles and Tawn (1991).

Value

Numeric n \times d matrix of simulations of the multivariate Pareto distribution.

References

Coles S, Tawn JA (1991). “Modelling extreme multivariate events.” J. R. Stat. Soc. Ser. B Stat. Methodol., 53, 377–392.

Dombry C, Engelke S, Oesting M (2016). “Exact simulation of max-stable processes.” Biometrika, 103, 303–317.

Engelke S, Hitz AS (2020). “Graphical models for extremes (with discussion).” J. R. Stat. Soc. Ser. B Stat. Methodol., 82, 871–932.

See Also

Other sampling functions: rmpareto_tree(), rmstable(), rmstable_tree()

Examples

## A 4-dimensional HR distribution
n <- 10
d <- 4
G <- cbind(
  c(0, 1.5, 1.5, 2),
  c(1.5, 0, 2, 1.5),
  c(1.5, 2, 0, 1.5),
  c(2, 1.5, 1.5, 0)
)

rmpareto(n, "HR", d = d, par = G)

## A 3-dimensional logistic distribution
n <- 10
d <- 3
theta <- .6
rmpareto(n, "logistic", d, par = theta)

## A 5-dimensional negative logistic distribution
n <- 10
d <- 5
theta <- 1.5
rmpareto(n, "neglogistic", d, par = theta)

## A 4-dimensional Dirichlet distribution
n <- 10
d <- 4
alpha <- c(.8, 1, .5, 2)
rmpareto(n, "dirichlet", d, par = alpha)

Sampling of a multivariate Pareto distribution on a tree

Description

Simulates exact samples of a multivariate Pareto distribution that is an extremal graphical model on a tree as defined in Engelke and Hitz (2020).

Usage

rmpareto_tree(n, model = c("HR", "logistic", "dirichlet")[1], tree, par)

Arguments

Details

The simulation follows the algorithm in Engelke and Hitz (2020). For details on the parameters of the Huesler-Reiss, logistic and negative logistic distributions see Dombry et al. (2016), and for the Dirichlet distribution see Coles and Tawn (1991).

Value

Numeric n \times d matrix of simulations of the multivariate Pareto distribution.

References

Coles S, Tawn JA (1991). “Modelling extreme multivariate events.” J. R. Stat. Soc. Ser. B Stat. Methodol., 53, 377–392.

Dombry C, Engelke S, Oesting M (2016). “Exact simulation of max-stable processes.” Biometrika, 103, 303–317.

Engelke S, Hitz AS (2020). “Graphical models for extremes (with discussion).” J. R. Stat. Soc. Ser. B Stat. Methodol., 82, 871–932.

See Also

Other sampling functions: rmpareto(), rmstable(), rmstable_tree()

Examples

## A 4-dimensional HR tree model

my_tree <- igraph::graph_from_adjacency_matrix(rbind(
  c(0, 1, 0, 0),
  c(1, 0, 1, 1),
  c(0, 1, 0, 0),
  c(0, 1, 0, 0)
),
mode = "undirected"
)
n <- 10
Gamma_vec <- c(.5, 1.4, .8)
set.seed(123)
rmpareto_tree(n, "HR", tree = my_tree, par = Gamma_vec)

## A 4-dimensional Dirichlet model with asymmetric edge distributions

alpha <- cbind(c(.2, 1, .5), c(1.5, .6, .8))
rmpareto_tree(n, model = "dirichlet", tree = my_tree, par = alpha)

Sampling of a multivariate max-stable distribution

Description

Simulates exact samples of a multivariate max-stable distribution.

Usage

rmstable(n, model = c("HR", "logistic", "neglogistic", "dirichlet")[1], d, par)

Arguments

Details

The simulation follows the extremal function algorithm in Dombry et al. (2016). For details on the parameters of the Huesler-Reiss, logistic and negative logistic distributions see Dombry et al. (2016), and for the Dirichlet distribution see Coles and Tawn (1991).

Value

Numeric n \times d matrix of simulations of the multivariate max-stable distribution.

References

Coles S, Tawn JA (1991). “Modelling extreme multivariate events.” J. R. Stat. Soc. Ser. B Stat. Methodol., 53, 377–392.

Dombry C, Engelke S, Oesting M (2016). “Exact simulation of max-stable processes.” Biometrika, 103, 303–317.

See Also

Other sampling functions: rmpareto(), rmpareto_tree(), rmstable_tree()

Examples

## A 4-dimensional HR distribution
n <- 10
d <- 4
G <- cbind(
  c(0, 1.5, 1.5, 2),
  c(1.5, 0, 2, 1.5),
  c(1.5, 2, 0, 1.5),
  c(2, 1.5, 1.5, 0)
)

rmstable(n, "HR", d = d, par = G)

## A 3-dimensional logistic distribution
n <- 10
d <- 3
theta <- .6
rmstable(n, "logistic", d, par = theta)

## A 5-dimensional negative logistic distribution
n <- 10
d <- 5
theta <- 1.5
rmstable(n, "neglogistic", d, par = theta)

## A 4-dimensional Dirichlet distribution
n <- 10
d <- 4
alpha <- c(.8, 1, .5, 2)
rmstable(n, "dirichlet", d, par = alpha)

Sampling of a multivariate max-stable distribution on a tree

Description

Simulates exact samples of a multivariate max-stable distribution that is an extremal graphical model on a tree as defined in Engelke and Hitz (2020).

Usage

rmstable_tree(n, model = c("HR", "logistic", "dirichlet")[1], tree, par)

Arguments

Details

The simulation follows a combination of the extremal function algorithm in Dombry et al. (2016) and the theory in Engelke and Hitz (2020) to sample from a single extremal function. For details on the parameters of the Huesler-Reiss, logistic and negative logistic distributions see Dombry et al. (2016), and for the Dirichlet distribution see Coles and Tawn (1991).

Value

Numeric n \times d matrix of simulations of the multivariate max-stable distribution.

References

Coles S, Tawn JA (1991). “Modelling extreme multivariate events.” J. R. Stat. Soc. Ser. B Stat. Methodol., 53, 377–392.

Dombry C, Engelke S, Oesting M (2016). “Exact simulation of max-stable processes.” Biometrika, 103, 303–317.

Engelke S, Hitz AS (2020). “Graphical models for extremes (with discussion).” J. R. Stat. Soc. Ser. B Stat. Methodol., 82, 871–932.

See Also

Other sampling functions: rmpareto(), rmpareto_tree(), rmstable()

Examples

## A 4-dimensional HR tree model

my_tree <- igraph::graph_from_adjacency_matrix(rbind(
  c(0, 1, 0, 0),
  c(1, 0, 1, 1),
  c(0, 1, 0, 0),
  c(0, 1, 0, 0)
),
mode = "undirected"
)
n <- 10
Gamma_vec <- c(.5, 1.4, .8)
rmstable_tree(n, "HR", tree = my_tree, par = Gamma_vec)

## A 4-dimensional Dirichlet model with asymmetric edge distributions

alpha <- cbind(c(.2, 1, .5), c(1.5, .6, .8))
rmstable_tree(n, model = "dirichlet", tree = my_tree, par = alpha)

Simulate HR extremal functions

Description

Simulates the Huesler-Reiss extremal functions

Usage

simu_px_HR(n, d, idx, trend, chol_mat)

Arguments

Value

Numeric n \times d matrix. Simulated data.

Simulate Dirichlet extremal functions

Description

Simulates Dirichlet extremal functions

Usage

simu_px_dirichlet(n, d, idx, alpha)

Arguments

Value

Numeric n \times d matrix. Simulated data.

Simulate logistic extremal functions

Description

Simulates logistic extremal functions

Usage

simu_px_logistic(n, d, idx, theta)

Arguments

Value

Numeric n \times d matrix. Simulated data.

Simulate negative logistic extremal functions

Description

Simulates negative logistic extremal functions

Usage

simu_px_neglogistic(n, d, idx, theta)

Arguments

Value

Numeric n \times d matrix. Simulated data.

Simulate HR extremal functions on a tree

Description

Simulates the Huesler-Reiss extremal functions on a tree

Usage

simu_px_tree_HR(n, Gamma_vec, A)

Arguments

Value

Numeric n \times d matrix. Simulated data.

Simulate Dirichlet extremal functions on a tree

Description

Simulates Dirichlet extremal functions on a tree

Usage

simu_px_tree_dirichlet(n, alpha.start, alpha.end, A)

Arguments

Value

Numeric n \times d matrix. Simulated data.

Simulate logistic extremal functions on a tree

Description

Simulates logistic extremal functions on a tree

Usage

simu_px_tree_logistic(n, theta, A)

Arguments

Value

Numeric n \times d matrix. Simulated data.

Split graph into invariant subgraphs

Description

Split graph into invariant subgraphs

Usage

split_graph(g)

Arguments

Value

List of invariant subgraphs

Uniform margin

Description

Rescale the vector x empirically to uniform margin.

Usage

unif(x, na.rm = FALSE)

Arguments

Value

Numeric vector with entries rescaled to uniform margins