Skip to contents

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

data

Numeric \(n \times d\) matrix, where n is the number of observations and d is the dimension.

p

Numeric between 0 and 1 or NULL. If NULL (default), it is assumed that the data are already on multivariate Pareto scale. Else, p is used as the probability in the function data2mpareto() to standardize the data.

method

One of "vario", "ML", "chi". Default is method = "vario".

cens

Logical. This argument is considered only if method = "ML". If TRUE, then censored likelihood contributions are used for components below the threshold. By default, cens = FALSE.

Value

List consisting of:

graph

An igraph::graph object. The fitted minimum spanning tree.

Gamma

Numeric \(d \times d\) estimated variogram matrix \(\Gamma\) corresponding to the fitted minimum spanning tree.

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
#>      [,1] [,2] [,3] [,4]
#> [1,]  0.0  0.5  1.9  1.3
#> [2,]  0.5  0.0  1.4  0.8
#> [3,]  1.9  1.4  0.0  2.2
#> [4,]  1.3  0.8  2.2  0.0

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)