List of Figures

Apprentissage de structures dans les valeurs extrêmes en grande dimension

Maël Chiapino RESUME : Nous présentons et étudions des méthodes d'apprentissage non-supervisé de phénomènes extrêmes multivariés en grande dimension. Dans le cas où chacune des distributions marginales d'un vecteur aléatoire est à queue lourde, l'étude de son comportement dans les régions extrêmes (i.e. loin de l'origine) ne peut plus se faire via les méthodes usuelles qui supposent une moyenne et une variance finies. La théorie des valeurs extrêmes offre alors un cadre adapté à cette étude, en donnant notamment une base théorique à la réduction de dimension à travers la mesure angulaire. La thèse s'articule autour de deux grandes étapes : -Réduire la dimension du problème en trouvant un résumé de la structure de dépendance dans les régions extrêmes. Cette étape vise en particulier à trouver les sous-groupes de composantes étant susceptible de dépasser un seuil élevé de façon simultané. -Modéliser la mesure angulaire par une densité de mélange qui suit une structure de dépendance déterminée à l'avance. Ces deux étapes permettent notamment de développer des méthodes de classification non-supervisée à travers la construction d'une matrice de similarité pour les points extrêmes.

MOTS-CLEFS :

Théorie des valeurs extrêmes, apprentissage non-supervisé, réduction de dimension.

ABSTRACT :

We present and study unsupervised learning methods of multivariate extreme phenomena in high-dimension. Considering a random vector on which each marginal is heavy-tailed, the study of its behavior in extreme regions is no longer possible via usual methods that involve finite means and variances. Multivariate extreme value theory provides an adapted framework to this study. In particular it gives theoretical basis to dimension reduction through the angular measure.

KEY-WORDS :

Extreme value theory, unsupervised learning, dimension reduction, clustering.

Résumé 1.INTRODUCTION

La théorie des valeurs extrêmes répond à la nécessité de construire un modèle d'extrapolation de phénomènes hors des limites au sein desquelles ils sont habituellement observés. Étant donné un certain nombre d'observations, l'objectif est de prévoir dans quelle mesure les futures réalisations du phénomène sont susceptibles de franchir des seuils du même ordre de grandeurs que les valeurs maximales précédemment atteintes. Et plus particulièrement, de caractériser leur comportement au-delà de ces limites. Ces considérations ne prennent pleinement leur légitimité que lorsque la distribution des quantitées étudiées est à queue lourde, c'est à dire qu'elle ne concentre pas toute sa masse autour d'un point central et donc que la probabilité qu'une valeur soit extrême est non négligeable. Lorsque le processus considéré est univarié, la théorie des valeurs extrêmes peut se résumer à l'étude des valeurs maximales que prend ce processus. Le comportement asymptotique du maximum d'une telle variable aléatoire est bien connu et est entièrement décrit par une famille de distributions paramétriques, la loi de valeurs extrêmes généralisée [START_REF] De Haan | On regular variation and its application to the weak convergence of sample extremes[END_REF]). Dans un contexte multivarié, contrairement à R, l'espace ne possède plus de relation d'ordre et il n'est donc plus possible d'ordonner les points du processus de sorte à pouvoir en définir les points maximaux. Cependant, via les dépassements de seuils, le maximum composante par composante ou plus généralement la norme du vecteur aléatoire, il est possible de donner un sens au caractère extrémal de certains points.

Plusieurs difficultés propres aux extrêmes surviennent avec la dimensionnalité du problème. Les méthodes habituelles de réduction de la dimension, comme l'analyse en composantes principales, se basent sur l'étude de la matrice de covariance afin d'en extraire un modèle dont la structure est simplifiée. Dans le cas de distributions à queue lourde la covariance des composantes n'est plus définie, il n'est donc plus possible d'utiliser de tels outils. S'ajoute à cela un problème d'ordre statistique inhérent à l'objectif poursuivi d'extrapoler au-delà des limites déjà observées du phénomène. Plus précisémment, le caractère extrémal du processus que nous voulons estimer nous oblige à ne considérer, pour l'inférence statistique, qu'une partie des observations données, i.e. les plus extrêmes d'entre elles. Et cela ne peut représenter, par nature, qu'une faible proportion du jeu de données.

Ce résumé est articulé de la manière suivante. Tout d'abord dans la Section 1 1.2. THÉORIE DES VALEURS EXTRÊMES 1.2, le cadre général de la théorie des extrêmes est présenté dans le cas univarié 1.2.1 puis multivarié 1.2.2. En particulier y est introduite la mesure angulaire Φ, qui caractérise la structure de dépendance de toute loi d'extrêmes multivariés, après une certaine normalisation des distributions marginales. Une représentation parcimonieuse de cette mesure est décrite dans la Section 1.3. Sur ces bases, une methode générale pour trouver le support de Φ est détaillé dans la section 1.4. En particulier, l'heuristique CLEF (CLustering of Extreme Features) est présentée dans 1.4.2 suivi de versions alternatives de l'algorithme basées sur une série de tests statistiques 1.4.3. Une modélisation paramétrique de la mesure angulaire est enfin proposée dans la Section 1.5.

THÉORIE DES VALEURS EXTRÊMES

Théorie des valeurs extrêmes univariées

Dans le cas univarié, la théorie des valeurs extrêmes caractérise le comportement asymptotique du maximum de variables aléatoires à valeurs dans R. Soient X 1 , . . . , X n , n copies indépendantes et identiquement distribuées (i.i.d. ) d'une variable aléatoire X ∈ R de distribution F , telle que 0 < F (x) < 1 pour tout x ∈ R.

Admettons que l'on se donne pour but d'estimer de larges quantiles situés dans la queue de distribution, i.e. P[X > x p ] = p avec x p un seuil élevé. Dans le cas où aucune des données observées X i ne se trouve dans la région [x p , ∞), la seule solution est de construire un modèle d'extrapolation. L'approche classique porte donc sur l'étude du maximum M n := max{X 1 , . . . , X n } et la problématique se traduit par la limite suivante, soit x ∈ R:

P[M n ≤ x] = F n (x) -→ n→∞ 0, (1.1)
de telle sorte que la distribution limite est nécessairement dégénérée. Il est clair que les méthodes usuelles qui portent sur le comportement 'moyen' de la variable d'étude ainsi que toute approche proprement empirique ne peuvent ici être appliquées. Afin de contourner (1.1) il est nécessaire d'appliquer une certaine normalisation au maximum M n permettant l'émergence d'une limite non-dégénérée. La propriété suivante, fondatrice de la théorie des extrêmes, est en quelque sorte l'analogue du théorème central limite (TCL) et de la limite de la somme normalisée de variables aléatoires indépendantes et identiquement distribuées. Bien que l'existence d'une telle limite non-dégénérée ne soit pas garantie, à supposer qu'elle existe, sa forme est quant à elle entièrement caractérisée. C'est une différence majeure avec le TCL qui stipule l'existence de la limite pour toute variable dont les deux premiers moments sont finis. Supposons l'existence de deux suites (a n ) n≥1 et (b n ) n≥1 avec a i > 0 telles que:

lim n→∞ P M n -b n a n ≤ x = G(x), (1.2) CHAPTER 1. RÉSUMÉ 1
où G est une distribution non-dégénérée, alors cette distribution limite a nécessairement la forme suivante:

G(x) = exp -1 + ξ x -µ σ -1 ξ + .
(1.3)

On dit que F appartient au domaine d'attraction de la loi de valeurs extrêmes G, noté F ∈ DA (G). Lorsque ξ = 0, G est une loi extremum de type I (Gumbel), ξ > 0 correspond à une loi extremum de type II (Fréchet) et ξ < 0 une loi extremum de type III (Weibull).

Le comportement asymptotique du maximum de telles variables aléatoires étant entièrement caractérisé au sein d'une famille de distributions paramétriques, les problématiques liées à l'inférence de tels modèles concernent principalement la proportion de points extrémaux du jeu de données à prendre en compte.

Théorie des valeurs extrêmes multivariées

Dans un contexte multivarié, la définition d'un point extrême n'est plus aussi naturelle. En effet, dès lors que d ≥ 2, R d n'est plus un ensemble ordonné, il n'est plus possible de classer les points du plus au moins grand par une relation d'ordre. Néanmoins, considérant x ∈ R d il est approprié de définir un tel point 'extrême' lorsque max j∈{1,...,d} x j > t ou plus généralement x > t pour t > 0 un seuil élevé et • une norme quelconque. De façon analogue au cas univarié, la théorie des valeurs extrêmes multivariées se fonde sur la distribution limite du maximum composante par composante normalisé.

Soient X 1 , . . . , X n , n variables aléatoires i.i.d. dans R d de distribution F . Les distributions marginales de F sont notées F j pour j = 1, . . . , d. Soit M n := ( max i=1,...,n X i,1 , . . . , max i=1,...,n X i,d ) le n-ème maximum composante par composante. Supposons l'existence de deux suites (a n ) n≥1 dans (0, ∞) d et (b n ) n≥1 dans R d , telles que:

lim n→∞ P
M n,j -b n,j a n,j ≤ x j , j = 1, . . . , d = G 0 (x), (1.4) où G 0 est une distribution multivariée dont les distributions marginales G 0,j , j = 1, . . . , d sont non-dégénérées. Alors G 0 est caractérisée par l'expression suivante [START_REF] Resnick | Extreme Values, Regular Variation, and Point Processes[END_REF]), pour un certain x 0 ∈ R d :

G 0 (x) =     
exp -µ 0 [x 0 , x] c , pour x ≥ x 0 , 0 sinon, (1.5) où µ 0 est la dite mesure exposant, une mesure de Radon définie sur

[x 0 , ∞] \ {x 0 } et où [x 0 , x] c = [x 0 , ∞] \ [x 0 , x].
Chaque marginale G 0,j est une loi de valeurs extrêmes univariée et donc pour tout j ∈ {1, . . . , d}, F j ∈ DA(G 0,j ). Cependant, 1 1.2. THÉORIE DES VALEURS EXTRÊMES la mesure exposant µ 0 , et donc G 0 , ne peuvent être caractérisées par une famille de distributions paramétriques. C'est une différence majeure avec le cas univarié. Soit X = (X 1 , . . . , X d ) une variable aléatoire à valeurs dans R d de distribution F . L'étude de la queue de distribution de F se sépare en deux parties distinctes. D'une part, la caractérisation des distributions marginales, qui relève du cas univarié. D'autre part, l'étude de la structure de dépendance, qui constitue l'essentiel du problème. De façon similaire à l'étude des copules, une étape préliminaire naturelle est de standardiser chaque marginale de sorte à ce qu'elles soient de distributions identiques. Bien que le choix de la distribution soit arbitraire, il est commun de transformer chaque marginale en Pareto unitaire [START_REF] Resnick | Extreme Values, Regular Variation, and Point Processes[END_REF]); une standardisation que l'on obtient par le procédé suivant:

V j := (1 -F j (X j )) -1 , pour j = 1, . . . , d.

(1.6)

Remarque 1. Soit X 1 , . . . , X n , n copies i.i.d. de X. Dans la pratique, les distributions marginales F j n'étant pas connues, la transformation (1.6) est effectuée par le biais des distributions empiriques Fj (x) := 1 n n i=1 1 {X i,j <x} où x ∈ [x 0,j , ∞). Notons rang(X i ,j ) = n -n i=1 1 {X i,j <X i ,j } . La transformation devient donc, pour tout j ∈ {1, . . . , d}: (X i,j ) .

Vi,j = (1 -Fj (X i,j )) -1 ) = n rang
Nous travaillons désormais avec le vecteur transformé V = (V 1 , . . . , V d ), à valeurs dans (1, ∞) d , et de distributions marginales Pareto unitaire. L'existence de la distribution asymptotique G 0 dans (1.4) peut alors se reformuler de façon approprié à travers le concept de variation régulière. La variable aléatoire V est dite à variation régulière d'indice -1 sur (1, ∞) d et de mesure limite µ définie sur [0, ∞] d \ {0} si (voir par ex. [START_REF] Resnick | Extreme values, regular variation and point processes[END_REF]):

tP t -1 V ∈ A ---→ t→∞ µ(A),
(1.7) pour tout borélien A dans [0, ∞] d \ {0} tel que µ(∂A) = 0 et 0 / ∈ ∂A. Autrement dit, µ est la mesure exposant associée à la loi de valeurs extrêmes G de V de telle sorte que pour tout v ∈ [0, ∞) d \ {0}:

G(v) = exp[-µ([0, v] c )].
(1.8)

Le lien entre les distributions asymptotiques G 0 et G associées respectivement à X et V se traduit par:

G(v) = G 0 (G ← 0,1 (e -1/v 1 ), . . . , G ← 0,d (e -1/v d )),
où G ← 0,i est l'inverse généralisé de G 0,i .

CHAPTER 1. RÉSUMÉ 1 Remarque 2. Soit t > 0 un seuil élevé, sachant que pour une coordonnée particulière j 0 ∈ {1, . . . , d}, V j 0 est extrême, i.e. V j 0 > t, la fonction x → µ [0, x] c est approximativement, pour un v ∈ [0, ∞) d \ {0} donné, la probabilité conditionnelle suivante:

P V 1 > tv 1 or . . . or V d > tv d |V j 0 > t = P V ∈ t[0, v] c |V j 0 > t = tP V ∈ t[0, v] c ≈ µ [0, v] c = -log G(v)
L'intérêt de la transformation (1.6) ainsi que de l'hypothèse de variation régulière (1.7) est bien l'accent porté sur la structure de dépendance asymptotique caractérisée par la mesure µ. Il convient de rappeler qu'il n'existe pas de famille paramétrique englobant complètement la classe de telles mesures. Néanmoins, en conséquence de (1.6) et du caractère max-stable de toute loi de valeurs extrêmes généralisée [START_REF] Resnick | Extreme Values, Regular Variation, and Point Processes[END_REF]), la fonction µ est homogène de degré -1:

µ(tA) = t -1 µ(A),
(1.9) pour tout t > 0 et pour tout borélien A tel que 0 / ∈ ∂A. On peut déduire de (1.9) que les distributions marginales de G sont de loi Fréchet unitaire.

Remarque 3. Il suit de la propriété d'homogénéité (1.9) que les régions de [0, ∞] d \ {0}, minimales pour l'inclusion, sur lesquelles µ possède une masse non nulle sont nécessairement des cônes, i.e. C ⊂ [0, ∞] d \ {0} tel que x ∈ C ⇒ tx ∈ C pour tout t > 0.

Cette propriété est essentielle afin de mieux caractériser la classe des mesures exponents suite à la standardisation des marginales. En effet, (1.9) induit une décomposition pseudo-polaire de µ. Pour tout v ∈ [0, ∞) d \ {0}, considérons l'application bijective:

T : v → ( v , v v ),
où • est une norme quelconque sur R d . Il est alors possible de définir une mesure Φ sur le quadrant positif de la sphère unité S d = {x ≥ 0 : x = 1}, la dite mesure angulaire, caractérisant la structure de dépendance de V :

Φ(A) := µ({v : v > 1, v v ∈ A}),
(1.10) pour tout A ⊂ S d . La propriété d'homogénéité de µ (1.9) se traduit alors par:

µ({v : T (v) ∈ (t, ∞) × A}) = µ({v : v > t, v v ∈ A}) = t -1 Φ(A).

1.2. THÉORIE DES VALEURS EXTRÊMES

Finalement, cette dernière équation ainsi que la bijectivité de T impliquent la factorisation de µ en deux mesures indépendantes (de Haan and [START_REF] De Haan | Limit theory for multivariate sample extremes[END_REF]): µ • T -1 (dr, dw) = r -2 drΦ(dw).

(1.11)

Remarque 4. Lorsque V est grand, la mesure angulaire Φ jauge la part relative que chacune des coordonnées V j prend dans le caractère extrême de V , autrement dit, la direction dans laquelle V sera susceptible d'être extrême. Posons R = V et W = R -1 V , une réécriture de la limite (1.7) en fonction de Φ est:

P (W ∈ A, R > rt | R > t) ---→ t→∞ r -1 Φ(S d ) -1 Φ(A),
(1.12) pour tout borélien A de S d tel que Φ(∂A) = 0 et r > 1. En d'autres termes, lorsque la composante radiale R est grande, R et le pseudo-angle W sont approximativement indépendants. La distribution de W est alors approximativement (et à une normalisation près) la mesure angulaire et R suit approximativement une loi de Pareto unitaire. Une étude plus précise de ces caractéristiques est développée dans la section 1.3. Tout comme pour la mesure exposant, il n'existe pas de famille paramétrique décrivant entièrement l'espace des mesures angulaires. La seule condition sur Φ, pour être une mesure angulaire valide, est la contrainte des moments: S d w j Φ(dw) = 1, ∀j ∈ {1, . . . , d}.

(1.13) C'est une conséquence de la transformation de chaque marginale en Pareto unitaire (1.6). En effet nous avons, pour tout j ∈ {1, . . . , d} et pour tout x j ∈ (0, ∞):

G(∞, . . . , x j , . . . , ∞) = e -1

x j .

De plus, en utilisant la décomposition (1.11), nous pouvons réécrire x → µ([0, x] c ) en fonction de la mesure angulaire. Soit x dans [0, ∞] d \ {0}:

µ([0, x] c ) = 1 {∃j:u j >x j } µ(du) = S d ∞ 0 1 {r>min j x j w j } r -2 drΦ(w) = S d (min j x j w j ) -1 Φ(w) = S d max j w j x j Φ(w).
Ainsi, pour x = (∞, . . . , x j , . . . , ∞), avec x j ∈ (0, ∞) nous avons:

1 x j = -log G(x) = µ([0, x] c ) = S d w j
x j Φ(dw).

Travaux antérieurs et problèmes

D'un point de vue statistique, de nombreux modèles paramétriques ont été proposés pour la loi de valeurs extrêmes G à travers la modélisation de la fonction x → µ[0, x] c , i.e. G(x) = exp [-µ[0, x] c ] (voir par ex. [START_REF] Coles | Modeling extreme multivariate events[END_REF]). Les deux stratégies principales pour l'inférence d'un tel modèle sont les méthodes de maximum composante par composante et des excès de seuils, l'étude porte alors respectivement sur l'estimation de G et sur l'estimation de x → µ[0, x] c . Par le biais du maximum composante par composante, l'estimation bute sur le calcul de la vraisemblance lorsque la dimension augmente (d > 10). En effet, il est nécessaire de calculer:

∂ d ∂x 1 . . . ∂x d e -µ[0,x] c
qui est une somme dont le nombre de termes explose en fonction de d, (voir la revue [START_REF] Huser | Likelihood estimators for multivariate extremes[END_REF]). Des simplifications ont été proposées (par ex. [START_REF] Wadsworth | On the occurrence times of componentwise maxima and bias in likelihood inference for multivariate max-stable distributions[END_REF], [START_REF] Stephenson | Exploiting occurrence times in likelihood inference for componentwise maxima[END_REF]) mais le calcul reste infaisable lorsque d > 10.

Que ce soit pour le maximum composante par composante ou pour les excès de seuils, la plupart de ces approches ne considère que le cas de dépendance totale entre les variables, c'est-à-dire que le support de la mesure exposant µ est confiné à l'intérieur de (0, ∞) d . L'estimation non-paramétrique n'échappe pas, quant à elle, au problème plus général du fléau de la dimension. Vient s'ajouter un problème inhérent aux extrêmes: le faible nombre de points disponibles pour effectuer l'estimation.

Estimation de la mesure angulaire Diverses approches paramétriques et non-paramétriques ont été proposées pour l'estimation de la mesure angulaire. (Boldi and Davison (2007b)) et [START_REF] Sabourin | Bayesian model averaging for multivariate extremes[END_REF]) considèrent la classe des mélanges de Dirichlet sur le simplexe. Cette classe de modèles est dense (au sens faible) dans l'espace des mesures angulaires. D'autres modèles non-paramétriques 1 1.3. MESURE ANGULAIRE PARCIMONIEUSE ont été proposés [START_REF] Guillotte | Non-parametric bayesian inference on bivariate extremes[END_REF]), Fougeres et al. (2013)). Cependant les expérimentations ne vont pas au-delà de dimensions moyennes (d ≈ 5). De récentes approches tentent de dépasser ce cap à travers des formes de clustering [START_REF] Chautru | Dimension reduction in multivariate extreme value analysis[END_REF]) et un algorithme permettant de retrouver le support de la mesure angulaire en grande dimension (d ≈ 50) est proposé dans (Goix et al. (2015a), Goix et al. (2016a)).

MESURE ANGULAIRE PARCIMONIEUSE

En grande dimension, lorsque d ≈ 100, il est raisonnable de supposer qu'un phénomène extrême, décrit par une variable aléatoire V = (V 1 , . . . , V d ), ne soit pas dû à l'ensemble de ses composantes simultanément, i.e. V j > t, ∀j ∈ {1, . . . , d} où t > 0 est un seuil élevé. Plus précisément, l'hypothèse fondamentale que nous faisons est qu'un tel phénomène ne prend son caractère extrémal que par le biais de certains sous-groupes précis de composantes. De telle sorte que lorsque V est grand, il existe un sous-groupe α ⊂ {1, . . . , d} tel que

V ≈ V α , où |α| d et V α = (V j ) j∈α .
Nous entendons mettre en avant le caractère multimodal du comportement extrême de V et donc la pluralité de tels α, qui identifient les sous-groupes de composantes dépendantes lors des différentes réalisations de V au delà d'un large seuil. Il est naturel de supposer que le nombre de sous-groupes de composantes pouvant vérifier la propriété précédente est bien moindre que le nombre total de sous-ensembles de {1, . . . , d}, à savoir 2 d -1. De façon plus formelle, soit Φ la mesure angulaire associée à V , la structure de dépendance asymptotique de V est caractérisée par la répartition de la masse de Φ sur l'orthant positif de la sphère unité. Dans le cas d'une structure de dépendance parcimonieuse, la masse de Φ n'est répartie que sur l'intérieur de certains bords de S d , correspondant aux différents groupes de composantes asymptotiquement dépendantes. En effet, supposons qu'au moins deux composantes V i et V j soient asymptotiquement indépendantes. Cela signifie:

tP (V i > t, V j > t) = tP t -1 V ∈ {v : v i > 1, v j > 1} ---→ t→∞ 0. Or pour tout w ∈ Sd , il existe v ∈ [0, ∞) d \ {0} tel que v v = w, v i > 1 et v j > 1, donc: Φ( Sd ) = µ {v : v > 1, v v ∈ Sd } ≤ µ {v : v i > 1, v j > 1} = 0.
À travers cet exemple simple, nous voyons que l'indépendence asymptotique entre certaines des composantes implique que la masse de la mesure angulaire est confinée sur les bords de S d . Nous dénommons S α l'intérieur du bord de S d associée aux coordonnées α ⊂ {1, . . . , d}: (1.14) et M l'ensemble des sous-groupes de composantes de {1, . . . , d} pour lesquels la mesure angulaire attribue une masse non nulle au sous-espace associé:

S α := {w ∈ S d : ∀j ∈ α w j > 0, ∀j / ∈ α w j = 0},
M := {α ⊂ {1, . . . , d} : Φ(S α ) > 0}.
(1.15) L'hypothèse de parcimonie se traduit donc par la double inégalité suivante:

   max α∈M |α| d, |M| 2 d -1.
(1.16)

Problème de l'estimation de M.

Soient V la variable aléatoire à valeurs dans (1, ∞) d de marginales Pareto unitaire et de mesure angulaire Φ, et (R, W ) la décomposition pseudo-polaire associée où

R = V et W = V V .
Il est primordial de rappeler que l'équation limite (1.12), où Φ est la mesure limite de W pour R > t avec t → ∞, n'est vérifiée que sur les boréliens A de S d tels que Φ(∂A) = 0. En particulier, cette limite n'est valide que pour les sous-espaces de S d de mesure de Lebesgue non nulle. Cela ne concerne donc pas les sous-espaces S α chargés, puisque dans R d , S α ⊂ ∂S α , et donc pour α ∈ M:

0 < Φ(S α ) ≤ Φ(∂S α ).
(1.17) [START_REF] Coles | Modeling extreme multivariate events[END_REF]), qui avec le modèle logistique asymmétrique, constitue une classe de modèles paramétriques très flexible en dimension d > 1 quelconque. Pour tout v dans (0, ∞) d :

µ lgtc,θ ([0, v] c ) = d j=1 v -1 θ j θ , (1.18)
où θ > 0. L'indépendance asymptotique est atteinte pour θ = 1 et dans ce cas, la mesure angulaire associée place toute sa masse sur les axes, i.e. Φ(e j ) = 1 pour j ∈ {1, . . . , d} où e j ∈ R d tel que e j i = 1 pour i = j et e j i = 0 sinon. Pourtant, il est clair que tout point distribué selon la loi de valeurs extrêmes G(v) = expµ lgtc,1 ([0, v] c ) est presque sûrement à valeurs dans (0, ∞) d et donc le pseudo-angle associé est, quant à lui, presque sûrement à valeurs dans l'intérieur de S d .

1.4. ESTIMATION DU SUPPORT DE LA MESURE ANGULAIRE

Remarque 5. Le lien entre la loi de valeurs extrêmes G et la mesure exposant µ n'existe qu'à travers l'équation (1.8), ainsi, le support de µ (et donc de la mesure angulaire Φ) n'est pas le support de G. Autrement dit, la mesure angulaire ne représente que l'angle asymptotique de V , et le fait que le support de Φ ne soit inclus que dans certains sous-espaces S α ne signifie pas que le support de G est inclus dans des hyperplans de R d .

L'avantage de modéliser Φ plutôt que G est que cela permet d'appliquer naturellement un modèle parcimonieux à la structure de dépendance. Le revers, est que le support de la distribution que l'on modélise n'est jamais atteint dans l'espace au sein duquel les points sont observés (Remarque 5). Nous considérons par la suite deux voies principales afin de répondre à ce problème. La première, correspondant aux Chapitres 3 et 4, revient à considérer des sous-espaces proches des bords de S d mais de mesure de Lebesgue non nulle, et d'estimer la mesure angulaire (ou la mesure exposant) de façon non-paramétrique par le biais d'une mesure de comptage sur ces sous-espaces. La deuxième approche est d'appliquer un modèle de mélange paramétrique pour Φ, où chaque composante du mélange correspond à un groupe de composantes asymptotiquement dépendantes, puis de projeter V sur les sous-espaces correspondants, et de considérer le résidu comme un bruit orthogonal. Cette dernière approche est développée dans le Chapitre (5).

ESTIMATION DU SUPPORT DE LA MESURE ANGULAIRE

D'un point de vue pratique, lorsqu'il s'agit simplement d'estimer le support de la mesure angulaire sur S d , et non de modélisation à proprement parler, il est plus naturel de se placer directement dans l'espace de définition de la mesure exposant µ, i.e. [0, ∞] d \ {0}. Par définition de Φ à l'équation (1.10), la contrepartie d'un sous-espace S α où α ⊂ {1, . . . , d}, en terme de sous-espace de [0, ∞] d \ {0}, est le cône tronqué:

C α = {x ≥ 0 : x > 1, ∀j ∈ α x j > 0, ∀j / ∈ α x j = 0}. (1.19)
En effet, pour tout α dans {1, . . . , d}:

Φ(S α ) = µ(v : v > 1, v v ∈ S α ) = µ(C α ).

Inférence

Considérons à nouveau la variable aléatoire V , à variation régulière dans (1, ∞) d de mesure exposant µ et de marginales Pareto unitaires. À l'instar de Φ pour le pseudo-angle V V à valeurs dans l'intérieur de S d , V étant à valeurs (1, ∞) d , pour tout t > 0 et α ⊂ {1, . . . , d} un sous-ensemble strict:

P V t ∈ C α = 0. (1.20)
De plus, les cônes C α étant d'intérieur vide nous avons µ(∂C α ) > 0, de façon analogue à (1.17) pour les sous-espaces S α .

Pour estimer µ(C α ) à partir d'un échantillon i.i.d. de copies de V , il est nécessaire d'élargir les cônes tronqués C α dans l'intérieur de [0, ∞] d \ {0}. L'idée est de trouver vers quels bords de l'orthant positif (i.e. C α ), V se concentre lorsque V est grand. Considérons donc les cônes élargis suivant, pour > 0 et α ⊂ {1, . . . , d}:

R α = {x ≥ 0 : x > 1, ∀j ∈ α x j > , ∀j / ∈ α x j ≤ }.
(1.21) De façon plus formelle, et conséquemment à la propriété de passage à la limite monotone décroissante de toute mesure:

µ(C α ) = µ(∩ >0, ∈Q R α ) = lim →0 µ(R α ).
Afin de pouvoir appliquer la limite (1.7) aux ensembles R α , ces derniers doivent vérifier µ(∂R α ) = 0. Le bord de R α dans [0, ∞] d \ {0} est composé d'hyperplans parallèles au axes et disjoints pour des valeurs différentes de ∈ (0, 1). Il n'existe de ce fait qu'un nombre au plus dénombrable de tel que µ(∂R α ) > 0, sinon µ serait infinie sur une région compacte de [0, ∞] d \ {0}. Ainsi nous pouvons choisir un seuil > 0, arbitrairement petit et tel que µ(∂R α ) = 0, pour estimer la masse des cônes C α via la limite suivante:

µ(C α ) = lim →0 lim t→∞ tP (V ∈ tR α )
De récents travaux (Goix et al. (2016a)) proposent un algorithme nommé DAMEX permettant d'estimer la masse µ(C α ) par une méthode de comptage sur les sous-espaces tR α pour > 0 petit et un large seuil t > 0. Plus précisément, soient V 1 , . . . , V n , n copies i.i.d. de la variable V , alors l'estimateur de µ(C α ) est: (1.22) où l'on a remplacé le seuil t par n k , avec k ≤ n un ordre de grandeur du nombre de points extrêmes parmi l'échantillon considéré. Un seuil minimal µ min > 0 est ensuite choisit de telle sorte que l'on décide:

μα = 1 k n i=1 1 {V i ∈ n k R α } ,
µ(C α ) > 0 si μα ≥ µ min .
L'algorithme est testé sur des données réelles (données de directions de vagues dans la mer du nord fournies par Shell) et parvient à réduire nettement la dimension du problème en rassemblant l'essentiel de la masse de µ dans un nombre limité de sous-cônes C α de dimension moyenne. De la même façon, cette méthode montre de bons résultats sur des données simulées à partir d'un modèle paramétrique pour G, une version parcimonieuse du modèle logistique asymmétrique (Tawn (1990a)):

G(z) = exp - α∈M j∈α (|A(j)|z j ) -1/θα θα , (1.23)
où A(j) = {α ∈ M : j ∈ α} et où θ α > 0 est un paramètre de dépendance tel que θ α = 1 correspond à l'indépendance asymptotique et θ α ↓ 0 à la dépendance totale. Néanmoins, les expériences sont faites pour θ α = 0.1 (dépendance forte) et lorsque l'on augmente la dimension d, θ α et le nombre |M| de sous-groupes de composantes dépendantes, la méthode ne parvient plus à retrouver M. De façon intrinsèque, en raison de la partition de l'espace en un nombre considérable de sous-cônes -(2 d -1), le modèle montre une forte sensibilité à la moindre variation des sous-groupes de composantes impliquées dans le comportement extrême de V . En effet, soit α 0 tel que µ(C α 0 ) > 0, alors il est possible que pour un certain nombre de composantes j ∈ {1, . . . , d} \ α 0 les grandeurs estimées μα 0 ∪{j} soient non-négligeables:

μα 0 ∪{j} = 1 k n i=1 1 {V i ∈ n k R α 0 ∪{j} } ≥ µ min . (1.24)
Il en résulte un éparpillement de la masse au sein de sous-groupes de cônes proches les uns des autres, rendant impossible l'émergence d'un modèle simplifié de la structure de dépendance. Dans le cas d'une dépendance asymptotique moyenne, voire faible (i.e. θ α ∈ [0.5, 1)) entre les composantes ou plus généralement de la présence de bruit dans le jeu de données, il est nécessaire de contruire une méthode plus robuste pour estimer le support de µ.

Les figures suivantes (1.2) montrent la répartition des points extrêmes dans les cônes n k R α où k est choisit de telle sorte que le nombre de points extrêmes

n i=1 1 { V i > n k } représente 5% des données.
Pour chaque figure le coefficient de dépendance θ α est fixé à une valeur particulière. Les points (V 1 , . . . , V n ), avec n = 1e5, sont générés à partir d'une distribution logistique asymmétrique en dimension d = 50, où l'on a généré au préalable 50 sous-groupes de composantes de façon aléatoire pour former M. L'axe des abscisses représente les sous-groupes de composantes α pour lesquels l'estimateur μα est non nul, ordonné de façon décroissante. L'axe des ordonnés représente la valeurs des estimateurs μα . En rouge sont représentés les éléments de M au nombre de 50. Ainsi, afin de retrouver M, le seuil minimal µ min > 0 doit être tel que:

   μα ≥ µ min , pour tout α ∈ M μα < µ min , sinon.
(1.25)

Pour θ α = 0.1 (1.2 en haut à gauche) l'écart entre les 'vrais' sous-groupes de composantes (i.e. appartenant à M) et les autres est net, et il est possible de choisir un seuil µ min tel que (1.25) soit vérifiée. Lorsque θ α augmente, correspondant à une dépendence asymptotique plus faible, la répartition de la masse des points extrêmes s'éparpille et pour θ α = 0.7, 0.9 (1.2 en bas gauche (resp. droite)), il n'est plus possible de retrouver M.

.2 Algorithme CLustering Extreme Feature

Plutôt que de considérer la partition en sous-cônes C α qui conduit à différencier des événements extrêmes proches, l'idée est de rechercher seulement les sous-groupes α ⊂ {1, . . . , d} de taille |α| ≥ 2 tels que pour tout j ∈ α, V j soit grand, quelles que soient les contreparties V j pour j ∈ {1, . . . , d} \ α. Nous définissons les sous-espaces imbriqués suivant:

Γ α = {x : ∀j ∈ α x j > 1}.
(1.26)

De façon analogue au sous-cônes R α , nous avons µ(∂Γ α ) = 0, et il est possible d'estimer leur masse en utilisant la limite (1.7). L'avantage de cette approche est qu'elle permet d'éviter la situation (1.24) d'éparpillement de la masse lors de l'estimation du support de µ. Pour tout α ⊂ {1, . . . , d}, Γ α correspond à l'évènement:

{ Toutes les variables (V j ) j∈α sont grandes (quelles que soient les variables complémentaires (V j ) j / ∈α ). } Soit t > 0 un seuil élevé, chaque rectangle tΓ α rassemble les points extrêmes possédant de large coordonnées en α indépendamment de leurs coordonnées complémentaires {1, . . . , d} \ α. De cette façon, les points appartenant à tΓ α , appartiennent nécessairement aux sous-rectangles tΓ β pour tout β ⊂ α. L'objectif est alors de rechercher, de façon incrémentale, les sous-groupes α les plus grands tel que le nombre de points dans le rectangle tΓ α soit non-négligeable.

En termes d'estimation du support de la mesure exposant µ, ce cadre est équivalent à la partition de l'espace en cônes C α lorsque l'on considère les sous-groupes α ⊂ {1, . . . , d} maximaux pour l'inclusion tels que µ(C α ) > 0. En effet, admettons que µ donne une masse non nulle au cône C α . Suivant la remarque 3, la région C ⊂ C α , minimale pour l'inclusion, où se concentre la masse de µ, est elle-même un cône. Or tout cône inclus dans C α possède une intersection non vide avec Γ α , dès lors

µ(C α ) > 0 implique µ(Γ α ) > 0. À l'inverse, supposons que µ(Γ α ) > 0, alors si pour tout β tel que α ⊂ β, µ(Γ β ) = 0, nous avons nécessairement µ(C α ) > 0. Par contraposition, supposons µ(Γ α ) > 0 et µ(C α ) = 0, alors il existe j C ∈ {1, . . . , d} \ α et un cône C inclu dans Γ α∪{j C } tel que µ(C) > 0, autrement dit, il existe β tel que α ⊂ β et µ(Γ β ) > 0.
De ce fait, les sous-groupes α ⊂ {1, . . . , d} maximaux pour l'inclusion tels que µ(C α ) > 0, sont les sous-groupes maximaux tels que µ(Γ α ) > 0. Soit M l'ensemble des sous-groupes α ⊂ {1, . . . , d} maximaux tels que µ(Γ α ) > 0, alors nous avons la propriété suivante [START_REF] Chiapino | Feature clustering for extreme events analysis, with application to extreme stream-flow data[END_REF]), pour tout α ⊂ {1, . . . , d}: α ∈ M ⇔ α est maximal dans M.

(1.27) [START_REF] Agrawal | Fast algorithms for mining association rules[END_REF]). L'algorithme permet en particulier de réduire l'exploration totale des sous-groupes de T par le biais d'une exploration incrémentale. Appliquée à notre cadre, la procédure parcourt les sous-groupes α ⊂ {1, . . . , d}, |α| ≥ 2 de taille croissante en ne gardant que ceux qui vérifient un certain critère de décision C(α) équivalent à µ(Γ α ) > 0. La sortie de l'algorithme est donc:

M 0 := α ⊂ {1, . . . , d} : |α| ≥ 2, C(α) ,
(1.28) dont on ne garde que les éléments maximaux pour l'inclusion afin d'obtenir M.

Cette procédure peut être visualisée à travers le diagramme de Hasse (1.3), que l'on explore de haut en bas, en coupant toutes les branches partant d'un noeud lorsque celui-ci ne vérifie pas la propriété C. Cette méthode réduit drastiquement le nombre de sous-groupes qu'il aurait potentiellement fallu tester. (1.29) À l'étape s, étant donné tous les sous-groupes α de taille s tels que C(α) est vérifiée, l'algorithme construit les sous-groupes β ⊂ {1, . . . , d} de taille s + 1 susceptibles de vérifier C(β). Soit M s 0 = α ⊂ {1, . . . , d} : |α| = s et C(α) l'ensemble des sous-groupes de composantes de taille s et vérifiant le critère C. L'ensemble des sous-groupes candidats de taille s + 1 pour s ≥ 3 est donc:

A s = β ⊂ {1, . . . , d} : |β| = s et β \ {j} ∈ M s-1 0 pour tout j ∈ β .
(1.30)

γ α = 1 k n i=1 1 {V i ∈ n k Γα} ≥ γ min . (1.31)
Cependant, du fait de l'imbrication des rectangles Γ α entre eux, la valeur de tP[V ∈ tΓ α ] ≈ µ(Γ α ), pour t grand, décroit nécessairement pour des tailles croissantes de α. En effet, soient

α 1 ⊂ α 2 ⊂ {1, . . . , d}, alors Γ α 2 ⊂ Γ α 1 et donc tP[V ∈ tΓ α 2 ] ≤ tP[V ∈ tΓ α 1 ] pour tout t > 0.
Le choix du seuil γ min > 0 devrait donc dépendre de la taille des sous-groupes α considérés.

De ce fait, un résumé numérique alternatif du degré de dépendance des composantes α ⊂ {1, . . . , d}, dénommé κ α , est proposé:

κ α := µ(Γ α ) µ( β⊂α |β|=|α|-1 Γ β )
.

(1.32)

La quantité κ α ne décroit pas nécessairement avec la taille des α, et κ α > 0 si et seulement si µ(Γ α ) > 0. Il est donc possible de choisir un seuil κ min > 0, indépendant de la taille de α, afin d'avoir un critère {κ α ≥ κ min } pour décider de {κ α > 0}. L'idée est de faire dépendre l'acceptation d'un sous-groupe α dans M 0 en fonction des sous-groupes de taille inférieure le composant: • Générer les sous-groupes candidats de taille s:

κ α = P[ Tous les V j avec j ∈ α sont
A s = {α ⊂ {1, . . . , d} : |α| = s et α \ j ∈ M s-1 0 pour tout j ∈ α}.
• Retenir les sous-groupes vérifiant le critère, M s 0 = α ∈ A s : κα > κ min .

• Si M s 0 = ∅, S = s.

Sortie: M 0 = ∅ si S = 1 et M 0 = S s=2 M s 0 si S ≥ 2. PHASE 2: Ne retenir que les α maximaux. Si S = 1, alors M = ∅. Sinon:

Initialisation: M ← M s 0 . pour s = (S -1) : 2, pour α ∈ M s 0 , S'il n'existe pas de β ∈ M tel que α ⊂ β, alors M ← M ∪ {α}. Sortie: M 1.4.2.

Résultats

Dans la pratique, le choix du résumé du degré de dépendance κ α est d'autant plus justifié lorsque du bruit est observé sur les composantes. Afin de reproduire un jeu de données pouvant être assimilé à un jeu de données réelles bruité, nous reprenons le modèle logistique décrit dans la section 1.4.1, pour d = 20, |M| = 15 et θ α = 0.5, sur lequel nous perturbons artificiellement la structure de dépendance. Plus précisément, chacun des points V i ∈ (0, ∞) d , i = 1, . . . , n, du jeu de données est généré par le biais d'une version des sous-groupes de composantes dépendantes M i , légèrement perturbé par rapport à M. Pour i = 1, . . . , n, chaque sous-groupe de composantes α i k dans M i est modifié d'une composante par rapport à α k , de telle sorte qu'il existe j 0 ∈ {1, . . . , d} tel que:

{j 0 ∈ α k et j 0 / ∈ α i k } ou {j 0 / ∈ α k et j 0 ∈ α i k }. (1.35)
L'objectif est de retrouver M, à partir de ce jeu de données perturbé. Il est donc nécessaire, à chaque étape s = 2, . . . , S, avec S = max α∈M |α|, de discerner les sous-groupes de composantes suivants:

κα ≥ κ min pour tout α ∈ M s 0 κα < κ min pour tout α ∈ A s \ M s 0 .
Nous dénommons M F = ∪ s=1,...,S

A s \ M 0 l'ensemble des sous-groupes de composantes générés comme candidats mais n'appartenant pas à M 0 . Notons que M 0 = ∪ S s=2 M s 0 . De cette façon, la condition permettant l'existence d'un seuil κ min > 0 tel que M = M est:

max α∈M F κα < min α∈M 0 κα . ( 1

Coefficient de dépendance de queue

En vue de construire une série de tests statistiques qui ne reposent pas sur l'utilisation d'un seuil arbitraire, il est nécessaire d'introduire un nouveau critère, autre que {κ α > 0}, en prenant en compte la variance, afin de décider de l'hypothèse µ(Γ α ) > 0. En effet, lorsque µ(Γ α ) = 0 la distribution limite des statistiques √ k( κ α -κ α ) est dégénérée en zéro. De ce fait, nous n'avons pas de contrôle sur le niveau asymptotique des tests basés sur ces statistiques sous l'hypothèse H 0 : κ α = 0. Considérons le coefficient de dépendance de queue de distribution η α ∈ (0, 1], introduit dans le cas bivarié dans [START_REF] Ledford | Statistics for near independence in multivariate extreme values[END_REF]) et étendu en dimension d ≥ 3 dans (De Haan and Zhou (2011)) et [START_REF] Eastoe | Modelling the distribution of the cluster maxima of exceedances of subasymptotic thresholds[END_REF]). L'hypothèse

P[V ∈ tΓ α ] = t -1/ηα L α (t).
(1.37)

Supposons que la limite lim t→∞ tP[V ∈ tΓ α ] = µ(Γ α ) existe et que l'hypothèse précédente (1.37) est valide, alors µ(Γ α ) > 0 implique η α = 1. À l'inverse, supposons (1.37) et lim inf t→∞ L α (t) > 0 alors η α = 1 implique µ(Γ α ) > 0. L'hypothèse nulle µ(Γ α ) > 0 correspond donc à l'hypothèse η α = 1, à une condition peu coûteuse sur L α près. Dès lors, si η α = 1 la limite √ k( η α -η α ) est non-dégénérée et il est ainsi possible de contrôler le niveau asymptotique du test associé.

Un nouveau critère C(α) pour l'algorithme CLEF est construit à partir de l'estimateur du coefficient de dépendance de queue η α . Cette serie de tests statistiques est développée dans la partie 4 pour différents estimateurs non-paramétriques de η α , à savoir une extension multivariée de l'estimateur de Peng [START_REF] Peng | Estimation of the coefficient of tail dependence in bivariate extremes[END_REF]) et l'estimateur de Hill [START_REF] Draisma | Tail dependence in independence[END_REF], [START_REF] Draisma | Bivariate tail estimation: dependence in asymptotic independence[END_REF]). Une version non-dégénérée du test H 0 : κ α = 0 est aussi développée, en considérant la variance de κ α ainsi qu'un seuil κ min > 0, donnant la nouvelle hypothèse H 0 : κ α > κ min .

En reprenant la méthode employée pour comparer les résultats des critères {κ α ≥ κ min } et {γ α ≥ γ min } sur les données artificiellement bruitées (1.4.2.3), nous testons l'estimateur de Hill par le biais d'un nouveau critère d'acceptation de sous-groupe de composantes α ⊂ {1, . . . , d}, dont nous contrôlons le niveau asymptotique: où q 1-δ est le (1-δ)-quantile de la loi normale centrée réduite, σα la variance asymptotique de l'estimateur et k un ordre de grandeur du nombre de points extrêmes. La figure 1.9 montre les résultats du critère pour les mêmes paramètres initiaux d = 20, |M| = 15 et θ α = 0.5 que pour les figures 1.4 et 1.5 et pour une très faible valeur δ = 1e -7. La valeur mesurée en ordonnée est ηα -(1 -q 1-δ σα √ k ) et est strictement positive pour tous les sous-groupes devant être acceptés (en rouge) et strictement négative pour tous les sous-groupes devant être rejetés (en bleu).

ηα > 1 -q 1-δ σα √ k

MODÈLE PARAMÉTRIQUE POUR LA MESURE ANGULAIRE

Cette dernière partie est motivée par un problème de détection d'anomalies sur des données aéronautiques fournit par Airbus. L'approche générale de la détection d'anomalies est de construire un modèle décrivant le comportement 'normal' d'un phénomène puis de définir comme 'anormaux' des évènements en fonction de l'écart qu'ils ont avec ce modèle. Dans le cadre de la théorie des extrêmes multivariés, nous proposons par une approche différente, d'assimiler les évènements 'anormaux' aux évènements extrêmes et de construire un modèle permettant la classification non supervisée de ces évènements. Pour ce faire, après avoir estimé le support de la structure de dépendance M, un modèle de mélange paramétrique pour la mesure angulaire Φ est proposé. L'intérêt de cette modélisation est double. Premièrement, cela permet d'attribuer de manière probabiliste des points à des sous-groupes de composantes. Deuxièmement, il devient possible de construire une matrice de similarité S = (s i,j ) i,j=1,...,n , où n est le nombre de points extrêmes, telle que: s i,j = P[ les points i et j sont attribués au même sous-groupe de composantes ] (1.38) L'intérêt de la matrice de similarité S est que l'on peut lui appliquer un algorithme de partitionnement spectral et ainsi faire du clustering sur les points extrêmes, non plus simplement en les attribuant à un sous-groupe de composante particulier.

Par souci de simplicité nous supposons ici min α∈M |α| > 1 ansi que M = M, le cas général étant développé dans la partie 5. Supposons que la fonction z → µ([0, z] c ) soit différentiable d fois sur [0, ∞) d \ {0}, alors Φ possède des densités sur l'intérieur du simplexe S d ( • étant la norme L 1 dans cette partie) ainsi que sur chacun des sous-simplexes S α (voir le Théorème 1 dans [START_REF] Coles | Modeling extreme multivariate events[END_REF]

φ α (w|m α , ν α ) = Γ(ν α ) i∈α Γ(ν α m α,i ) i∈α w ναm α,i -1 i .
(1.41)

       α∈M π α m α,j = 1 d , ∀j ∈ {1, . . . , d} j∈α m α,j = 1, ∀α ∈ M α∈M π α = 1 (1.43)
et la maximisation de la log-vraisemblance devient de ce fait non convexe. Cette difficulté est atténuée par le changement de variable suivant, pour tout α ∈ M:

ρ α,j = π j m α,j , ∀j ∈ α.
(1.44)

Ainsi, les trois contraintes précédentes sont réduites à:

α∈M ρ α,j = 1 d , ∀j ∈ {1, . . . , d}.
(1.45)

En effet, en posant π α := j∈α ρ α et m α,j := ρ α,j j∈α ρα , pour tout α ∈ M et pour tout j ∈ α, il est aisé de voir que (1.45) implique j∈α m α,j = 1, ∀α ∈ M et α∈M π α = 1, et donc (1.45) est équivalent à (1.42). 1 1.6. CONTRIBUTIONS La deuxième difficulté vient de la nature asymptotique du modèle à estimer, conséquence directe du problème explicité dans la section 1.3.1. Pour tout élément V , obtenu après transformation du jeu de données initial (1.6), le pseudo-angle V V appartient à l'intérieur du simplexe central S d . Pour palier ce problème, l'idée est de séparer V entre sa contribution en α, V α = (V j ) j∈α et ce que l'on pourrait nommer son résidu non-asymptotique, ε α c = (V j ) j / ∈α (figure (1.10)). L'estimation se fait via l'ajout d'une variable cachée Z à valeurs dans {0, 1} K telle que

π α k = P(Z k = 1) où α k ∈ M = {α 1 , . . . , α K } et K k=1 Z k = 1.
Un modèle sous-asymptotique est alors proposé, comprenant le modèle de mélange de Dirichlet auquel s'ajoute un modèle pour les résidus non-asymptotiques. Ainsi, pour α k ∈ M et pour v ∈ (1, ∞) d , la vraisemblance est donnée par:

p(v|z k = 1) = r -|α k |-1 k φ α k (w k |m α k , ν α k ) j / ∈α f ε (v j |λ k ), (1.46) où f ε (•|λ k ) est une distribution exponentielle de paramètre λ k , r k = j∈α k v j et w k = ( v j r k ) j∈α k . L'inférence

CONTRIBUTIONS

Articles de conférence avec actes

• Feature clustering for extreme events analysis, with application to extreme stream-flow data. (ECML 2016) Authors: M. Chiapino, A. Sabourin. Dans un premier temps, une méthode d'estimation robuste, non paramétrique, du support de la mesure limite est développée. Le procédé considère séquentiellement des régions emboitées de l'espace de sorte que le plus petit de ces ensembles satisfaisant un certain critère correspond à un ensemble maximal de composantes asymptotiquement dépendantes. Suivant les principes de l'algorithme 'Apriori', nous répondons à l'explosion combinatoire induite par les régions emboitées, par l'algorithme CLEF (1.4.2) qui prouve son efficacité à estimer le support de la structure de dépendance sur des données bruitées et en grande dimension. L'algorithme ouvre une voie générale pour la réduction de dimension dans un contexte d'extrêmes multivariés et permet l'adaptation de différents résumés du degré de dépendance asymptotique (κ α , η α , etc.) (1.4.3) afin de déterminer la dépendance des différents groupes de composantes α.

Dans un deuxième temps, une modélisation paramétrique de la structure de dépendance parcimonieuse a été proposée. L'écart (en loi) entre les données observées et la mesure limite est modélisé par un bruit (à queue légère), ce qui permet d'appliquer des méthodes de vraisemblances classiques fondées sur les modèles de mélange. De multiples façons de modéliser cet écart sont possibles et leur étude constitue un vaste champ de recherche. Le choix de la décomposition V observé = W R + bruit, où W est utilisé pour estimer la mesure angulaire, est motivé par sa relative simplicité d'adaptation à l'algorithme EM, et prouve son efficacité pour la classification non supervisée de données extrêmes.

Introduction

INTRODUCTION

Extreme value theory meets the need to build a model for extrapolation over a particular kind of rare phenomena. Given a set of observations, the goal is to foresee how likely future realisations will exceed the ranges previoulsy recorded. More precisely, our goal is to characterize the behavior of random variables above some large thresholds. These considerations are all the more justified when the distribution is heavy-tailed, so that the probability for a point to be extreme is significant. In the univariate case, extreme value theory can be summed up to the study of the maximum values of a process. The asymptotic behavior of the maximum of a random variable is well-known and is characterized within a parametric familly of distributions, the generalized extreme value distribution [START_REF] De Haan | On regular variation and its application to the weak convergence of sample extremes[END_REF]). In the multivariate case, unlike R, the space is not ordered and therefore it is not possible to define the maximum points of a given dataset. However, by means of thresholds exceedances, componentwise maximum or more generally by using the norm of the vectors, it is possible to provide a meaning to the extreme aspect of some points.

Several difficulties specific to the extreme value theory arise when the dimensionality increases. Usual methods of dimension reduction like principal component analysis (PCA) rely on the study of the covariance matrix. Yet, the covariance between features of heavy-tailed distribution is not defined and therefore such tools are of no use. In addition, on a statistic point of view extreme value theory has an inherent weakness: extreme points are rare by nature thus only a small proportion of the data should be used for the inference.

This introduction is organized as follow. In Section 2.2, the general framework of extreme value theory is exposed, in the univariate 2.2.1 and multivariate cases 2.2.2. In particular, we introduce the angular measure Φ which characterizes the dependence structure of any extreme value distribution after standardization of the marginals. A sparse representation of this measure is described in Section 2.3. On this basis, a general method to find the support of Φ is detailed in Section 2.4. In particular, the algorithm CLEF (CLustering Extreme Feature) is developed in 2.4.2, followed by alternative versions of the algorithm based on test statistics 2.4.3. Finally, a parametric model for the angular measure is proposed in Section 2.5.

2.2. EXTREME VALUE THEORY

EXTREME VALUE THEORY

Univariate extreme value theory

In the univariate case, extreme value theory characterizes the asymptotic behavior of the maximum of random variables valued in R. Let X ∈ R be a random variable distributed according to F and X 1 , . . . , X n n ≥ 1 i.i.d. (independent and identically distributed) copies of X. Suppose that our goal is to estimate large quantiles in the distribution tail, i.e. P[X > x p ] = p with p small. In the case where none of the observed variables lie in [x p , ∞), the only solution is to build a model based on extrapolation. The standard approach relies on the study of the maximum M n := max{X 1 , . . . , X n }. Given that the variable X is not bounded, the challenge is reflected by the following limit. For any x ∈ R:

P[M n ≤ x] = F n (x) -→ n→∞ 0, (2.1)
so the limiting distribution is necessarily degenerate. Usual methods based on the average behavior of the variable and naive empirical approaches are of no use. In order to tackle (2.1) one has to apply some normalization on the maximum M n to make a non-degenerate limit arise.

The following property is the basis of extreme value theory. It applies to the normalized maximum in a similar way than the central limit theorem (CLT) does with the normalized sum of i.i.d. random variables. The main difference is that the existence of the limit in the CLT is guaranteed as soon as the first two moments of the variable are finite. However, assuming the existence of two sequences (a n ) n≥1 and (b n ) n≥1 with a i > 0 such that:

lim n→∞ P M n -b n a n ≤ x = G(x), (2.2)
with G a non-degenerate distribution, then the limiting distribution G has necessarily the form:

G(x) = exp -1 + ξ x -µ σ -1 ξ + . (2.3)
F is said to be in the max-domain of attraction of the generalized extreme value (GEV) distribution G, written F ∈ DA(G). When ξ > 0, G is a heavy-tailed Fréchet, ξ = 0 corresponds to a light-tailed Gumbel and ξ < 0 to a bounded-tail Weibull. The asymptotic behavior of the maximum is therefore entirely characterized within a parametric familly. Therefore, inferential issues on such models are mostly limited to the proportion of upper data points that has to be considered.

Multivariate extreme value theory

In the multivariate case, identifying extreme points lacks of a proper definition. Indeed when d ≥ 2, R d is not an ordered set. Nevertheless, let x be a vector valued

CHAPTER 2. INTRODUCTION 2
in R d then one can consider x as extreme as soon as max j∈{1,...,d} x j > t or more generaly when x > t with t > 0 a large threshold and • a given norm. In a similar way than for the univariate case, multivariate extreme value theory is based on the limiting distribution of the standardized componentwise maximum of a random variable.

Let X 1 , . . . , X n be n ≥ 1 i.i.d. random variables valued in R d distributed according to F . Marginal distributions of F are named F j for j = 1, . . . , d and are assumed to be continuous. Let

M n := ( max i=1,...,n X i,1 , . . . , max i=1,...,n X i,d ) be the n-th component- wise maximum. We assume that there exists two sequences (a n ) n≥1 in (0, ∞) d and (b n ) n≥1 in R d , such that: lim n→∞ P M n,j -b n,j a n,j ≤ x j , j = 1, . . . , d = G 0 (x), (2.4)
where G 0 is a multivariate distribution with non-degenerate marginals G 0,j , j = 1, . . . , d. Then G 0 is characterized by the following formula [START_REF] Resnick | Extreme Values, Regular Variation, and Point Processes[END_REF]), for some

x 0 ∈ R d : G 0 (x) =      exp -µ 0 [x 0 , x] c , for x ≥ x 0 , 0 otherwise, (2.5)
where µ 0 is the so-called exponent measure, a Radon measure defined on

[x 0 , ∞] \ {x 0 } and with [x 0 , x] c = ∪ d j=1 y ∈ [x 0 , ∞] \ {x 0 } : y j > x j .
The multivariate distribution F is said to be in the domain of attraction of G 0 , i.e. F ∈ DA(G 0 ). Also, each marginal G 0,j is a univariate extreme value distribution, so that for all j ∈ {1, . . . , d}, F j ∈ DA(G 0,j ). However, the exponent measure µ 0 and therefore G 0 cannot be characterized within a unique familly of parametric distributions. This is a major difference with the univariate case.

Let X = (X 1 , . . . , X d ) be a random variable valued in R d according to F . The study of the tail of the distribution F splits into two separate parts. On the one hand, the characterization of the marginals which is related to the univariate case. On the other hand, the study of the dependence structure which is the main purpose of the multivariate case. In an analogous manner than the study of copulas, a natural preliminary step is to standardize each marginal to the same distribution. The choice of the distribution is somewhat arbitrary, though it is common to transform each marginal to a unit Pareto [START_REF] Resnick | Extreme Values, Regular Variation, and Point Processes[END_REF]). The standardization is given by the following transformation:

V j := (1 -F j (X j )) -1 , for j = 1, . . . , d. (2.6) Remark 2.1. Let X 1 , . . . , X n be n ≥ 1 i.i.d. copies of X.
In practice the marginals F j are usually not known. Thefore, we apply the transformation through the empirical distributions Fj (x) :=

1 n n i=1 1 {X i,j <x} where x ∈ [x 0,j , ∞). Let us call 2 2.2. EXTREME VALUE THEORY rank(X i ,j ) = n -n i=1 1 {X i,j <X i ,j } .
The transformation then becomes, for all j ∈ {1, . . . , d}:

Vi,j = (1 -Fj (X i,j )) -1 ) = n rank(X i,j )
.

From now on, we work with the transformed random vector

V = (V 1 , . . . , V d ) valued in (1, ∞) d
with unit Pareto marginals. The existence of the asymptotic distribution G 0 in (2.4) can be reformulated in a suitable way through regular variation. The random variable V is said to be regularly varying on (1, ∞) d with index -1 and with limiting measure µ defined on [0, ∞] d \ {0} if (e.g.Resnick ( 2013)):

tP t -1 V ∈ A ---→ t→∞ µ(A),
(2.7)

for any Borel set A in [0, ∞] d \{0} such that µ(∂A) = 0 and 0 / ∈ ∂A. The distribution of V is in the domain of attraction of the multivariate extreme value distribution G, such that for all v ∈ [0, ∞] d \ {0}: G(v) = exp[-µ([0, v] c )].
(2.8)

The link between the asymptotic distribution G 0 and G can be written:

G(v) = G 0 (G ← 0,1 (e -1/v 1 ), . . . , G ← 0,d (e -1/v d )),
where G ← 0,i is the inverse function of G 0,i . Remark 2.2. Let t > 0 be a large threshold, given that for some j 0 ∈ {1, . . . , d},

V j 0 is extreme, i.e. V j 0 > t, the function x → µ [0, x] c is approximatively, for v ∈ [0, ∞) d \ {0}
, the following conditionnal probability:

P V 1 > tv 1 or . . . or V d > tv d |V j 0 > t = P V ∈ t[0, v] c |V j 0 > t = tP V ∈ t[0, v] c ≈ µ [0, v] c = -log G(v)
The benefit of the transformation (2.6) along with the regular variation assumption (2.7) is to focus on the dependence structure that is entirely characterized by the exponent measure µ. We recall that the class of such measures is not embedded within a parametric familly. Nevertheless, as a consequence of (2.6) and the maxstable nature of any extreme value distribution [START_REF] Resnick | Extreme Values, Regular Variation, and Point Processes[END_REF]), the measure µ is homogeneous of degree -1:

µ(tA) = t -1 µ(A),
(2.9) for any t > 0 and for all Borel set A bounded away from the origin, i.e. 0 / ∈ ∂A. One can deduce from (2.9) that the marginals of G are unit Fréchet. Indeed, for any j ∈ {1, . . . , d} and x j ∈ (0, ∞):

G(∞, . . . , x j , . . . , ∞) = exp -µ ([0, ∞] × . . . × [0, x j ] × . . . [0, ∞]) c = exp - 1 x j µ A , with A = ([0, ∞] × . . . × [0, 1] × . . . [0, ∞]) c
. And we get µ(A) = 1 using the limit (2.7) along with the fact that V has unit Pareto marginals:

µ(A) = lim t→∞ tP t -1 V ∈ A = lim t→∞ tP (V j > t) = 1.
Remark 2.3. It follows from the homogeneity property (2.9) that the inclusion-

wise mininal sets of [0, ∞] d \ {0} with non-zero µ-mass are necessarily cones, i.e. C ⊂ [0, ∞] d \ {0} such that x ∈ C ⇒ tx ∈ C for all t > 0.
The homogeneity property is fundamental in order to characterize the class of exponent measures resulting from the standardization of the marginals. Indeed (2.9) leads to a pseudo-polar decomposition of µ. For all v ∈ [0, ∞) d \ {0}, let us consider the bijective function:

T : v → ( v , v v ),
where • is a given norm on R d . It is then possible to define a measure Φ on the positive orthant of the unit sphere S d = {x ≥ 0 : x = 1}, the so-called angular measure, which characterizes the dependence structure of V :

Φ(A) := µ({v : v > 1, v v ∈ A}), (2.10)
for any set A ⊂ S d . Considering the homogeneity property of µ (2.9) one can write:

µ({v : T (v) ∈ (t, ∞) × A}) = µ({v : v > t, v v ∈ A}) = t -1 Φ(A).
Finally, this last equation along with the bijectivity of T imply a factorization of µ in two independent measure [START_REF] De Haan | Limit theory for multivariate sample extremes[END_REF]):

µ • T -1 (dr, dw) = r -2 drΦ(dw).
(2.11)

Remark 2.4. For large values of V , the angular measure Φ quantifies the relative contribution of each coordinate V j in the extremal behavior of V . One can say that Φ measures the directions of (1, ∞) d in which V is likely to be extreme. Let us consider the pseudo-polar decomposition R = V and W = R -1 V . One can re-write the limit (2.7) depending on Φ:

P (W ∈ A, R > rt | R > t) ---→ t→∞ r -1 Φ(S d ) -1 Φ(A), (2.12)
for all Borel set A in S d such that Φ(∂A) = 0 and r > 1. In other words, when the radial part R is large, then R and the pseudo-angle W are approximately independent. The distribution of W is approximately (and up to a normalization) the angular measure and R is approximately distributed according to a unit Pareto.

As for the exponent measure, there is no parametric familly that can entirely describe the class of angular measures. The only condition on Φ in order to be a proper angular measure, is the so-called moment constraint:

S d w j Φ(dw) = 1, ∀j ∈ {1, . . . , d}.
(2.13) Indeed, using the decomposition (2.11), one can rewrite

x → µ([0, x] c ) depending on Φ. Let x be in [0, ∞] d \ {0}: µ([0, x] c ) = 1 {∃j:u j >x j } µ(du) = S d ∞ 0 1 {r>min j x j w j } r -2 drΦ(w) = S d (min j x j w j ) -1 Φ(w) = S d max j w j x j Φ(w).
Therefore, considering x = (∞, . . . , x j , . . . , ∞), with x j ∈ (0, ∞) we have:

1 x j = -log G(x) = µ([0, x] c ) = S d w j x j Φ(dw).
The constraint (2.13) is thus a necessary condition on Φ to be a proper angular measure. The proposition 5.11 in [START_REF] Resnick | Extreme Values, Regular Variation, and Point Processes[END_REF] demonstrates that it is also a sufficient condition. As a straight consequence of the moment constraint (2.13), in the case where the L 1 norm is chosen for

• , i.e. x = |x 1 | + . . . + |x d |, we have Φ(S d ) = d. Indeed: Φ(S d ) = S d Φ(dw) = S d d j=1 w j Φ(dw) = d j=1 S d w j Φ(dw) = d.

Previous works and issues

On a statistic point of view, various models have been proposed for the extreme value distribution G through the modelization of the function [START_REF] Coles | Modeling extreme multivariate events[END_REF]). Two main strategies arise to infer such models. Methods based on the componentwise maximum are used to estimate G while threshold exceedence methods are used to estimate µ. Componentwise maximum methods face a major issue to estimate the likelihood when the dimensionality increase (d > 10). Indeed, one has to compute:

x → µ[0, x] c , i.e. G(x) = exp[-µ[0, x] c ] (see e.g.
∂ d ∂x 1 . . . ∂x d e -µ([0,x] c ) ,
which is a sum whose number of terms explodes with the dimension (see e.g. [START_REF] Huser | Likelihood estimators for multivariate extremes[END_REF]). Some simplifications have been proposed (e.g. [START_REF] Wadsworth | On the occurrence times of componentwise maxima and bias in likelihood inference for multivariate max-stable distributions[END_REF], [START_REF] Stephenson | Exploiting occurrence times in likelihood inference for componentwise maxima[END_REF]), but the computation remains infeasible when d > 10.

Either for componentwise maximum methods or for threshold exceedence, most of these approaches only consider the case of total dependence between the variables. In other words, the case where the support of µ is infinite on any hyperplane of [0, ∞] d \{0}. Also, in high dimension, non-parametric estimations of these measures are not spared by the general issue of the curse of dimentionality which is all the more burdened by the small proportion of points (i.e. the extreme points) that should be used for inference. The modelization of the angular measure, which is in theory equivalent to the modelization of the exponent measure, lends itself well to methods aiming at reducing the dimension (see Section 2.3).

Estimation of the angular measure Several parametric and non-parametric approaches have been proposed for the estimation of the angular measure. The class of Dirichlet mixture distributions on the simplex has been proposed in (Boldi and Davison (2007b)) and [START_REF] Sabourin | Bayesian model averaging for multivariate extremes[END_REF]). This class of models is dense (in the weak sense) in the class of angular measures. Non-parametric models have been proposed in [START_REF] Guillotte | Non-parametric bayesian inference on bivariate extremes[END_REF]), Fougeres et al. (2013)). However, experiments are only made in moderate dimensions (d ≈ 5). Recently proposed works intend to 2

SPARSE ANGULAR MEASURE

go through this issue with adapted clustering methods [START_REF] Chautru | Dimension reduction in multivariate extreme value analysis[END_REF]). Finally, an algorithm developped in (Goix et al. (2015a), Goix et al. (2016a)) is used to bring out a sparse support of the exponent measure (and equivalently of the angular measure) in greater dimension (d ≈ 50).

SPARSE ANGULAR MEASURE

In high dimension (d ≈ 100), it is reasonable to assume that an extreme phenomenon described by a random variable V = (V 1 , . . . , V d ) is not due to all its features simultaneously, i.e. V j > t for all j ∈ {1, . . . , d} where t > 0 is a large threshold. More precisely, we make the assumption that the extreme nature of the phenomenon is only due to some particular subgroups of coordinates. Let us assume that V is large, then there exists a subgroup α ⊂ {1, . . . , d} such that V ≈ V α , where |α| d and V α = (V j ) j∈α . The aim is to bring to the fore the multimodal nature of the extremal behavior of V and thus the plurality of subgroups α that identify the asymptotically dependent coordinates. It is also reasonable to assume that the number of such subgroups is small compared to the total number of subsets of {1, . . . , d}, that is 2 d -1. More formally, let Φ be the angular measure associated with V , then the dependence structure of V is characterized through the distribution of the mass of Φ over the positive orthant of the unit sphere. In the case of a sparse dependence structure, the mass of Φ is only distributed over some particular subspaces of the boundary of S d , which correspond to asymptotic dependent subgroups of features. Indeed, let us assume that at least two coordinates V i and V j are asymptotically independent, which write:

tP (V i > t, V j > t) = tP t -1 V ∈ {v : v i > 1, v j > 1} ---→ t→∞ 0. Now for all w ∈ Sd , there exists v ∈ [0, ∞) d \ {0} such that v v = w, v i > 1 and v j > 1, therefore: Φ( Sd ) = µ {v : v > 1, v v ∈ Sd } ≤ µ {v : v i > 1, v j > 1} = 0.
This toy example shows that asymptotical independence between variables implies that the mass of the angular measure is confined on the boundary of S d . Let us call S α ⊂ S d the subspace associated with the coordinates α ⊂ {1, . . . , d}:

S α := {w ∈ S d : ∀j ∈ α w j > 0, ∀j / ∈ α w j = 0}, (2.14)
and M the set of all subgroups of coordinates of {1, . . . , d} that correspond to subspaces with non-zero Φ-mass:

M := {α ⊂ {1, . . . , d} : Φ(S α ) > 0}.
(2.15)

The sparsity assumption can then be stated through the double inequality:

   max α∈M |α| d, |M| 2 d -1.
(2.16)

Issue on the estimation of M.

Let V be the random variable valued in (1, ∞) d with unit Pareto marginals and associated with the angular measure Φ. We apply the pseudo-polar decomposition (R, W ) where R = V and W = V V . It is essential to point out that the limit (2.12), where Φ is the limiting measure of W for R > t with t → ∞, is only satisfied on continuity sets of Φ, i.e. any Borel set A of S d such that Φ(∂A) = 0. In particular, the limit is only established on subspace of S d that has a non-zero Lebesgue measure. Therefore it does not apply on subspaces S α such that Φ(S α ) > 0. Indeed, in S d we have S α ⊂ ∂S α , and thus for α ∈ M:

0 < Φ(S α ) ≤ Φ(∂S α ).
(2.17)

Furthermore, the sparcity assumption (2.16) implies that the angular measure assigns no mass on the interior of S d , i.e. Φ( Sd ) = 0. On the other side, W is always valued in Sd i.e. P(W ∈ Sd ) = 1. In this situation, a naive empirical estimation of Φ using an i.i.d. sample (W i ) i=1,...,n is no longer possible as it was suggested by the remark 2.4.

Sparse angular measure for the logistic model

Let us consider the logistic model [START_REF] Coles | Modeling extreme multivariate events[END_REF]) which constitutes, along with the asymmetric logistic model, a very flexible class of parametric extreme value distributions for general dimensions. For all v in (0, ∞) d :

µ lgtc,θ ([0, v] c ) = d j=1 v -1 θ j θ , (2.18)
where θ ∈ (0, 1]. Asymptotic independence is reached for θ = 1 and in this case the mass of the angular measure is only located on the axes, i.e. Φ(e j ) = 1 for j ∈ {1, . . . , d} where e j ∈ R d such that e j i = 1 for i = j and e j i = 0 otherwise. Nonetheless, it is obvious that any point distributed according to the extreme value

distribution G(v) = exp -µ lgtc,1 ([0, v] c ) is almost surely valued in (0, ∞) d and,
as well, the pseudo-angle almost surely lies in the interior of S d .

Remark 2.5. The extreme value distribution G and the exponent measure µ are only linked through the equation (2.8), so the support of µ (as well as the support of the angular measure) is not the support of G. Therefore, the fact that the support

ESTIMATION OF THE SUPPORT OF THE ANGULAR MEASURE

of Φ is included in sub-dimensional spaces of S d does not mean that the support of G is included in hyperplans of [0, ∞] d \ {0}. In other words, even if V was directly distributed according to G, the angular measure only appears as a limit of the conditional probability given that V goes to infinity.

Modeling Φ instead of G has the advantage to exploit a natural representation of sparse dependence structures. On the other hand, the support we aim to estimate is never reached in the space of the observations (Remark 2.5). Two main ways are considered in the following to tackle this issue. In Chapter 3 and Chapter 4, we approach the sub-parts of interest of the boundary with subspaces of non-zero Lebesgue measure, and we apply a non-parametric estimation of the support of the angular measure (or equivalently of the exponent measure) with a counting measure on these subspaces. In Chapter 5, we propose a parametric mixture model for the angular measure in which each mode corresponds to a subgroup of features asymptotically dependent. In order to infer such a model, we project the observations on the corresponding subspaces and we consider the residual as an orthogonal noise.

ESTIMATION OF THE SUPPORT OF THE ANGULAR MEASURE

In practice, when the only matter is to estimate the support of the angular measure on S d , it is more suitable to work directly in the definition space of the exponent measure µ, i.e. [0, ∞] d \ {0}. By definition of Φ in equation (2.10), the counterpart of the subspace S α where α ⊂ {1, . . . , d}, expressed in terms of subspaces of [0, ∞] d \ {0}, is the truncated cone:

C α = {x ≥ 0 : x > 1, ∀j ∈ α x j > 0, ∀j / ∈ α x j = 0}. (2.19)
Indeed, for all α in {1, . . . , d}:

Φ(S α ) = µ(v : v > 1, v v ∈ S α ) = µ(C α ).
So that finding the support of Φ amounts to find the support of µ. Note that the set of all the cones C α forms a partition of {v ∈ R d

+ : v > 1}. The figure (2.1) shows a partition of R 3 + \ [0, 1] 3 by the cones C α . Remark 2.6. By considering the partition ∪ α⊂{1,...,d}
C α we browse over all possible associations of asymptotically dependent coordinates of {1, . . . , d}. Indeed, let α ⊂ {1, . . . , d} be a subgroup such that a non-null probability is attributed to the event: { All the variables (V j ) j∈α are large while all the complementary variables (V j ) j / ∈α are small. } Then there exists a cone C ⊂ C α such that µ(C) > 0. By extension, let C tot ⊂ [0, ∞] d \ {0} be the cone containing the whole mass of the exponent measure, then C tot is necessarily a union of sub-cones themselves included in some cones C α .

Inference

The inferential issue over the mass of Φ, explained in Section 2.3.1, transposes to the estimation of the mass of µ on the subspaces of [0, ∞] d \ {0}. We consider again the random variable V regularly varying in (1, ∞) d with limiting exponent measure µ and unit Pareto marginals. As V is almost surely valued in (1, ∞) d , for all t > 0 and a strict subset α ⊂ {1, . . . , d}:

P V t ∈ C α = 0. (2.20)
Moreover, the cones C α are of empty interior and therefore µ(∂C α ) > 0. In order to estimate µ(C α ) using i.i.d. copies of V , it is necessary to thicken the truncated cones C α in the interior of [0, ∞] d \ {0}. When V is large, the idea is to locate the subspaces of the positive orthant (i.e. C α ) toward which V settle down. We therefore consider the following thickened cones, for > 0 and α ⊂ {1, . . . , d}:

R α = {x ≥ 0 : x > 1, ∀j ∈ α x j > , ∀j / ∈ α x j ≤ }.
(2.21)

Remark 2.7. Note that for ∈ (0, 1) the set of all cones R α remains a partition of

[0, 1] c (for the L ∞ -norm). Assuming that µ(C α ) > 0, the inclusion-wise minimal set C ⊂ C α such that µ(C) > 0, is a cone (Remark 2.3). Yet, any cone included in C α is necessarily of non-empty intersection with R α . Therefore, the thickened cone R α satisfies µ(C α ) > 0 ⇒ µ(R α ) ≥ µ(C) > 0 and is furthermore of non-zero Lebesgue measure.
More formally, as a consequence of the continuity from above of any measure:

µ(C α ) = µ(∩ >0, ∈Q R α ) = lim →0 µ(R α ).
However, in order to apply the limit (2.7) to the sets R α , they must verify µ(∂R α ) = 0. The boundary of R α in [0, ∞] d \ {0} is composed of hyperplans parallel to the axes. For different values of ∈ (0, 1) the boundaries are disjoint. Therefore, there is at most a countable number of such that µ(∂R α ) > 0, otherwise µ would be infinite on a compact set of [0, ∞] d \ {0}. It is thus possible to choose a threshold > 0 arbitrarily close to zero such that µ(∂R α ) = 0, and estimate the quantity µ(C α ) through the limit:

µ(C α ) = lim →0 lim t→∞ tP (V ∈ tR α ) .
An algorithm named DAMEX (Goix et al. (2016a)) have been developed to estimate the mass µ(C α ) by a counting method on the subspaces tR α for a small > 0 and a large threshold t > 0. More precisely, let V 1 , . . . , V n be n ≥ 1 i.i.d. copies of the random variable V . The estimator of µ(C α ) is defined as:

μα = 1 k n i=1 1 {V i ∈ n k R α } , (2.22)
where t is replaced by n k , with k ≤ n an order of magnitude of the number of extreme points considered in the sample. A low threshold µ min > 0 is then chosen in order to decide:

µ(C α ) > 0 if μα ≥ µ min .
The algorithm is tested on real data (a wave directions dataset from the North Sea provided by Shell) and achieves to reduce drastically the dimensionality of the problem by gathering most of the mass of µ in a limited number of cones C α of low dimensions. Similarly, the method shows good results on simulated data. A dataset is generated from a sparse version of the asymmetric logistic model (Tawn (1990a)):

G(z) = exp - α∈M j∈α (|A(j)|z j ) -1/θα θα , (2.23)
where A(j) = {α ∈ M : j ∈ α} and with θ α ∈ (0, 1] a dependence parameter such that θ α = 1 corresponds to asymptotic independence and θ α ↓ 0 to total dependence. The algorithm easily recovers the support M of the dependence structure for strong asymptotic dependence, i.e. θ α = 0.1. However, when the dimension d, θ α and the number |M| of dependent subgroups increase, the method is no longer able to recover M. Due to the partition of the space in a huge number of subcones (2 d -1), the model is inherently sensible to any variability of the subgroups of coordinates involved in the extremal behavior of V . Indeed, let α 0 be a subgroup of coordinates such that µ(C α 0 ) > 0. Then it is possible that for several coordinates j ∈ {1, . . . , d} \ α 0 the quantity μα 0 ∪{j} is significant:

μα 0 ∪{j} = 1 k n i=1 1 {V i ∈ n k R α 0 ∪{j} } ≥ µ min . (2.24)
As a result, the mass is scattered into clusters of sub-cones that are close to each other and it is no longer possible to recover the sparse dependence structure. It is therefore necessary to build a more robust method to estimate the support of µ in the cases of moderate and low dependence between the components (i.e. θ α ∈ [0.5, 1)) or more generaly in order to tackle noisy datasets. In order to exhibit this scattering of the mass, we generate several dataset (V 1 , . . . , V n ) with n = 1e5, from an asymmetric logistic distribution for different levels of dependency θ α . We set d = 50 and we randomly generate 50 subgroups of coordinates to form M. The following figures (2.2) correspond to the different datasets and show the repartition of extreme points in the cones n k R α where k is chosen such that the proportion of extreme points

n i=1 1 { V i > n k }
represents 5% of the dataset. The x-axis represents the subgroups of coordinates α such that μα > 0. The y-axis represents the values of the estimators μα . The elements of M are represented in red. In order to recover M the threshold µ min > 0 must verify:

   μα ≥ µ min , for all α ∈ M μα < µ min , otherwise.
(2.25)

For θ α = 0.1 (2.2 up-left) there is a big gap between the 'true' subgroups (i.e. in M) and the others. Thus, it is possible to chose µ min > 0 such that (2.25) is verified. When θ α increases, corresponding to lower asymptotic dependences, we observe a scattering of the extreme points. For θ α = 0.7, 0.9 (2.2 bottom-left (resp. right)) it is no longer possible to recover M.

CLustering Extreme Feature Algorithm

The inner weakness of the partition of the space into cones is that it leads to differenciate close extreme events. Instead, we decide to only look for subgroups α ⊂ {1, . . . , d} such that V j is large for all j ∈ α, regardless of the counterparts V j for j ∈ {1, . . . , d} \ α. We define the nested rectangles:

Γ α = {x : ∀j ∈ α, x j > 1}.
(2.26)

One can show that µ(∂Γ α ) = 0, just as for the cones R α , and thus we can use the limit (4.1) to estimate their mass. The advantage of this approach is to avoid the scattering of the empirical mass of µ as in (2.24). For all α ⊂ {1, . . . , d}, Γ α corresponds to the event: { All the variables (V j ) j∈α are large (regardless of their counterpart (V j ) j / ∈α ). } Let t > 0 be a high threshold, each rectangle tΓ α gathers extreme points with large coordinates in α regardless of the complementary coordinates {1, . . . , d} \ α. One simple but important property is that any point in tΓ α also belongs to tΓ β for all β ⊂ α. The goal is to find, in an incremental way, the larger subgroups α such that the number of points in tΓ α is still significant. When we consider the inclusion-wise maximal subgroups α ⊂ {1, . . . , d} such that µ(C α ) > 0, the new framework becomes equivalent to the partition of the space into cones. On the one hand, let us assume that µ has a non-zero mass on C α . The inclusion-wise minimal region C ⊂ C α such that µ(C) > 0 is a cone (2.3). As any cone included in C α has a non-empty intersection with Γ α , µ(C α ) > 0 implies µ(Γ α ) > 0. On the other hand, let us assume µ(Γ α ) > 0. If for all β such that α ⊂ β, µ(Γ β ) = 0 then necessarily µ(C α ) > 0. By contraposition, if we assume µ(Γ α ) > 0 and µ(C α ) = 0, then there exists j C ∈ {1, . . . , d}\α and a cone C included in Γ α∪{j C } such that µ(C) > 0. In other words, there exists β such that α ⊂ β and µ(Γ β ) > 0.

Hence, maximal subgroups α ⊂ {1, . . . , d} such that µ(C α ) > 0, are also maximal subgroups such that µ(Γ α ) > 0. Let M be the set of maximal subgroups α ⊂ {1, . . . , d} such that µ(Γ α ) > 0, then we have the following property:

α ∈ M ⇔ α is maximal in M.
(2.27)

The counterpart of this approach is the implied combinatorial explosion due to the nested rectangles. Indeed, the maximal number of cones C α to estimate is at most the number of extreme points. On the contrary, the potential number of rectangles Γ α that might be browsed is much bigger, i.e. O(2 d 0 -1) for d 0 ≤ d. This issue arises from the fact that ∪ α⊂{1,...,d} Γ α is no longer a partition of the space.

'Apriori' Algorithm

More generaly this exploratory issue can be reformulate in the following manner: we want to find all subgroups I of a set of indexes T that verify a condition C(I). The essential property is that for all I 1 ⊂ I 2 then C(I 2 ) implies C(I 1 ). This problem is directly linked to the Apriori algorithm introduced in [START_REF] Agrawal | Fast algorithms for mining association rules[END_REF]).

The algorithm reduces drastically the browse of all subset of T with an incremental procedure. Applied to our framework, the procedure browses subsets α ⊂ {1, . . . , d}, |α| ≥ 2, of increasing sizes and only keeps those that verify a criterion C(α) equivalent to µ(Γ α ) > 0. The output of the algorithm is thus:

M 0 := α ⊂ {1, . . . , d} : |α| ≥ 2, C(α) .
(2.28)

In order to get M, a postprocessing step is done to only keep the maximal elements of M 0 for inclusion. This process can be visualized through the Hasse diagram (2.3).

We explore the diagram top to bottom and prune all the edges of a node as soon as the node does not verify the criterion C. This method reduces notably the potential number of subgroups to be tested. For the criterion µ(Γ α ) > 0, the pruning of the graph follows the property:

Let α be such that µ(Γ α ) = 0, then for all β ⊃ α, µ(Γ β ) = 0.

(2.29)

A s = β ⊂ {1, . . . , d} : |β| = s and β \ {j} ∈ M s-1 0 for all j ∈ β .
(2.30)

For s = 2, we put A 2 = α ⊂ {1, . . . , d} : |α| = 2 .

Stopping criteria

Let (V 1 , . . . , V n ) be an i.i.d. sample of random variables valued in (1, ∞) d with limiting angular measure Φ. In order to decide whether or not α ⊂ {1, . . . , d} belongs to M 0 , one natural criterion would be:

γ α = 1 k n i=1 1 {V i ∈ n k Γα} ≥ γ min , (2.31)
with γ min > 0 a small threshold and k ≤ n an order of magnitude of the number of extreme points considered. However, as the rectangles are nested, the value of

tP[V ∈ tΓ α ] ≈ µ(Γ α ), with t > 0, is necessarily decreasing for increasing sizes of α. Let us consider α 1 ⊂ α 2 ⊂ {1, . . . , d}, so that Γ α 2 ⊂ Γ α 1 and thus tP[V ∈ tΓ α 2 ] ≤ tP[V ∈ tΓ α 1 ] for all t > 0.
The choice of the threshold should therefore depends on the size of the subsets considered. To tackle this issue, an alternative summary of the dependency degree between the coordinates α ⊂ {1, . . . , d} is proposed:

κ α := µ(Γ α ) µ( β⊂α |β|=|α|-1 Γ β )
.

(2.32)

The quantity κ α does not necessarily decrease with the size of α, and we have κ α > 0 if and only if µ(Γ α ) > 0. It is then possible to chose a threshold κ min > 0, that does not depend on the size of α, in order to get a criterion {κ α ≥ κ min } to decide {κ α > 0}. The idea is to make the acceptation of a group α in M 0 dependent on the subgroups of α:

κ α = P[ All the variables V j for j ∈ α are large | All variables are large except at most one ]
The estimator κα of the quantity (2.32) is then:

κα = γα μ( β⊂α |β|=|α|-1 Γ β )
.

(2.33)

The algorithm that estimates M is described in (2).

Algorithm 2 CLEF (CLustering Extreme Features)

Input: Threshold κ min > 0.

STAGE 1: Construct the set M 0 of tail-dependent groups.

Step 1:

M 2 0 = α ⊂ {1, . . . , d} : |α| = 2, κα > κ min , S = 2. Step s = 3, . . . , d: If M s-1 0 = ∅, end PHASE 1. Otherwise:
• Generate candidate of size s:

A s = {α ⊂ {1, . . . , d} : |α| = s and α \ j ∈ M s-1 0 for all j ∈ α}.

• Keep all subgroups that verify the criterion, M s 0 = α ∈ A s : κα > κ min .

• Si M s 0 = ∅, S = s.

Sortie: M 0 = ∅ if S = 1 and M 0 = S s=2 M s 0 if S ≥ 2. STAGE 2: Only keep maximal α. If S = 1, then M = ∅. Otherwise: Initialization: M ← M s 0 . for s = (S -1) : 2, for α ∈ M s 0 , If there is no β ∈ M such that α ⊂ β, then M ← M ∪ {α}. Output: M 2.4.2

.3 Results

In practice the choice of a criterion based on κ α is all the more justified when noise is observed on the components. In order to reproduce a dataset that can be assimilated to a real-world noisy dataset, we consider again the asymmetric logistic model described in 2.4.1, for d = 20, |M| = 15 and θ α = 0.5, on which we artificially disturbe the dependence structure. More precisely, each points V i ∈ (0, ∞) d , i = 1, . . . , n, is generated with a slightly modified version M i of the subgroups of dependent coordinates M. For i = 1, . . . , n, each subgroup of coordinates α i k in M i is randomly altered from the original α k , so that there exists j 0 ∈ {1, . . . , d} such that:

{j 0 ∈ α k and j 0 / ∈ α i k } or {j 0 / ∈ α k and j 0 ∈ α i k }.
(2.34)

The goal is to recover M from the disturbed dataset. At each step s = 2, . . . , S, with S = max α∈M |α|, it is therefore necessary to distinguish the following subgroups, using a threshold κ min > 0:

κα ≥ κ min for all α ∈ M s 0 κα < κ min for all α ∈ A s \ M s 0 . We call M F = ∪ s=1,...,S
A s \ M 0 the set of all candidates that will not verify the criterion. Note also that M 0 = ∪ S s=2 M s 0 . The condition for the existence of a threshold κ min > 0 such that M = M is therefore:

max α∈M F κα < min α∈M 0 κα . (2.35) If (2.35) is verified, it is possible to choose κ min in the interval ( max α∈M F κα , min α∈M 0 κα ) to recover M. Figures (2.4) (resp.
(2.5)) displays the value of the empirical mass γα (resp. κα ) on the elements of M 0 (in red) and M F (in blue). The blue line represents max α∈M F γα (resp. max α∈M F κα ) while the red line represents min α∈M 0 γα (resp. min α∈M 0 κα ). The red subgroups are ranked by increasing sizes and it is clear that the empirical mass (2.31) decreases with the size of the subsets. Moreover, as max α∈M F γα > min α∈M 0 γα , it is not possible to find a threshold γ min > 0 whereby we can recover M.

We then compare the algorithm DAMEX along with the algorithm CLEF for the different criterions C(α) = {γ α > γ min } and C(α) = {κ α > κ min } on a dataset of daily recorded river flows. The dataset comes from d = 92 stations spread accross the different rivers of France, recorded from january 1 st 1969 to december 31 st 2008.

It is composed of n = 14610 vectors (X 1 , . . . , X n ) valued in R d + .
The goal is to find the subgroups of stations that are likely to be concomitantly flooded. After the transformation (2.6), we consider the extreme points ( V1 , . . . , Vnextr ) where

n extr = n i=1 1 Vi > n k with k chosen such that n extr = n • 5%. Figures (2.

Coefficient of tail dependence

The purpose of this section is to build a serie of test statistics based on the hypothesis µ(Γ α ) > 0. Using the criterion κ α > 0 raises an issue. Indeed, as soon as µ(Γ α ) = 0 the limiting distribution of the statistics √ k( κ α -κ α ) is degenerate. Thus, we have no control of the asymptotic levels of the tests. Let us consider the tail dependence coefficient η α ∈ (0, 1] introduced in the bivariate case in [START_REF] Ledford | Statistics for near independence in multivariate extreme values[END_REF]) and extended to general dimensions d ≥ 3 in (De Haan and Zhou (2011)) et [START_REF] Eastoe | Modelling the distribution of the cluster maxima of exceedances of subasymptotic thresholds[END_REF]). The fundamental assumption stipulates the existence of η α ∈ (0, 1] and a slowly varying function L α such that:

P[V ∈ tΓ α ] = t -1/ηα L α (t).
(2.36)

Under the assumption that both the limit lim t→∞ tP[V ∈ tΓ α ] = µ(Γ α ) exists and (2.36) is verified, µ(Γ α ) > 0 implies η α = 1. On the contrary, suppose (2.36) and lim inf t→∞ L α (t) > 0 then η α = 1 implies µ(Γ α ) > 0. In other words, the null hypothesis µ(Γ α ) > 0 corresponds to the hypothesis η α = 1 under mild conditions on L α . Hence, if η α = 1 the limit √ k( η α -η α ) is non-degenerate and it is possible to control the asymptotic levels of the associated test.

A new criterion C(α) based on the estimator of the coefficient of tail dependence η α is used in the algorithm CLEF. This series of test statistics is developed in Chapter 4 for different non-parametric estimators of η α , namely a multivariate extention of the Peng estimator [START_REF] Peng | Estimation of the coefficient of tail dependence in bivariate extremes[END_REF]) along with the Hill estimator [START_REF] Draisma | Tail dependence in independence[END_REF], [START_REF] Draisma | Bivariate tail estimation: dependence in asymptotic independence[END_REF]). A non-degenerate version of the test H 0 : κ α = 0 is also developed by adding a threshold κ min > 0, to get the new hypothesis H 0 : κ α > κ min .

ηα > 1 -q 1-δ σα √ k ,
where q 1-δ is the (1 -δ)-quantile of the normal distribution, σα is the variance of the estimator and k is an order of magnitude of the number of extreme points. Figure 2.9 displays the values of the quantity ηα -(1 -q 1-δ σα √ k ) over elements of the sets M F and M 0 for a low value δ = 1e -7. We see that the quantity of interest is strictly positive for any elements of M 0 (in red) and strictly negative on the elements of M F (in blue).

PARAMETRIC MODELING OF THE ANGULAR MEASURE

The motivation of this part is to detect and classify anomalies of a flight dataset provided by Airbus. The general approach of anomaly detection is to build a model that describes the 'normal' behavior of a phenomenon. On this basis, some events are defined as 'anomalies' based on the deviation they have with the estimated model of the 'normal' behavior. We propose a different approach, using the extreme value theory framework, by assimilating the 'anomalies' with the extreme events. This method allows to apply clustering on the 'anomalies'. To that end, once we have estimated the support of the dependence structure M, a parametric mixture model is proposed for the angular measure Φ. This modelization has two benefits. First, it allows to assign the points to subgroups of features in a probabilistic manner. Secondly, it is possible to build a similarity matrix S = (s i,j ) i,j=1,...,n , where n is the number of extreme points, such that: s i,j = P[ Points i and j are assign to the same subgroup of features. ]

(2.37)

Methods of spectral clustering are then applied on the similarity matrix S in order to cluster the extreme points.

For the sake of simplicity we assume here that min α∈M |α| > 1 along with M = M, the general case being developed in Section 5. Suppose that

z → µ([0, z] c ) is d-times differentiable on [0, ∞) d \ {0}
, then Φ has densities on the interior of the simplex S d ( • being the L 1 in this part) and on each sub-simplex S α (see Theorem 1 in [START_REF] Coles | Modeling extreme multivariate events[END_REF]). We can thus apply the decomposition:

1 Φ(S d ) Φ(dw) = α∈M π α φ α (w α )dw α (2.38)
where π α ∈ (0, 1) is a weight, such that α π α = 1, and φ α is a non zero density on S α . In order to correspond to a proper angular measure the mixture model has to verify the moment constraint (5.8). The densities φ α along with the weights π α have to verify:

2.5. PARAMETRIC MODELING OF THE ANGULAR MEASURE

α∈M j∈α π α Sα w j φ α (w α )dw α = 1 d , ∀j ∈ {1, . . . , d}. (2.39)
As the Dirichlet distribution is a natural distribution on the simplex we propose a Dirichlet mixture for the parametric modelling of Φ. The distribution can be parametrized by a mean parameter m α ∈ S α and a concentration parameter ν α > 0. We have, for all w ∈ S α :

φ α (w|m α , ν α ) = Γ(ν α ) i∈α Γ(ν α m α,i ) i∈α w ναm α,i -1 i . (2.40)
Therefore, given that Sα w φ α (w|m α , ν α )dw = m α , the moment constraint becomes: Two main issues arise with the inference of the model. First, the estimators m and π have to verify the constraints:

α∈M π α m α,j = 1 d , ∀j ∈ {1, . . . ,
       α∈M π α m α,j = 1 d , ∀j ∈ {1, . . . , d} j∈α m α,j = 1, ∀α ∈ M α∈M π α = 1 (2.42)
and therefore the log-likelihood maximization is non-convex. This difficulty is eased by the following substitution, for all α ∈ M:

ρ α,j = π j m α,j , ∀j ∈ α.
(2.43) Thus, the three constraints become:

α∈M ρ α,j = 1 d , ∀j ∈ {1, . . . , d}.
(2.44) Indeed, considering π α := j∈α ρ α and m α,j := ρ α,j j∈α ρα , for all α ∈ M and for all j ∈ α, it is straightforward that (2.44) implies j∈α m α,j = 1, ∀α ∈ M and α∈M π α = 1, so that (2.44) is equivalent to (2.41). The second difficulty comes from the asymptotic nature of the model. This is a straight consequence of the issue explained in the section 2.3.1. For all V , obtained after transformation of the initial dataset (2.6), the pseudo-angle V V lies in the interior of the central simplex S d . In order to tackle this issue, the idea is to split V between its contribution in α, V α = (V j ) j∈α and its non-asymptotic residual ε α c = (V j ) j / ∈α (figure (2.10)). A hidden variable Z valued in {0, 1} K is added for the inference such that π α k = P(Z k = 1), where α k ∈ M = {α 1 , . . . , α K } and

K k=1 Z k = 1.
A sub-asymptotic model is then proposed, that comprehends the Dirichlet mixture, on which is added a model for the non-asymptotic residuals. Thus, for α k ∈ M and for v ∈ (1, ∞) d , the likelihood is given by:

p(v|z k = 1) = r -|α k |-1 k φ α k (w k |m α k , ν α k ) j / ∈α f ε (v j |λ k ),
(2.45)

where

f ε (•|λ k ) is an exponential law of parameter λ k , r k = j∈α k v j and w k = ( v j r k ) j∈α k .
The inference on the model is then made by using the pseudo-angles W α = Vα Vα and the residuals ε α c through an adapted Expectation-Maximization (EM) algorithm.

CONTRIBUTIONS

Conference articles with proceedings

• Feature clustering for extreme events analysis, with application to extreme stream-flow data. (ECML 2016)

Figure 2.10: Pseudo-angle and residual.

Authors: M. Chiapino, A. Sabourin.

• A multivariate extreme value theory approach to anomaly clustering and visualization. (submitted) Authors: M. Chiapino, A. Sabourin, S. Clémençon, V. Feuillard.

Journal article

• Identifying groups of variables with the potential of being large simultaneously. (submitted) Authors: M. Chiapino, A. Sabourin, J. Segers.

OPEN PROBLEMS

In this thesis, we proposed new methods for infering the dependence structure of high dimensional extreme phenomena. As usual methods for dimension reduction do not apply in this context, we found a way to tackle high dimension through a sparse version of the angular measure Φ. Nevertheless, as Φ is a limit measure, the observations cannot be used straightforwardly for inference. In particular, we have to deal with the different supports between the law of the observed data and the limit model. Under some sparsity assumption, the asymptotic tests carried out in Chapter 4 allow to recover the asymptotic dependence structure in the extremes, in relation to the support of Φ. The decomposition V observed = W R + ε residual is then used to infer the model (Chapter 5). The consistency of the model comes from the fact that lim R→∞ ε residual = 0. However, we did not investigate the rate of convergence of such a limit, and this should constitute new studies. Also, other decompositions are possible such as V observed = (W + ε residual )R, that would lead to new models and algorithms.

Clustering Extreme Features

Abstract The dependence structure of extreme events of multivariate nature plays a special role for risk management applications, in particular in hydrology (flood risk). In a high dimensional context (d > 50), a natural first step is dimension reduction. Analyzing the tails of a dataset requires specific approaches: earlier works have proposed a definition of sparsity adapted for extremes, together with an algorithm detecting such a pattern under strong sparsity assumptions. Given a dataset that exhibits no clear sparsity pattern we propose a clustering algorithm allowing to group together the features that are 'dependent at extreme level', i.e.that are likely to take extreme values simultaneously. To bypass the computational issues that arise when it comes to dealing with possibly O(2 d ) subsets of features, our algorithm exploits the graphical structure stemming from the definition of the clusters, similarly to the Apriori algorithm, which reduces drastically the number of subsets to be screened. Results on simulated and real data show that our method allows a fast recovery of a meaningful summary of the dependence structure of extremes.

INTRODUCTION

Extreme value analysis is of primarily interest in many contexts. One example is the machine learning problem of anomaly detection, where one needs to control the false positive rate in the most remote regions of the sample space (Clifton et al. (2011a); [START_REF] Lee | On-line novelty detection using the kalman filter and extreme value theory[END_REF]; Goix et al. (2015a[START_REF] Goix | Sparse representation of multivariate extremes with applications to anomaly ranking[END_REF]). Another example is the field of environmental sciences, where extreme events (floods, droughts, heavy rainfall, . . . ) are of particular concern to risk management, considering the disastrous impact these events may have. Using Extreme Value Theory (EVT) as a general setting to understand or predict extreme events has a long history [START_REF] Katz | Statistics of extremes in hydrology[END_REF]).

In spatial problems, exhibiting areas (groups of weather stations) which may be concomitantly impacted by severe events is of direct interest for risk management policies. Identifying these groups may also serve as a preliminary dimensionality reduction step before more precise modeling. Before proceeding further, we emphasize that standard dimension reduction techniques such as PCA do not apply to extremes as these methods essentially focus on the data around the mean by analyzing their covariance structure, which does not characterize the behavior of extremes (i.e. data

INTRODUCTION

far away in the tails of the distribution). In the present paper, the quantity of interest is river water-flow recorded at several locations of the French river system. The features of the experiment are thus the stream-flow records at different gauging stations, and the goal is to recover maximal groups of stations where extreme discharge may occur simultaneously. Our dataset consists of daily stream-flow recorded at 92 gauging stations scattered over the French river system, from 1969, January 1 st to 2008, December 31 st . It is the same dataset as in [START_REF] Giuntoli | Low flows in france and their relationship to large-scale climate indices[END_REF] Goix et al. (2015a[START_REF] Goix | Sparse representation of multivariate extremes with applications to anomaly ranking[END_REF] proposed an alternative algorithm to do so with a reduced computational complexity of order O(nd log n). To the best of our knowledge, these are currently the only available examples in the literature to handle the recovery of groups of features which are representative of the extremal dependence structure. (See Section 3.2.2 for a precise definition of the latter). In Goix et al. (2015a[START_REF] Goix | Sparse representation of multivariate extremes with applications to anomaly ranking[END_REF], the extremal dependence structure is called sparse if the number of such groups is small compared with 2 d -1, the total number of groups. The output of Goix et al. (2015a[START_REF] Goix | Sparse representation of multivariate extremes with applications to anomaly ranking[END_REF]'s DAMEX algorithm is a (hopefully sparse) vector M = (μ α , α ⊂ {1, . . . , d}) of size 2 d -1, where μα is a summary of the dependence strength at extreme levels between features j ∈ α. The fact that μα is positive means that the probability that all features in α be large while all others stay small, is not negligible. Various datasets have been analyzed in Goix et al. (2015a[START_REF] Goix | Sparse representation of multivariate extremes with applications to anomaly ranking[END_REF] (wave data from the north sea, standard anomaly detection datasets, simulated data) for which the DAMEX algorithm does exhibit a sparsity pattern, thus pointing to a relatively small number of groups of features α (each being of relatively small size |α| compared to the original dimension of the problem) which could be jointly extreme. However, DAMEX becomes unusable in situations where the subsets of features impacted by extreme events vary from one event to another: DAMEX then finds a very large number of subsets to be dependent, but not significantly so, (i.e.0 < μα 1), so that no sparsity pattern emerges. This is precisely the case with the river flow dataset analyzed in the present paper (see Section 3.5). Contributions. One remarkable aspect of the preliminary analysis of the river flow dataset using DAMEX is the tendency of those many subsets α's such that μα > 0, to form clusters, whose members differ from each other by a single or two features only. In practice, this means that several distinct events have impacted 'almost' the same group (cluster) of stations. The aim of this paper is to propose a methodol-ogy enabling to gather together such 'close-by' feature subsets into feature clusters. This is done by relaxing the constraint that 'features not in α take small values' when constructing the representation of the dependence structure. The output of the CLEF algorithm (CLustering Extreme Features) proposed in the present work (Section 3.4) is an alternative representation which remains usable in this 'weakly sparse' context. This representation can still be explained and understood in the multivariate EVT framework (Section 3.3), as in [START_REF] Chautru | Dimension reduction in multivariate extreme value analysis[END_REF]; Goix et al. (2015a[START_REF] Goix | Sparse representation of multivariate extremes with applications to anomaly ranking[END_REF]. We emphasize that the scope of CLEF algorithm concerns situations similar to the hydrological problem considered here, where the DAMEX algorithm does not yield a readable output. In the opposite case (e.g. with the wave dataset or the anomaly detection datasets analyzed in [START_REF] Goix | Sparse representation of multivariate extremes with applications to anomaly ranking[END_REF]Goix et al. ( , 2015a)), DAMEX remains a better option than CLEF in view of its computational simplicity.

Relationships with Apriori. The dimension reduction problem considered here (determining for which subgroups of features concomitant large values are frequent) is closely related to the problem of frequent itemsets mining, specifically to the well known Apriori algorithm introduced by [START_REF] Agrawal | Fast algorithms for mining association rules[END_REF], see also [START_REF] Gunopulos | Discovering all most specific sentences[END_REF]. Indeed, the present problem can be recast as follows: encoding as a '1' any value above a specified threshold and as a '0' any value below this threshold, one obtains a binary dataset. The goal is now to recover the groups items (features) for which concomitant '1' values are frequent, which is precisely the frequent itemsets mining problem. The combinatorial issue that arises with possibly 2 d -1 subsets is circumvented in Apriori by considering subsets of increasing sizes, letting a subset 'grow' until its frequency in the database is not significant anymore. This incremental principle is also related to a subset clustering method proposed in [START_REF] Agrawal | Automatic subspace clustering of high dimensional data[END_REF]. CLEF proceeds in a similar way to Apriori, the main difference being that CLEF comes with a natural interpretation in terms of multivariate EVT. Also, in practice, the stopping criterion used to decide whether incrementing a feature subset is different in CLEF and in Apriori, allowing CLEF to detect larger groups, as discussed in Sections 3.3 and 4.1.

The paper is organized as follows. Section 3.2 sets up the extremal feature clustering problem and establishes connections with multivariate EVT. The dimension reduction method that we promote is explained in Section 3.3: Section 3.3.1 recalls existing work and points out some limitations, Section 3.3.2 makes explicit the links between the considered problem and the Apriori algorithm. The CLEF algorithm is described in Section 3.4. Section 3.5 gathers results: the output of CLEF is compared with that of DAMEX and Apriori. Section 3.6 concludes. The Python code for CLEF, the scripts and the dataset used for our hydrological case study are available at https://bitbucket.org/mchiapino/clef_algo.

PROBLEM STATEMENT AND MULTIVARIATE EVT VIEWPOINT

PROBLEM STATEMENT AND MULTIVARIATE EVT VIEWPOINT

Formal statement of the problem

Consider a random vector X = (X 1 , . . . , X d ) in R d (here, X j is the water discharge recorded at location j). The first step when it comes to learning dependence properties of X is to standardize the features, in the same spirit as in the copula framework, which allows one to focus only on the dependence structure of X. One popular standardization choice in multivariate EVT is the probability integral transform: Denote by F the joint cumulative distribution function (c.d.f.) of X and by F j the marginal c.d.f. of X j . For simplicity, let us assume that each F j is continuous (no point masses), so that with probability one, 0 < F j (X j ) < 1. The standardized variable used for dependence analysis are

V j = (1 -F j (X j )) -1 , j = 1, . . . , d and V = (V 1 , . . . , V d ).
Doing so, the V j 's are identically distributed according to standard Pareto distribution, P(V j > t) = 1/t, t ≥ 1.

Our goal here is to recover all the maximal subsets of features (stations) α ⊂ {1, . . . , d} which 'may be large together' with non negligible probability. In more formal terms, define the extremal joint excess coefficient,

ρ α := lim t→∞ tP ∀j ∈ α, V j > t = lim t→∞ P ∀j ∈ α, V j > t | V α 1 > t ∈ [0, 1]. (3.1)
The variable t plays the role of a high threshold above which the standardized feature V j is considered as extreme. In practice, estimation will be done by fixing a large t and assuming that the limit in (3.1) is approximately reached. An advantage of the standardization procedure is that a single threshold t is needed to define an extreme event, not d thresholds, since all the features share the same scale. The limit in (3.1) exists under the regularity property (3.3) in the next paragraph. Notice already that the second equality also comes from our standardization choice ensuring that for any j ≤ d, t -1 = P(V j > t) = P(V α 1 > t), which justifies the scaling factor t in the definition. The coefficient ρ α ∈ [0, 1] may be seen as a 'correlation' coefficient for the features X j , j ∈ α at extreme levels. We say that the features {V j , j ∈ α} 'may be large together' if ρ α > 0. One relevant summary of the dependence structure of extremes is thus the set of subgroups

M = {α ⊂ {1, . . . , d} : ρ α > 0}. (3.2)
More precisely, we would like to recover those subgroups α ∈ M which are maximal for inclusion in M, i.e. ∀β such that α β, β / ∈ M. A maximal set of features α ∈ M may be viewed as a cluster, in the sense that every subset β ⊂ α is dependent at extreme level (i.e.ρ β > 0), and that α 'gathers' all of them together. In this paper, a 'cluster' of features is understood as a maximal element α ∈ M.

Connections with multivariate EVT

The working hypothesis in EVT is that, up to marginal standardization, the distribution of X is 'approximately homogeneous' on extreme regions. As pointed out above, if the margins F j are continuous, then the V j 's have the homogeneity property: tP

V j t ≥ x = 1/x, for 1 ≤ j ≤ d, t > 1, x > 0.
The key assumption is that the latter property holds jointly at extreme levels, i.e. that V is jointly regularly varying (see e.g. [START_REF] Resnick | Extreme values, regular variation and point processes[END_REF]), which writes

tP V t ∈ A ---→ t→∞ µ(A), (3.3)
where µ is the so-called exponent measure and where A is any set in R d which is bounded away from 0 and such that µ(∂A) = 0. The exponent measure is finite on any such set A and satisfies, for t > 0,

A ⊂ R d + , tµ(tA) = µ(A)
, where tA = {tx : x ∈ A}. Notice that many commonly used textbook multivariate distributions (e.g. multivariate Gaussian or Student distributions) satisfy (3.3), after standardization to V variables. The measure µ characterizes the distribution of V at extreme levels, since for t large enough (so that the region tA is an 'extreme region' of interest), one may use the approximation P(V ∈ tA) t -1 µ(A). The connection between µ and the ρ α 's is as follows: consider the 'rectangle'

Γ α := {x ∈ R d + : ∀j ∈ α, x j > 1} (3.4)
From the definitions (3.1) and ( 3.3), it follows that ρ α = µ(Γ α ). Thus the family of subset M in (3.2) writes M = {α : µ(Γ α ) > 0}.

Non parametric estimation. In a word, non parametric estimation of extremal characteristics based on i.i.d. data X 1 , . . . , X n (distributed as X) is performed by replacing probability distributions with their empirical counterparts, and by proceeding as if the limit in (3.3) were reached above some large threshold t. Since the F j 's are unknown, set V j i = 1/(1 -F j (X j i )), = 1, . . . , n, j = 1, . . . , d, where F j (x) = n -1 n i=1 1{X j i < x}. Thus V j i ∈ {1, n/(n -1), n/(n -2), . . . , n/2, n} and for each fixed j, and t ≤ n, the number of examples i such that V j i > t is equal to n/t . This suggests that t should be chosen as a function of the sample size and indeed, theoretical guarantees on the estimators are obtained for t = o(n) and t → ∞, e.g. t ∼ √ n , see [START_REF] Beirlant | Statistics of extremes: theory and applications[END_REF], Chapter 3 for more details. After this data preprocessing step, the exponent measure µ of any region

A ⊂ R d + \ {0} is approximated by µ n (A) = t Pn (tA),
where

Pn (A) = n -1 n i=1 δ Vi (A), (3.5)
where δ denotes the Dirac mass. Statistical properties of µ n (or of other functional summaries of it) have been investigated by many authors, see e.g. [START_REF] Qi | Almost sure convergence of the stable tail empirical dependence function in multivariate extreme statistics[END_REF]; Einmahl and Segers (2009); [START_REF] Fougeres | Bias correction in multivariate extremes[END_REF] for the asymptotic behavior, [START_REF] Goix | Learning the dependence structure of rare events: a non-asymptotic study[END_REF] for finite sample error bounds.

DIMENSION REDUCTION FOR MULTIVARIATE EXTREMES

DIMENSION REDUCTION FOR MULTIVARIATE EXTREMES

Existing work

Numerous modeling strategies for low dimensional multivariate extremes (say d ≤ 10) have been proposed, see e.g. [START_REF] Coles | Modeling extreme multivariate events[END_REF] 2013) for semior non-parametric ones. For higher dimensional problems, to this date, the only available dimensionality reduction methods are (to our best knowledge) the recent works [START_REF] Chautru | Dimension reduction in multivariate extreme value analysis[END_REF]; Goix et al. (2015a[START_REF] Goix | Sparse representation of multivariate extremes with applications to anomaly ranking[END_REF]. These three references share the common idea of recovering the sub-cones of R d + on which the exponent measure µ concentrates. The seminal paper on this subject appears to be [START_REF] Chautru | Dimension reduction in multivariate extreme value analysis[END_REF]. It relies on principle nested spheres and spherical k-means and is designed for moderate dimensional problems only (d ≤ 20 in their simulation experiments and d = 4 in their case study) with a relatively simple dependence structure (at most 4 groups of features with extremal dependence, only two for d = 20 in their simulation experiments). The computational burden significantly increases for larger dimensions or more elaborate dependence structures, as discussed in Section 4.4 of the cited reference. In particular the dimensionality of the problem considered in the present paper (d = 92 with up to 53 dependent groups of features) is outside the scope of Chautru ( 2015)'s algorithm.

The present work is mainly related to Goix et al. (2015a[START_REF] Goix | Sparse representation of multivariate extremes with applications to anomaly ranking[END_REF]) insofar as it relies on a simple counting procedure on rectangular regions as in (3.4). As a comparison, Goix et al. (2015a[START_REF] Goix | Sparse representation of multivariate extremes with applications to anomaly ranking[END_REF] consider the truncated cones

C α = {x : x ∞ ≥ 1, x j > 0 for j ∈ α ; x j = 0 for j / ∈ α}. (3.6)
The importance of such cones in the analysis comes from the homogeneity property of µ. More precisely, a subset of features α may take large values together while the others take small values, if and only if µ assigns a positive mass to C α . The approach proposed in [START_REF] Goix | Sparse representation of multivariate extremes with applications to anomaly ranking[END_REF] consists in 'thickening' the cones C α , i.e. defining for some small > 0 (typically, = 0.1),

C α, = {x ∈ R d + : x ∞ ≥ 1 ; x -1 ∞ x j > for j ∈ α ; x -1 ∞ x j ≤ for j / ∈ α}. (3.7) The quantity µ α := µ(C α ) is approximated by its empirical counterpart on C α, , μ(C α ) = µ n (C α, )
, where µ n is the empirical estimator defined in (3.5). In practice a tolerance parameter µ min has to be chosen: for any α such that µ n (C α, ) < µ min , one sets μ(C α ) = 0. The final output of [START_REF] Goix | Sparse representation of multivariate extremes with applications to anomaly ranking[END_REF]'s DAMEX algorithm is the potentially sparse 2 d -1-vector M = (μ α ) α⊂{1,...,d} mentioned in the introduction, with μα := μ(C α ).

One shortcoming of DAMEX is that no sparsity pattern is produced in case of 'noise'. Here, noise is understood as a small variability affecting the groups of features concomitantly impacted by an extreme event. As an example, for the hydrological dataset considered here, geophysics determines the main underlying dependence patterns, i.e. the groups of stations where floods tend to occur simultaneously (such as, say, group α 0 = {1, 2, 3, 4}); however due to meteorological variability, the actual observed floods sometimes affect some neighboring stations 5, 6, so that in the dataset, the observed groups would be e.g. {1, 2, 3, 4, 5}, {1, 2, 3, 4, 6}, {1, 2, 3, 4} . In such a case, the empirical mass is scattered over many sub-cones C α, (three instead of one). This example suggests an alternative approach allowing to gather together those C α, 's that are 'close', as detailed next.

Gathering together 'close-by' cones, incremental strategy

One way to gather different C α, 's together is to relax the condition that 'all the features V j for j / ∈ α take small values' in the definition of C α, . This yields the rectangular region Γ α defined in (3.4). Unlike the regions C α, 's, the Γ α 's do not form a partition of the positive orthant of R d , and indeed the fact that a point V i belongs to Γ α does not tell anything about its features V j i for j / ∈ α. The problem addressed in [START_REF] Goix | Sparse representation of multivariate extremes with applications to anomaly ranking[END_REF] (recovering M := {α : µ(C α ) > 0}) and the relaxed problem considered here (recovering M := {α : ρ α > 0} = {α : µ(Γ α ) > 0}) are different but however related through the maximal elements of M and M, as stated in the following lemma. Recall that α is said to be maximal in M (resp. M) if there is no superset α α in M (resp. M).

Lemma 3.1. For α ⊂ {1, . . . , d}, α is maximal in M ⇔ α is maximal in M.

(3.8)

The proof is deferred to the Appendix. Another important property from an algorithmic perspective is the following: Lemma 3.2. For α ⊂ {1, . . . , d}, if ρ α = 0 then also for all α ⊃ α, ρ α = 0.

The proof is immediate: remind that ρ α = µ(Γ α ) and notice that for α ⊂ α , Γ α ⊂ Γ α .

Apriori-like incremental strategy Lemma 3.2 suggests searching for α's satisfying ρ α > 0 following the Hasse diagram, among α's of increasing size, and stopping the search along a given path as soon as ρ α = 0 for some α. This incremental strategy is also the main ingredient of the Apriori algorithm [START_REF] Agrawal | Fast algorithms for mining association rules[END_REF]), which we recall for convenience: Let I = {item 1 , . . . , item d } be set of items and let T = {t 1 , . . . , t n } be a set of transactions with t i ⊂ I, ∀i ∈ {1, . . . , n}. The frequency of occurrence of the list of items (itemset) α ⊂ I is defined as f α := 1 n 1≤i≤n 1 α⊂t i . Apriori returns the set {α : f α > f min } with f min > 0. It begins with pairs of items and then increments the size of the itemsets at each step. Indeed if f α ≤ f min then all supersets α ⊃ α verify f α ≤ f min as well, which reduces drastically the number of subsets to be tested.

The CLEF algorithm described next proceeds similarly: a concomitant occurrence of threshold excesses {V j i > t for features j ∈ α} can be identified with a 3

EMPIRICAL CRITERION AND IMPLEMENTATION

transaction and the dependence parameter ρ α can be seen as a (rescaled) theoretical frequency. The main difference between CLEF and Apriori concerns the stopping criterion used by CLEF, which involves a ratio between frequencies for a group α and for subgroups β ⊂ α. The idea behind is to allow detection of larger groups, as described in the following section.

EMPIRICAL CRITERION AND IMPLEMENTATION

Conditional criterion for extremal dependence

Considering the relaxed framework where the goal is to recover the set M defined in (3.2), one needs an empirical criterion for testing the condition 'ρ α (= µ(Γ α )) > 0'. One option would be to consider the empirical estimator ρα = µ n (Γ α ) where µ n is defined in (3.5) which would be the (rescaled) counterpart of the empirical frequency f α used in Apriori. Then the stopping criterion would be 'ρ α ≤ ρ min ', with ρ min a user-defined tolerance level. However, since the Γ α 's (for increasing α's ) are nested, the ρ α 's can only decrease with increasing sizes of α. In other words larger groups tend to be less frequent than smaller groups, even dependent ones. Thus in principle, detecting larger groups as well as smaller ones would require the tolerance level ρ min to depend on the size |α| of the considered subgroup, which would result in d -1 tuning parameters instead of one.

The alternative chosen in the present paper is to consider a conditional frequency, the conditioning event for a group α of size s being such that at least s -1 features are large among the s considered ones. Now, there is no reason why the conditional frequency of occurrence should decrease with |α|, so that a single tuning parameter needs to be chosen, without preventing the detection of large groups. In practice, computing conditional frequencies amounts to compare µ n (Γ α ) with µ n (Γ β ), with β ⊂ α. More precisely, let α ⊂ {1, . . . , d} be such that for some j ∈ α, ρ α\{j} > 0. Consider the probability that all the features in α be large given that all of them but at most one are large and call κ α the limiting conditional probability, namely

κ α = lim t→∞ P ∀j ∈ α, V j i > t
P for all but at most one j ∈ α, V j i > t .

(3.9)

In the sequel, κ α is referred to as the conditional dependence coefficient of α. Notice that the limit in (3.9) does exist: indeed, let

∆ α = ∪ j∈α Γ α\{j} = {x ∈ R d + : x ∞ > 1, j∈α 1 x j ≥1 ≥ |α| -1}.
Since by assumption on α, for some j µ(Γ α\{j} ) = ρ α\{j} > 0, in view of (3.3), we have µ(∆ α ) = lim t→∞ tP( for all but at most one j ∈ α, V j i > t) > 0, so that an equivalent definition of κ α is

κ α = lim t→∞ tP ∀j ∈ α, V j i > t lim t→∞ tP for all but at most one j ∈ α, V j i > t = µ(Γ α ) µ(∆ α ) . (3.10)
The idea is now to compare empirical counterparts of κ α -using µ n instead of µ, see (3.5)-with a single fixed tolerance parameter κ min > 0. This amounts to decide that µ n (Γ α ) results from noise if µ n (Γ α )

µ n (∆ α ). Notice that Γ α ⊂ ∆ α , so that the empirical version of κ α is again a conditional probability and thus belongs to [0, 1] whenever µ n (Γ β ) > 0 for some β ⊂ α such that |α \ β| = 1, which is another argument in favor of an incremental strategy.

Algorithm

CLEF (summarized in Algorithm 3) uses the empirical counterpart of the conditional criterion κ α , which depends on a (high) threshold t as in (3.5):

κα,t := µ n (Γ α ) µ n (∆ α ) = n i=1 1 {#{j∈α: V j i >t} = |α|} n i=1 1 {#{j∈α: V j i >t} ≥ |α|-1}
.

(3.11)

For s ≥ 2, families Âs of subsets α of size s are constructed in an incremental way, among a set of candidates A s , as follows: Set Â1 = {{1}, . . . , {d}}, then

A s = α ⊂ {1, . . . , d} : |α| = s, ∀β ⊂ α s.t. |β| = s -1 : β ∈ Âs-1 } Âs = {α ∈ A s : κα,t > κ min .
(3.12)

The procedure stops at step S ≤ d-1 if ÂS+1 = ∅, at which point our estimator of the family M of dependent subsets is M = ∪ S s=1 Âs . Notice that restricting the search to the set of candidates A s ensures that the 'empirical counterpart' of Lemma 3.2 is satisfied, namely α / ∈ M ⇒ ∀β ⊃ α, β / ∈ M. It also avoids division by zero when computing (3.12). The final output of CLEF is the set Mmax of maximal elements of M. Remark 3.3 (Choice of the parameters t and κ min ). The choice of t is a classical bias/variance trade-off: according to standard good practice in EVT (see e.g. [START_REF] Coles | An introduction to statistical modeling of extreme values[END_REF]), t is chosen in a 'stability region' of relevant summaries of the output. Here we consider the cardinal of M and the mean cardinal of maximal subsets α ∈ M. When t is too small, the observed data may not have reached there ultimate regime (the extremal dependence structure characterized by µ in (3.3)), so that the bias

Algorithm 3 CLEF (CLustering Extreme Features)

INPUT: High threshold t, tolerance parameter κ min > 0. STAGE 1: constructing the Âs 's . Initialization: set S = d.

Step 1: Construct the family of extremal-dependent pairs: set Â2 = {i, j} ⊂ {1, . . . , d}, such that κ{i,j} > κ min .

Step 2: If Â2 = ∅, set S = 2; end STAGE 1. Otherwise

• generate candidate triplets A 3 = {i, j, k} ⊂ {1, . . . , d} s.t {i, j}, {i, k}, {j, k} ∈ Â2 }, • set Â3 = α ∈ A 3 s.t. κα > κ min .
. . .

Step s(s ≤ d):

If Âs = ∅, set S = s; end STAGE 1. Otherwise

• generate candidates of size s + 1:

A s+1 = {α ⊂ {1, . . . , d}, |α| = s + 1, α \ {j} ∈ Âs for all j ∈ α},

• set Âs+1 = α ∈ A s+1 such that κα > κ min .

Output: M = ∪ S s=1 Âs . STAGE 2: pruning (keeping maximal α's only) Initialization:

Mmax ← ÂS . for s = (S -1) : 2, for α ∈ Âs , If there is no β ∈ Mmax such that α ⊂ β, Mmax ← Mmax ∪ {α}. Output: Mmax of κα,t may be large. In contrast, for too large values of t, very few excesses are observed so that the sample size of the data used to compute κα,t is very small and the variance becomes too large. To wit, due to our standardization choice it holds that P(V j i > t) = 1/t. Thus for each j ∈ {1, . . . , d}, {i : V j i > t}| 1/t, so that the total number of data points for which at least one feature exceeds t is approximately within the interval [n/t, dn/t]. Results on real and simulated data (Section 3.5.3) bring out such a stability region for the above mentioned output summaries. It is empirically verified on simulated data that this region corresponds to near optimal values of t.

As for the tolerance parameter κ min , it should be chosen according to the context, keeping in mind that κα,t is an empirical conditional probability of a joint threshold excess of all features j ∈ α (given that at least |α| -1 excesses have occurred). κ min is the level above which this probability is considered as non negligible. The higher κ min , the more stringent the condition, the smaller and fewer the discovered groups α. In this work, we set κ min = 0.25. Remark 3.4. [Construction of the candidates A s+1 ] The graphical structure of the groups of features is exploited to construct candidate incremented groups of features. Namely, members of A s+1 are the maximal cliques of size s in the graph (A s , E s ), where

E s = {(α, α ) ∈ A s × A s : |α ∩ α | = s -1}.
The maximal clique problem is typically attacked via the max-clique algorithm [START_REF] Xie | Max-clique: a top-down graph-based approach to frequent pattern mining[END_REF]). In the present work, clique extraction is performed using the function find_clique of the Python package NetworkX, which uses the Bron & Kerbosch [START_REF] Bron | Algorithm 457: Finding all cliques of an undirected graph[END_REF], [START_REF] Tomita | The worst-case time complexity for generating all maximal cliques and computational experiments[END_REF]) algorithm.

RESULTS

The aim of our experiments is threefold. First, CLEF's output on the hydrological data is illustrated and compared with DAMEX's (Section 3.5.1). Second, the respective performances of CLEF, DAMEX and Apriori are compared quantitatively on simulated data (Section 3.5.2). Finally (Section 3.5.3), the question of the threshold choice is investigated: the goal is to verify whether a stability region such as the one mentioned in Remark 3.3 exists and whether it corresponds to optimal performances of CLEF.

Stream-flow data

The output of CLEF for the stream-flow data may be visualized in Figure 3.1 (Execution time: 0.09 s on a recent 4 cores laptop computer). Following the heuristic mentioned in Remark 3.3, the extremal threshold t was fixed to 320, yielding k = 1186 extreme events (time indexes i such that Vi ∞ ≥ t). The parameter κ min was fixed to 0.25. A total number of 53 clusters (elements of Mmax ) are returned by the CLEF algorithm, the size of which varies between 2 and 7. At first inspection, Figure 3.1 agrees with general climatologic facts: in the north-western part of France, the climate is driven by large scale oceanographic perturbations, so that extreme floods tend to impact a large number of gauging stations simultaneously. The south-eastern part of France is rather subject to localized events (e.g. the so-called 'orages Cévenols' in the vicinity of the Mediterranean coast). This yields smaller clusters, both in terms of number of stations and of spatial extent.

As a comparison, Table 3.1 shows the outcome of [START_REF] Goix | Sparse representation of multivariate extremes with applications to anomaly ranking[END_REF]'s DAMEX algorithm with the stream-flow data. These results show that no matter the choice of the thickening parameter in (3.7), the data do not concentrate on 'a few' thickened cones C α, , instead most of the empirical mass is spread onto many of them. In other words, there are too many subcones with positive mass, but not in a significant way.

RESULTS

Figure 3.1: Output of CLEF for the stream-flow dataset: Maximal groups of stations α ∈ M that are likely to be jointly impacted by an extreme event.

Clusters of stations are marked by colored edges between their members, the color scale indicates the number of stations forming the cluster.

Simulation experiments

In order to quantify the relative performances of CLEF, DAMEX and Apriori, we generate d-dimensional datasets under a model such that the exponent measure µ concentrates on p specified cones (C α 1 , . . . , C αp ). Notice that p, (α 1 , . . . , α p ) only determine the generating model, they are not used as inputs of any of the three algorithms compared here. The generated data are 'realistic' in the sense that all the features are positive (the points lie in the interior cone C {1,...,d} ), even though the furthest points in the tails concentrate near the subcones C α k 's. Namely, we use the asymmetric logistic extreme value model [START_REF] Tawn | Modelling multivariate extreme value distributions[END_REF]), from which data is simulated using Algorithm 2.2 in [START_REF] Stephenson | Simulating multivariate extreme value distributions of logistic type[END_REF]. 20 datasets of size n = 100 .10 3 , d = 100, are generated. For each dataset, p subsets α 1 , . . . , α p of {1, . . . , d} are randomly chosen, which sizes follow a truncated geometric distribution (the maximum cluster size is 8). We aim at reproducing the fact that different events associated with a single α usually impact a group of stations which differs from α by a few stations (the impacted area is not deterministic). To this end, for each step i = 1, . . . , n, and each subset α j , j = 1, . . . , p, one randomly chosen 'noisy' feature l i,j ∈ {1, . . . , d} \ α j is added to α j . For CLEF, DAMEX and Apriori algorithms, the extreme threshold parameter t is chosen so that #{i≤n: Vi ∞≥t} n ≈ 5%. Table 3.2 summarizes the average performance of the three algorithms, for p = 40, 50, 60, 70. In these experiments, the CLEF algorithm recovers most of the charged p subsets α 1 , . . . , α p in average, and significantly more than Apriori. In contrast, DAMEX does not recover the sparse structure of the data. It should be noted that in situations (not reported here) where no noisy feature is added, Apriori and DAMEX perform as well as CLEF.

Influence of the threshold choice

The high threshold t plays a decisive role in our framework as it determines which standardized features V j i are considered as extreme. Recall that the estimate V j i is discrete (see Section 3.2.2). A more convenient way to evaluate the influence of the threshold t is thus to consider instead k := # i ∈ {1, . . . , n} : || Vi || ∞ > t the total number of extreme points. Two significant summaries of CLEF output which are the number of clusters | M| and their average sizes 1

| M| α∈ M |α|, are plotted as a function of k. Figure 3.2 (hydrological data) and the first two panels of Figure 3.3 (simulated data) confirm the existence of stability regions (vertical red lines). The simulation experiments show that choosing the parameter in such regions ensures an optimal performance for CLEF, since both match exactly the one of lowest errors. The large width of the stability region for the simulated data (Figure 3.3 compared to Figure 3.2) may be explained by the fact that the generative model is a classical parametric extreme value model for which the asymptotic regime is nearly reached even for small thresholds, leading to large stability regions.

CONCLUSION

We propose a novel dimension reduction method for the analysis of extremes of multivariate datasets via feature clustering. This is done in adequacy with the framework of multivariate extreme value theory. The proposed algorithm makes use of the graphical structure of the problem, scanning the multiple possible subsets of features in a time efficient way. Results on a hydrological stream-flow data and on simulated data demonstrate the relevance of this approach on datasets which would not exhibit any sufficiently sparse structure when analyzed with existing algorithms. Future work will focus on the statistical properties of the empirical criteria κα,t involved in the algorithm, which would allow to analyze the output as a statistical test for independence at extreme levels. Upper panel: number of detected clusters, middle panel: average cluster size, lower panel: number of errors of CLEF (as in table 3.2). Vertical red lines (k ∈ {12000, . . . , 81000}): stability region.

APPENDIX: PROOF OF LEMMA 3.1

Step 1. As a first step we show that M ⊂ M, i.e. µ(C α ) > 0 ⇒ µ(Γ α ) > 0.

Proof. Write C α = >0, ∈Q R α, , where R α, = {x ∈ R d + : x ∞ ≥ 1; x j > (j ∈ α); x i = 0 (i / ∈ α)}. Assume µ(C α ) > 0.
Since µ(C α ) < ∞, by the monotonous limit property of the measure µ, we have µ(C α ) = lim →0 µ(R α, ). Also, from the definitions, R α, ⊂ Γ α . Thus,

µ(C α ) > 0 ⇒ ∃ ∈ (0, 1) : µ(R α, ) > 0 ⇒ µ( Γ α ) > 0 ⇒ ρ α = µ(Γ α ) = µ( Γ α ) > 0.
Step 2. We now prove the reverse inclusion for maximal elements of M, i.e.

α is maximal in M ⇒ α ∈ M. (3.13) Proof. Consider, for i / ∈ α, the set ∆ i, = Γ α ∩ {x ∈ R d + : x i > }, so that Γ α = i∈{1,...,d}\α ∈Q∩(0,1) ∆ i, ∪ R α,1 . Thus, α ∈ M ⇒ µ(Γ α ) > 0 ⇒ ∃i, µ(∆ i, ) > 0 or µ(R α,1 ) > 0 (3.14)

APPENDIX: PROOF OF LEMMA 3.1

To prove (3.13), it is enough to show that

α ∈ M ⇒ for i / ∈ α, µ(∆ i, ) = 0. (3.15)
Indeed if (3.15) is true, and if α ∈ M, then (3.14) implies that µ(R α,1 ) > 0, and the result follows from the inclusion R α,1 ⊂ C α . We show (3.15) by contradiction. If µ(∆ i, ) > 0 for some i / ∈ α, then

1 ∆ i, = 1 Γ α ∩ {x ∈ R d + : x i > 1} ⊂ Γ α∪{i} ,
thus µ(Γ α∪{i} ) > 0, which contradicts the maximality of α in M.

Step 3.

From (3.13), if α is maximal in M then α ∈ M. Now if α is maximal in M but not in M, there exists β α in M. Thus from Step 1, β ∈ M, a contradiction. Hence α is also maximal in M. Conversely, if α is maximal in M then (
Step 1) α ∈ M. If α was not maximal in M, there would exist β α maximal in M, and from (3.13), β ∈ M, contradicting the maximality of α in M.

Asymptotic Tests on the Coefficient of Tail Dependence

Abstract Identifying groups of variables that may be large simultaneously amounts to finding out which joint tail dependence coefficients of a multivariate distribution are positive. The asymptotic distribution of a vector of nonparametric, rank-based estimators of these coefficients justifies a stopping criterion in an algorithm that searches the collection of all possible groups of variables in a systematic way, from smaller groups to larger ones. The issue that the tolerance level in the stopping criterion should depend on the size of the groups is circumvented by the use of a conditional tail dependence coefficient. Alternatively, such stopping criteria can be based on limit distributions of rank-based estimators of the coefficient of tail dependence, quantifying the speed of decay of joint survival functions. Numerical experiments indicate that the algorithm's effectiveness for detecting tail-dependent groups of variables is highest when paired with a criterion based on a Hill-type estimator of the coefficient of tail dependence.

INTRODUCTION

A question that often arises when monitoring several variables is which groups of variables are prone to be large simultaneously. In food risk management, for instance, the variables under consideration may be the concentrations of different contaminants in blood samples of consumers. In environmental applications, one may be interested in several physical variables such as wind speed and precipitation recorded at several locations, with the purpose of setting off a regional warning when several of these variables exceed a high threshold. In the context of semi-supervised anomaly detection, when the training sample is mostly made of normal instances, identifying the groups of variables which are likely to be large together allows to label certain new instances as abnormal. The latter use case is the motivation behind the DAMEX algorithm [START_REF] Goix | Sparse representation of multivariate extremes with applications to anomaly ranking[END_REF][START_REF] Goix | Sparse representation of multivariate extremes with applications to anomaly detection[END_REF]. In a regular variation framework, identifying those groups among d variables that may be large simultaneously amounts to identifying the support of the exponent measure. The algorithm returns the list of groups of features α ⊂ {1, . . . , d} such that the mass of the empirical exponent measure on certain cones exceeds a user-defined threshold. However, when the empirical version of the exponent measure is scattered over a large number of such cones, the DAMEX

INTRODUCTION

algorithm does not discover a clear-cut structure. [START_REF] Chiapino | Feature clustering for extreme events analysis, with application to extreme stream-flow data[END_REF] encounter this difficulty for extreme streamflow data recorded at several locations of the French river system.

To overcome this issue, the same authors come up with the CLEF (CLustering Extreme Features) algorithm. Instead of partitioning the sample space, CLEF considers nested regions corresponding to increasing subsets of components. A group of variables is enlarged until there is no longer enough evidence that all features in it may be large together. In this respect, CLEF resembles the Apriori algorithm [START_REF] Agrawal | Fast algorithms for mining association rules[END_REF], which is a data-mining tool for discovering maximal sets of items among d available items that are frequently bought together by consumers. Apriori considers increasing itemsets that are made to grow until their frequency falls below a user-defined threshold. In CLEF, the stopping criterion concerns the relative frequency of simultaneous occurrences of large values of all components in a considered subset compared to the frequency of simultaneous occurrences of larges values of all but one component in this subset. [START_REF] Chiapino | Feature clustering for extreme events analysis, with application to extreme stream-flow data[END_REF] find the method to work well on real and simulated data but do not investigate the asymptotic properties of the statistic underlying the stopping criterion.

Our contributions are three-fold. First, we investigate the asymptotic behavior of the statistic underlying CLEF. In this way, the informal stopping criterion can be turned into a proper hypothesis test with controllable level. A second issue concerns the specification of the null hypothesis in the CLEF stopping criterion. Originally, a certain conditional tail dependence coefficient, κ α , related to a given group of variables α ⊂ {1, . . . , d} is supposed to be above a strictly positive, user-defined and therefore somewhat arbitrary threshold. We propose instead to base the stopping criterion on the hypothesis that a multivariate version of the coefficient of [START_REF] Ledford | Statistics for near independence in multivariate extreme values[END_REF] and [START_REF] Ramos | A new class of models for bivariate joint tails[END_REF] is equal to one. The test is based on the limit distributions of multivariate extensions of nonparametric estimators in [START_REF] Peng | Estimation of the coefficient of tail dependence in bivariate extremes[END_REF] and [START_REF] Draisma | Tail dependence in independence[END_REF][START_REF] Draisma | Bivariate tail estimation: dependence in asymptotic independence[END_REF]. Third, we conduct a numerical experiment to compare the finite-sample performance of the DAMEX algorithm and the CLEF algorithm with the various stopping criteria. We find that overall, the multivariate extension of the Hill-type estimator in [START_REF] Draisma | Bivariate tail estimation: dependence in asymptotic independence[END_REF] yields the most reliable procedure to detect maximal groups of asymptotically dependent variables.

Section 4.2 casts the problem in the language of regular variation and introduces the tail dependence coefficients upon which the CLEF stopping criteria will be based. Necessary background on empirical tail dependence functions and processes is reviewed in Section 4.3, including a new result for the empirical joint tail function. In Section 4.4, we derive the asymptotic distribution of the statistic used in CLEF and turn the heuristic stopping criterion implemented in [START_REF] Chiapino | Feature clustering for extreme events analysis, with application to extreme stream-flow data[END_REF] into a statistical test with asymptotically controllable level. Two alternative tests based on the asymptotic distributions of estimators of the Ledford-Tawn-Ramos coefficient of tail dependence are constructed in Sections 4.5 and 4.6. We report the results of our simulation experiments in Section 4.7. Section 4.8 concludes. Proofs are gathered in Appendix 4.9 while the pseudo-code for the CLEF algorithm and variations is provided in Appendix 4.10.

REGULAR VARIATION AND TAIL DEPENDENCE COEFFICIENTS

Bold letters denote vectors and binary operations between vectors are understood componentwise. The indicator function of a set A is denoted by 1 A . For t ∈ R∪{∞}, we let t α denote the constant vector of (R ∪ {∞}) α with all coordinates equal to t. In the special case α = {1, . . . , d}, the index α is usually omitted for brevity when clear from the context: for instance, 0 = 0 {1,...,d} = (0, . . . , 0) ∈ R d .

Let X = (X 1 , . . . , X d ) be a random vector in R d with cumulative distribution function F , whose margins F 1 , . . . , F d are continuous. We assume that the transformed vector V = (V 1 , . . . , V d ) with V j = 1/{1 -F j (X j )} for all j ∈ {1, . . . , d} is regularly varying on the cone [0, ∞] d \ {0} with (nonzero) limit or exponent measure µ. This means that µ is finite on Borel sets of [0, ∞] d \ {0} bounded away from the origin and that lim

t→∞ tP[V ∈ tA] = µ(A), (4.1)
for all Borel sets A ⊂ [0, ∞] d \ {0} such that 0 / ∈ ∂A and µ(∂A) = 0. The measure µ is homogeneous, i.e., µ(s • ) = s -1 µ( • ) for all 0 < s < ∞, and therefore assigns no mass to hyperplanes parallel to the coordinate axes. As a consequence, (4.1) applies to finite and infinite rectangles that are bounded away from the origin and whose sides are parallel to the coordinate axes. The measure µ characterizes the extremal dependence structure of X. The reader is referred to Resnick (2007a[START_REF] Resnick | Extreme values, regular variation and point processes[END_REF] for an introduction to regular variation.

Let ∅ = α ⊂ {1, . . . , d}. Particular instances of (4.1) include the extremal coefficient λ α [START_REF] Schlather | A dependence measure for multivariate and spatial extreme values: Properties and inference[END_REF] and the joint tail coefficient ρ α :

λ α = lim t→∞ tP[∃j ∈ α : V j > t] = µ({u ∈ [0, ∞) d | ∃j ∈ α : u j > 1}), (4.2) ρ α = lim t→∞ tP[∀j ∈ α : V j > t] = µ({u ∈ [0, ∞) d | ∀j ∈ α : u j > 1}). (4.3)
In the bivariate case, |α| = 2, and with our choice of Pareto margins, we have

ρ α = lim t→∞ P(V α 1 > t | V α 2 > t)
, the upper tail dependence coefficient denoted by χ in [START_REF] Coles | Dependence measures for extreme value analyses[END_REF].

Our general objective is to propose statistically sound procedures to recover maximal subgroups α of components that are likely to be concomitantly large. Our aim can thus be phrased as recovering the maximal subsets α ⊂ {1, . . . , d} such that ρ α > 0.

Since ρ α ≤ ρ β as soon as α ⊃ β, any positive tolerance level with which we would like to compare an estimate of ρ α should depend on α and in particular be decreasing as a function of the cardinality |α|. To circumvent this issue, [START_REF] Chiapino | Feature clustering for extreme events analysis, with application to extreme stream-flow data[END_REF] consider for α such that |α| ≥ 2 the conditional tail dependence

REGULAR VARIATION AND TAIL DEPENDENCE COEFFICIENTS

coefficient

κ α = lim t→∞ P ∀j ∈ α : V j > t j∈α 1{V j > t} ≥ |α| -1 , (4.4)
which is the limiting conditional probability that all variables in α exceed a large threshold given that all but at most one already do. In contrast to ρ α , the coefficient κ α has no particular reason to decrease as a function of |α|. Note that

ρ α = µ(Γ α ) while κ α = µ(Γ α )/µ(∆ α ) = ρ α /µ(∆ α ) where Γ α = {x ∈ [0, ∞) d | ∀j ∈ α : x j > 1} and ∆ α = {x ∈ [0, ∞) d | j∈α 1 {x j ≥1} ≥ |α| -1}, provided |α| ≥ 2. If µ(∆ α ) = 0, then µ(Γ β ) = 0 for all β ⊂ α with |β| = |α| -1; in that case, we define κ α = 0.
In the CLEF algorithm [START_REF] Chiapino | Feature clustering for extreme events analysis, with application to extreme stream-flow data[END_REF], the criterion to decide whether ρ α > 0 or not is that κ α ≥ C, where C is a user-defined tolerance level, κ α = µ(Γ α )/ µ(∆ α ), and µ is the empirical exponent measure in (4.8) below. The level C can be chosen independently of α. Still, its choice is somewhat arbitrary, and in particular, the user has no control of false positives. In Section 4.4, we will provide the asymptotic distribution of κ α and propose a test statistic with a guaranteed asymptotic level.

If ρ α = 0 (or κ α = 0), the limiting distributions of the statistics √ k( ρ α -ρ α ) and √ k( κ α -κ α ) are degenerate at zero. We therefore have no control on the asymptotic levels of tests based on those statistics under H 0 : κ 0 = 0. This is why will have to define a CLEF stopping criterion in terms of a test of H 0 : κ α ≥ κ min versus H 1 : κ α < κ min instead, in terms of a user-defined level κ min > 0. The choice of κ min is somewhat arbitrary; in the simulation experiments (Section 4.7), we choose κ min = 0.08.

In Sections 4.5 and 4.6, we consider alternative CLEF stopping criteria based on estimators of the coefficient of tail dependence η α ∈ (0, 1]. For bivariate distributions, the coefficient has been introduced by [START_REF] Ledford | Statistics for near independence in multivariate extreme values[END_REF] and extended by [START_REF] Ramos | A new class of models for bivariate joint tails[END_REF] in order to model situations in between asymptotic dependence (ρ {1,2} > 0) and full independence of X 1 and X 2 . De Haan and Zhou (2011) and [START_REF] Eastoe | Modelling the distribution of the cluster maxima of exceedances of subasymptotic thresholds[END_REF] proposed and studied a multivariate extension of η α for |α| ≥ 3. The model assumption is that there exist η α ∈ (0, 1] and a slowly varying function L α such that

P[∀j ∈ α : V j > t] = t -1/ηα L α (t).
(4.5)

Suppose that the limit ρ α in (4.3) exists and that (4.5) holds. Then ρ α > 0 implies η α = 1. The converse is true as well, provided lim inf t→∞ L α (t) > 0. Modulo this side condition, which we will take for granted, the null hypothesis ρ α > 0 corresponds to the simple hypothesis η α = 1.

We will test the null hypothesis η α = 1 via multivariate extensions of nonparametric estimators of η α in [START_REF] Peng | Estimation of the coefficient of tail dependence in bivariate extremes[END_REF] and [START_REF] Draisma | Bivariate tail estimation: dependence in asymptotic independence[END_REF]. The null limit of the test statistic is non-degenerate, so that the asymptotic level of the test can be controlled, with no need to introduce an additional tolerance parameter κ min . The estimators that we will study are related to the Pickands estimator and the Hill estimator for the extreme value index of T α = min j∈α V j , respectively. The maximum likelihood estimator, also considered in [START_REF] Draisma | Bivariate tail estimation: dependence in asymptotic independence[END_REF], is less suitable to our context due to its relative computational complexity, since the test is destined to be performed on a large number of subsets of {1, . . . , d}. See also the review [START_REF] Bacro | Measuring and modelling multivariate and spatial dependence of extremes[END_REF] and the references therein.

Remark 4.1. The DAMEX algorithm [START_REF] Goix | Sparse representation of multivariate extremes with applications to anomaly detection[END_REF] is designed to recover the family M of non-empty subsets α of {1, . . . , d} with the property that

µ x ∈ [0, ∞) d x ∞ ≥ 1; ∀j ∈ α , x j > 0 and ∀j / ∈ α, x j = 0 > 0.
In contrast, our focus is on

M = {α | ρ α > 0} = {α | κ α > 0}.
Still, the maximal elements of M for the inclusion order are also the maximal elements of M (Chiapino and Sabourin, 2016, Lemma 1). The two problems of finding the maximal elements of M or M are thus equivalent.

EMPIRICAL TAIL DEPENDENCE FUNCTIONS AND PROCESSES

To find the asymptotic distribution of nonparametric estimators of the various dependence coefficients, we rely on empirical tail processes. Let the random vector X ∼ F be as in Section 4.2; in particular, assume regular variation as in (4.1) with exponent measure µ. Let Λ be the push-forward measure of µ on [0, ∞] d \ {∞} induced by the transformation x → 1/x = (1/x 1 , . . . , 1/x d ), i.e., Λ(

• ) = µ({x ∈ [0, ∞] d \ {0} | 1/x ∈ • }).
For ∅ = α ⊂ {1, . . . , d}, consider the stable tail dependence function α : [0, ∞) α → [0, ∞) and the joint tail dependence function

r α : [0, ∞] α \{∞ α } → [0, ∞) given by α (x) = lim t→0 t -1 P[∃j ∈ α : F j (X j ) > 1 -tx j ] = Λ({y | ∃j ∈ α : y j < x j }), r α (x) = lim t→0 t -1 P[∀j ∈ α : F j (X j ) > 1 -tx j ] = Λ({y | ∀j ∈ α : y j < x j }). (4.6)
From (4.2) and (4.3), clearly λ α = α (1 α ) and ρ α = r α (1 α ). For brevity, we write = {1,...,d} and r = r {1,...,d} . Note that α

(x) = (xe α ) for x ∈ [0, ∞) α , where e α ∈ {0, 1} d has components e α,j = 1 α (j). Similarly, r α (x) = r(xι α ) for x ∈ [0, ∞] α \ {∞ α },
where ι α ∈ {1, ∞} d denotes the vector such that ι α,j = 1 if j ∈ α and ι α,j = +∞ otherwise. By the inclusion-exclusion formula, for x ∈ [0, ∞) α , writing x β = (x j ) j∈β , we have

r α (x) = ∅ =β⊂α (-1) |β|+1 β (x β ), α (x) = ∅ =β⊂α (-1) |β|+1 r β (x β ). (4.7)
Let X i = (X i,1 , . . . , X i,d ), for i ∈ {1, . . . , n}, be an independent random sample from F , having continuous margins and satisfying (4.1). Let k = k(n) → ∞ as

EMPIRICAL TAIL DEPENDENCE FUNCTIONS AND PROCESSES

n → ∞, while k(n) = o(n). Following for instance [START_REF] Einmahl | An m-estimator for tail dependence in arbitrary dimensions[END_REF]; [START_REF] Goix | Sparse representation of multivariate extremes with applications to anomaly detection[END_REF]; [START_REF] Qi | Almost sure convergence of the stable tail empirical dependence function in multivariate extreme statistics[END_REF], we rely on ranks to obtain an approximately Pareto-distributed sample

V i = ( V i,1 , . . . , V i,d ). Let F j (x) = n -1 n i=1 1 {X i,j
<x} be the (left-continuous) empirical distribution function of component j ∈ {1, . . . , d} and put V i,j = 1/{1 -F j (X i,j )} = n/(n + 1 -R i,j ), where R i,j is the rank of X i,j among X 1,j , . . . , X n,j . The empirical counterparts to µ and Λ are

µ( • ) = 1 k n i=1 δ (k/n) V i ( • ), Λ( • ) = 1 k n i=1 δ (n/k)/ V i ( • ), (4.8)
respectively, with δ a the Dirac measure at the point a. Replacing Λ by Λ in the definition of α and r α produces the empirical tail dependence function

α (x) = k -1 n i=1 1{∃j ∈ α : n + 1 -R i,j ≤ kx j } = k -1 n i=1 1{∃j ∈ α : X i,j ≥ X (n-kx j +1),j }
and the empirical joint tail function

r α (x) = k -1 n i=1 1{∀j ∈ α : n + 1 -R i,j ≤ kx j } (4.9) = k -1 n i=1 1{∀j ∈ α : X i,j ≥ X (n-kx j +1),j },
where X (1),j ≤ . . . ≤ X (n),j are the ascending order statistics of X 1,j , . . . , X n,j and • is the floor function. The identities (4.7) hold for α and r α as well. Einmahl et al. (2012, Theorem 4.6) find the weak limit of the empirical process √ k( -) on [0, T ] d for any T > 0. We leverage their theorem to show a similar result for √ k( r α -r α ), jointly in α. The following conditions stem from the cited article.

Condition 1 (Uniform tail convergence). There exists γ > 0 such that, uniformly in x ∈ [0, 1] d with d j=1 x j = 1, we have

t -1 P[∃j = 1, . . . , d : F j (X j ) > tx j ] -(x) = O(t γ ), t → ∞. Condition 2 (Moderate k). The sequence k = k(n) satisfies k = o(n 2γ/(1+2γ) ) as n → ∞, with γ > 0 as in Condition 1.
Condition 3 (Smoothness). For all j ∈ {1, . . . , d}, the partial derivative

∂ j = ∂ /∂x j exists and is continuous on the set {x ∈ [0, ∞) d | x j > 0}.
Since is convex, it is continuously differentiable Lebesgue almost everywhere (Rockafellar, 1970, Theorem 25.5). Condition 3 is satisfied for many popular max-stable models (logistic, asymmetric logistic, Brown-Resnick) but fails for max-linear models. Under Condition 3, the partial derivative ∂ j r α = ∂r α /∂x j (j ∈ α) exists and is continuous on {x ∈ [0, ∞) α | x j > 0} and satisfies ∂ j r α (x) = β:j∈β⊂α (-1) |β|+1 ∂ j β (x β ), where x β = (x s ) s∈β . Einmahl (1997) and [START_REF] Einmahl | An m-estimator for tail dependence in arbitrary dimensions[END_REF] consider a centered Gaussian process W indexed by the Borel sets of [0, ∞] d \ {∞} bounded away from ∞ with covariance function

E[W (A) W (B)] = Λ(A ∩ B).
(4.10)

Note that W (∅) = 0 almost surely. For ∅ = α ⊂ {1, . . . , d} and x ∈ [0, ∞) α , write

W α (x) = W ({y ∈ [0, ∞] d | ∀j ∈ α : y j < x j }).
We consider weak convergence as in van der Vaart (1998); van der Vaart and Wellner (1996); notation . We work in the metric space ∞ (S) of bounded, real functions f on an arbitrary set S, the metric being the one induced by the supremum norm, f ∞ = sup x∈S |f (x)|; the double use of the symbol should not give rise to any confusion. The proof of the following proposition and of other results in the paper is deferred to Appendix 4.9. Proposition 4.2. Let X i = (X i,1 , . . . , X i,d ), for i ∈ {1, . . . , n}, be an independent random sample from F , having continuous margins and satisfying (4.1).

Let k = k(n) → ∞ as n → ∞, while k(n) = o(n). If Conditions 1, 2 and 3 hold, then, for T > 0, in the product space ∅ =α⊂{1,...,d} ∞ ([0, T ] α ), we have, as n → ∞, the weak convergence √ k { r α (x) -r α (x)} W α (x) - j∈α ∂ j r α (x) W {j} (x j ) = Z α (x). (4.11)

ESTIMATING THE CONDITIONAL TAIL DEPENDENCE COEFFI-CIENT

This section investigates the asymptotic distribution of the empirical conditional dependence coefficient κ α based on the empirical exponent measure µ. This is achieved by re-writing κ α as a function of the empirical joint tail coefficients ρ α , the distribution of which follows from Proposition 4.2. We also propose consistent estimators of the asymptotic variance of κ α . Combining the two yields a test for the null hypothesis κ α ≥ κ min where κ min ∈ (0, 1) is a tolerance level fixed by the user, to be seen as the minimal limiting conditional probability that all components in a random vector exceed a threshold, given that all of them but at most one already do.

Let ∅ = α ⊂ {1, . . . , d} and recall the sets

Γ α = {x ∈ [0, ∞) d | ∀j ∈ α : x j > 1} and, provided α has at least two elements, ∆ α = {x ∈ [0, ∞) d | j∈α 1 {x j ≥1} ≥ |α| -1}. Write α \ j = α \ {j} for j ∈ α.
Since ∆ α is the disjoint union of the sets Γ α\j \ Γ α and Γ α , where j ∈ α, we find, for every Borel measure ν, the equality

ν(∆ α ) = j∈α ν(Γ α\j ) -(|α| -1) ν(Γ α ).
(4.12)

ESTIMATING THE CONDITIONAL TAIL DEPENDENCE COEFFICIENT

Recall ρ α = µ(Γ α ) and κ α = µ(Γ α )/µ(∆ α ) in (4.4). By (4.12) applied to ν = µ, we have

κ α = ρ α j∈α ρ α\j -(|α| -1)ρ α . (4.13)
Recall the joint tail function r α and its nonparametric estimator r α in (4.6) and (4.9), respectively. Since ρ α = r α (1 α ), we define the estimators

ρ α = µ(Γ α ) = r α (1 α ) and, provided |α| ≥ 2, κ α = µ(Γ α ) µ(∆ α ) = ρ α j∈α ρ α\j -(|α| -1) ρ α .
The asymptotic distribution of the vector of empirical joint tail coefficients follows immediately from Proposition 4.2. Write ρα,j = ∂ j r α (1 α ).

Corollary 4.3. In the setting of Proposition 4.2, we have, jointly in ∅ = α ⊂ {1, . . . , d}, the weak convergence

k n ( ρ α -ρ α ) Z α (1 α ) = G α , n → ∞. (4.14)
The limit distribution is centered Gaussian with covariance matrix

E[G α G α ] = ρ α∪α - j∈α ρj,α ρ α ∪{j} - j ∈α ρj ,α ρ α∪{j } + j∈α j ∈α ρj,α ρj ,α ρ {j,j } . (4.15)
The asymptotic distribution of κ α follows from the one of ( ρ β ) β via the delta method. The asymptotic variance involves the partial derivative ∂ j κ α = ∂κ α /∂x j of the function

κ α (x) = r α (x) j∈α r α\j (x α\j ) -(|α| -1)r α (x) (4.16) for x ∈ [0, ∞) α . Note that κ α (1 α ) = κ α . Write κj,α = ∂ j κ α (1 α ).
Proposition 4.4. In the setting of Corollary 4.3, we have, as n → ∞ and jointly in α ⊂ {1, . . . , d} such that |α| ≥ 2 and µ(∆ α ) > 0, the weak convergence

√ k ( κ α -κ α ) µ(∆ α ) -2 j∈α ρ α\j G α -ρ α j∈α G α\j . (4.17)
For a fixed such α, the limit distribution is N (0, σ 2 κ,α ) with

σ 2 κ,α = 1 -κ α )κ α µ(∆ α ) -1 -j∈α κj,α + i∈α j∈α κi,α κj,α ρ {i,j} + κ α j∈α κj,α 1 -µ(∆ α ) -1 ρ α\j . (4.18)
Following [START_REF] Peng | Estimation of the coefficient of tail dependence in bivariate extremes[END_REF], the asymptotic variance σ 2 κ,α in (4.18) can be estimated consistently by estimating the partial derivatives κi,α via finite differencing applied to the empirical version of κ α (x) in (4.16) obtained by replacing r α and r α\j by r α and r α\j , respectively:

κ α (x) = n i=1 1{∀j ∈ α : X i,j ≥ X (n-kx j +1),j } n i=1 1{∃m ∈ α : ∀j ∈ α \ m : X i,j ≥ X (n-kx j +1),j } Define κj,α,n = 1 2k -1/4 κ α (1 α + k -1/4 e j ) -κ α (1 α -k -1/4 e j ) , (4.19)
with e j the canonical unit vector of R α pointing in direction j ∈ α, and put

σ 2 κ,α = 1 -κ α ) κ α µ(∆ α ) -1 -j∈α κj,α,n + i,j∈α κi,α,n κj,α,n ρ {i,j} + κ α j∈α κj,α,n 1 -µ(∆ α ) -1 ρ α\j . (4.20)
Proposition 4.5. Under the conditions of Proposition 4.4, we have

σ 2 κ,α = σ 2 κ,α + o P (1) as n → ∞, so that √ k( κ α -κ α )/ σ κ,α N (0, 1), provided σ 2 κ,α > 0.
The proof relies on the weak convergence of the empirical process

√ k{ κ α ( • ) - κ α ( • )} on [0, T ] α
for any T > 0. This property follows in turn from Proposition 4.2 and the functional delta method.

We consider a tolerance level κ min ∈ (0, 1) under which the tail dependence between components j ∈ α is deemed negligible compared to the one between components j ∈ β α. In other words, we aim at testing H 0 : κ α ≥ κ min . Since κ α = ρ α /µ(∆ α ), the null hypothesis is that ρ α is greater than some level depending on α. Let 0 < δ < 1 be a (small) probability, and consider the test

τ α,n = 1 κ α < κ min + q δ k -1/2 σ κ,α (4.21)
where q δ is the δ-quantile of the standard normal distribution. By Proposition 4.5, if σ κ,α > 0, the test in (4.21) has asymptotic level δ for

H 0 against H 1 : κ α < κ min . If ρ α = 0, then, in Proposition 4.3, we have √ k( ρ α -ρ α ) = o P (1) as n → ∞: indeed, on the one hand, we have √ k( ρ α -ρ α ) = √ k ρ α ≥ 0,
and on the other hand, its limit distribution is centered Gaussian. Likewise, we have

√ k( κ α -κ α ) = o P (1) as n → ∞ in Proposition 4.4 if κ α = 0.
As a consequence, under the simple hypothesis H 0 : ρ α = 0, the asymptotic level of a test based on the asymptotic distribution of √ k( ρ α -ρ α ) or √ k( κ α -κ α ) cannot be controlled. This is why the test in (4.21) concerns the null hypothesis H 0 : κ α ≥ κ min for some κ min > 0 instead. Alternatively, we propose tests based on estimators of the coefficient of tail dependence η α in (4.5). In Sections 4.5 and 4.6, we consider two such estimators, extending the ones of [START_REF] Peng | Estimation of the coefficient of tail dependence in bivariate extremes[END_REF] and [START_REF] Draisma | Bivariate tail estimation: dependence in asymptotic independence[END_REF], respectively, to the multivariate setting.

COEFFICIENT OF TAIL DEPENDENCE: PENG'S ESTIMATOR

COEFFICIENT OF TAIL DEPENDENCE: PENG'S ESTIMATOR

For bivariate distributions, Peng's [START_REF] Peng | Estimation of the coefficient of tail dependence in bivariate extremes[END_REF] estimator of the coefficient of tail dependence η = η {1,2} is based on the property that the curve t → (log t, log

P[V 1 > t, V 2 > t]
) has an affine asymptote with slope -1/η. A similar idea motivates Pickands' [START_REF] Pickands | Statistical inference using extreme order statistics[END_REF] estimator for the extreme value index. Estimating the ordinate of the curve at t = n/k and t = n/(2k) allows to estimate that slope. Under a second-order regular variation condition, [START_REF] Peng | Estimation of the coefficient of tail dependence in bivariate extremes[END_REF] shows that his estimator is asymptotically normal, both if η = 1 and if η < 1. In the former case, the asymptotic variance depends on the tail dependence function and its partial derivatives, which are unknown but may be estimated consistently, thus leading to tests whose asymptotic levels can be controlled.

Let α ⊂ {1, . . . , d} have at least two elements. Recall the empirical joint tail function r α in (4.9). We define the multivariate extension of Peng's [START_REF] Peng | Estimation of the coefficient of tail dependence in bivariate extremes[END_REF] estimator of η α in (4.5) as

η P α = log(2)/ log{ r α (2 α )/ r α (1 α )}. (4.22)
The asymptotic normality of η P α follows from Proposition 4.2 and the delta method.

Proposition 4.6. In the setting of Proposition 4.2, we have, as n → ∞ and jointly in α ⊂ {1, . . . , d} such that |α| ≥ 2 and ρ α > 0, the weak convergence

√ k( η P α -1) -1 2ρ α log 2 {Z α (2 α ) -2Z α (1 α )} .
The right-hand side is a N (0, σ 2 α,P ) random variable with variance

σ 2 α,P = 1 2(ρ α log 2) 2 ρ α -4ρ 2 α + 2 j∈α ρj,α r α (2 α ∧ ι j ) + j∈α j ∈α ρj,α ρj ,α 3ρ {j,j } -2r {j,j } (2, 1) , (4.23)
where ρ {j,j } = r {j,j } (2, 1) = 1 if j = j and where ι j ∈ {1, ∞} α is the vector which all coordinates equal to 1 except for the j-th one which equals ∞, so that

(2 α ∧ ι j ) m = 1 if m ∈ α \ j and (2 α ∧ ι j ) m = 2 if m = j.
By extending the proof of (Peng, 1999, Theorem 2.1), it is also possible to obtain asymptotic normality of η P α in the case ρ α = 0 and η α < 1 in (4.5). This would require a multivariate extension of the second-order regular variation condition in [START_REF] Peng | Estimation of the coefficient of tail dependence in bivariate extremes[END_REF] in the style of Condition 4 below. For the application as a stopping criterion in the CLEF algorithm, we are only interested in the asymptotic distribution of η P α under the hypothesis ρ α > 0, so we do not pursue this idea any further.

As in Proposition 4.4, the asymptotic variance σ 2 α,P in (4.23) involves unknown quantities, all of which we can estimate consistently. For α ⊂ {1, . . . , d} and j ∈ α,

define ρj,α,n = 1 2k -1/4 r α (1 α + k -1/4 e j ) -r α (1 α -k -1/4 e j ) , (4.24)
where e j is the canonical unit vector in R α pointing in dimension j. Define

σ 2 α,P = 1 2( ρ α log 2) 2 ρ α + j∈α ρj,α,n {-4 ρ α + 2 r α (2 α ∧ ι j )} + j∈α j ∈α ρj,α,n ρj ,α,n 3 ρ {j,j } -2 r {j,j } (2, 1) . (4.25)
Proposition 4.7. In the setting of Proposition 4.2, we have

σ 2 α,P = σ 2 α,P + o P (1) as n → ∞, where α ⊂ {1, . . . , d} is such that |α| ≥ 2 and ρ α > 0. If σ 2 α,P > 0, then √ k( η P α -1)/ σ α,P N (0, 1) as n → ∞.
The proof parallels the one of Proposition 4.5 and is omitted for brevity. The main step is to verify that ρj,α,n = ρj,α + o P (1) as n → ∞, which follows from Proposition 4.2.

To test the hypothesis H 0 : ρ α > 0 at significance level δ ∈ (0, 1), we propose

τ α,η P ,n = 1 η P α < 1 -q 1-δ k -1/2 σ α,P , (4.26)
where q 1-δ is the (1 -δ)-quantile of the standard normal distribution. In the setting of Proposition 4.7, the test in (4.26) has asymptotic level δ for H 0 against H 1 : η α < 1.

COEFFICIENT OF TAIL DEPENDENCE: HILL ESTIMATOR

The coefficient of tail dependence η α in (4.5) is the tail index of the random variable

T α = min j∈α V j : the function t → P[T α > t]
is regularly varying at infinity with index -1/η α . A tractable alternative to Peng's estimator for η α is a Hill-type estimator as in [START_REF] Draisma | Tail dependence in independence[END_REF][START_REF] Draisma | Bivariate tail estimation: dependence in asymptotic independence[END_REF]. Replacing the unobservable Pareto variables V i,j by the rank-based versions

V i,j = n/(n + 1 -R ij ) in Section 4.3 yields an approximate sample T i,α = min j∈α V i,j , i = 1, . . . , n,
from the distribution of T α . Let T (1),α ≤ . . . ≤ T (n),α denote the order statistics of T 1,α , . . . , T n,α . The Hill estimator for η α is defined as

η H α = 1 k k i=1 log T (n-i+1),α T (n-k),α .
(4.27)

COEFFICIENT OF TAIL DEPENDENCE: HILL ESTIMATOR

Under the second-order regular variation conditions stated below, the asymptotic normality of η H α follows from [START_REF] Draisma | Bivariate tail estimation: dependence in asymptotic independence[END_REF], proof of Theorem 2.1). The results in the cited reference cover the bivariate case only. In this section, we verify that they remain valid in any dimension d ≥ 2, and we provide the general expression for the asymptotic variance. Put

E α = [0, ∞] α \ {∞ α }.
Condition 4. For each α ⊂ {1, . . . , d} with |α| ≥ 2, there exist functions c α , c 1,α : E α → [0, ∞) such that c 1,α is neither constant nor a multiple of c α , and there exists q 1,α : (0, ∞) → (0, ∞), with q 1,α (t) → 0 as t → 0, such that, for all x ∈ E α , we have

lim t→0 P[∀j ∈ α : 1 -F j (X j ) ≤ tx j ] P[∀j ∈ α : 1 -F j (X j ) ≤ t] -c α (x) q 1,α (t) = c 1,α (x).
Under Condition 4, the function q α (t) = P[∀j ∈ α : 1 -F j (X j ) ≤ t] is regularly varying at 0 with some index 1/η α . Condition 4 implies that the first-order condition (4.5) holds with the same index 1/η α . In addition, c α (1 α ) = 1 and c α is homogeneous of order 1/η α , i.e., c α (tx) = t 1/ηα c α (x) for t > 0, see [START_REF] Draisma | Tail dependence in independence[END_REF][START_REF] Draisma | Bivariate tail estimation: dependence in asymptotic independence[END_REF]. Under the regular variation assumption (4.1), we have ρ α = lim t→0 q α (t)/t, so that, under Condition 4, ρ α > 0 implies η α = 1, as in [START_REF] Draisma | Bivariate tail estimation: dependence in asymptotic independence[END_REF] for the bivariate case. Finally, if [START_REF] Draisma | Bivariate tail estimation: dependence in asymptotic independence[END_REF], our ρ α is denoted by l for α = {1, 2}.

ρ α > 0, then c α (x) = r α (x)/r α (1 α ) = r α (x)/ρ α . Note that in
The asymptotic variance of the Hill estimator (4.27) involves a Gaussian process whose distribution depends on whether ρ α = 0 or ρ α > 0. As in [START_REF] Draisma | Bivariate tail estimation: dependence in asymptotic independence[END_REF], introduce a centered Gaussian process W 1 on E α with covariance function E[W 1 (x) W 1 (y)] = c α (x∧y) for x, y ∈ E α . Recall the stochastic process Z α in (4.11) and the random variable G α = Z α (1 α ) in (4.14).

Proposition 4.8. Let X i = (X i,1 , . . . , X i,d ), for i ∈ {1, . . . , n}, be an independent random sample from F , having continuous margins and satisfying (4.1). , 2, 3, and 4 hold, then, as

Let k = k(n) → ∞ as n → ∞, while k(n) = o(n). If Conditions 1
n → ∞, √ k η H α -η α N (0, σ 2 α,H ), with σ 2 α,H = η 2 α Var{ W (1 α )}, where W (x) = W 1 (x) if ρ α = 0 and W (x) = ρ -1/2 α Z α (x) if ρ α > 0. In particular, if ρ α > 0, we have σ 2 α,H = ρ -1 α Var(G α ) = 1 -2ρ α + ρ -1 α j∈α j ∈α ρj,α ρj ,α ρ {j,j } . (4.28)
The proof of Proposition 4.8 is based on the arguments developed in the proofs of (Draisma et al., 2004, Theorem 2.1), (Drees, 1998b, Theorem 3.2), and (Drees, 1998a, Example 3.1), which we gather in Appendix 4.9.

Again, the unknown terms in (4.28) may be replaced by their empirical counterparts, leading to an asymptotically consistent test. Recall ρj,α,n in (4.24) and define

σ 2 α,H = 1 -2 ρ α + ρ -1 α j∈α j ∈α ρj,α,n ρj ,α,n ρ {j,j } .
The proof of the consistency of the variance estimator follows the same lines as the proofs of Propositions 4.5 and 4.7 and is omitted.

Corollary 4.9. Under the conditions of Proposition 4.8, if ρ α > 0, we have σ 2 α,H = σ 2 α,H + o P (1) as n → ∞ and thus √ k( η H α -1)/ σ α,P N (0, 1), provided σ 2 α,H > 0. We may exploit Corollary 4.9 to test H 0 : ρ α > 0 in the same way as we did by using Peng's estimator in (4.26): at significance level δ ∈ (0, 1), the null hypothesis is rejected in favour of

H 1 : η α < 1 when η H α < 1 -q 1-δ k -1/2 σ α,H . Remark 4.10. The condition σ 2
α,H > 0 in Corollary 4.9 is satisfied whenever 0 < ρ α < 1. Indeed, in (4.28), we have ρ {j,j } ≥ ρ α and ρj,α ρj ,α ≥ 0, whence σ 2

α,H ≥ 1 -2ρ α + (j,j )∈α 2 ρj,α ρj ,α = 1 -2ρ α + ρ 2 α = (1 -ρ α ) 2 .

SIMULATION STUDY

Our aim is to compare the finite sample performance of the various tests proposed in Sections 4.4, 4.5 and 4.6 within the framework of the CLEF algorithm, the pseudocode of which is given in Appendix 4.10. Three variants of the CLEF algorithm are obtained by varying the criterion according to which a subset α is declared as taildependent:

κ α > κ min -q δ σ κ,α / √ k for CLEF-asymptotic; η α,P > 1 -q δ σ α,P / √ k for CLEF-Peng; and η α,H > 1 -q δ σ α,H / √ k for CLEF-Hill.
The original CLEF criterion was κ α > C for some constant C chosen by the user. For completeness, the output of the DAMEX algorithm [START_REF] Goix | Sparse representation of multivariate extremes with applications to anomaly ranking[END_REF] is included in the comparison.

In practice, the dependence tests based on the tail dependence coefficient should not be carried out to the letter when the test statistic is not defined or when its estimated variance is infinite. Thus, in our experiments, CLEF-Peng and CLEF-Hill are modified so as to take into account additional, common-sense stopping criteria. A subset α will not be part of the list returned by the algorithms under the following conditions:

1. Concerning CLEF-Hill, when ρ α = 0, that is, no extreme record impacts all coordinates in α, the estimated variance of the Hill estimator of η α is infinite. Therefore, ρ α = 0 is considered as a stopping criterion in CLEF-Hill.

2. Concerning CLEF-Peng, when r α (2 α ) = r α (1 α ), the Peng estimator (4.22) is ill-defined. Such a case arises when there are very few points in the joint tail within the subspace generated by α. When the estimated derivatives ρj,α,n are close to zero, and when ρ α 1, the estimated variance σ 2 α,P in (4.25) becomes large, preventing rejection of the null hypothesis. To prevent these issues, each of the conditions ρ α < 0.05 and r α (2 α ) = r α (1 α ) are declared as a stopping criterion in CLEF-Peng.

Experimental setting. : CLEF Chiapino and Sabourin (2016) is designed to face situations where DAMEX [START_REF] Goix | Sparse representation of multivariate extremes with applications to anomaly ranking[END_REF] fails to exhibit a clear-cut dependence structure. A major issue reported in [START_REF] Chiapino | Feature clustering for extreme events analysis, with application to extreme stream-flow data[END_REF] for certain 4.7. SIMULATION STUDY hydrological data is the high variability of the groups of features for which large values occur simultaneously. Because of this, the empirical exponent measure µ assigns low mass to any sub-region partitioning the sample space, see Remark 4.1. The empirical finding motivating the latter work is that the various subsets α involved in simultaneous extreme records could nevertheless be clustered, meaning that many of them have a significant intersection, whereas many symmetric differences comprise just a single or at most a few features.

A natural assumption in this context is that a 'true' list of dependent subsets M = {α 1 , . . . , α K } exists such that µ(C α ) > 0 for α ∈ M and that noisy features are involved in each extreme event. Observed large records then concern groups of the kind α = α ∪ {j}, where α ∈ M and j ∈ {1, . . . , d} \ α.

In our experiments, datasets are generated as follows: The dimension is fixed to d = 100. A family of 'true' dependent subsets M = {α 1 , . . . , α K } of cardinality K = 80 is randomly chosen: the subset sizes |α| follow a truncated geometric distribution, with a maximum subset size set to 8. For simplicity, we forbid nested subsets, so α j ⊂ α k whenever j = k. The maximal elements of M = {α ⊂ {1, . . . , d} | ρ α > 0} are then precisely the elements of M, as explained in Remark 4.1. Finally, two different subsets may have at most two features in common.

Once the dependence structure M has been fixed, the data X 1 , . . . , X n are sampled independently from d-dimensional asymmetric logistic distributions Tawn (1990b), using Algorithm 2.2 in [START_REF] Stephenson | Simulating multivariate extreme value distributions of logistic type[END_REF]. The underlying 'true' distribution function is

G(x) = exp - K m=1 j∈αm (|A(j)|x j ) -1/wα m wα m , ( 4.29)
where A(j) = {α ∈ M | j ∈ α} and w αm is a dependence parameter which is set to 0.1 in our simulations. Actually, to mimic the noisy situation described above, each point X i is simulated according to a slightly different version, G i , of G. For each i = 1, . . . , n and k = 1, . . . , K, we randomly select an additional 'noisy feature'

j i,k ∈ {1, . . . , d} \ α k and set α i,k = α k ∪ {j i,k }. Then M i = {α i,1 , . . . , α i,K } is the collection of 'noisy subsets' for X i and G i (x) is as in (4.29) with A(j) replaced by A i (j) = {α ∈ M i | j ∈ α }.
Results. : We generate datasets of size n = 5e4 and n = 1e5. For each sample size, 50 independent datasets are simulated according to the procedure summarized in the preceding paragraph. We compare the average performance of the three proposed versions of CLEF, together with the original CLEF and DAMEX algorithms, for different choices of k and confidence level δ.

Tables 4.1 and 4.2 gather the results for a confidence level δ equal to 0.001 and 0.0001, respectively. In both tables, the results obtained with the original version of CLEF and DAMEX are included in the comparison with an identical choice of tuning parameters, so that the last two lines of the two tables are the same. In CLEF, the threshold C was chosen by trial and error in the interval (0, κ min ), namely C = 0.05. Imposing that C < κ min is intended to reproduce the effect of the variance term Confidence level for the tests: δ = 0.001. Standard deviations over the 50 samples in brackets. Bold face indicates the best performing algorithm on average for a given n and a given choice of k/n, the proportion of extreme data used.

upon the stopping criterion in CLEF-asymptotic. In DAMEX, the 80 subsets with highest empirical mass are retained and the subspace thickening parameter is set to the default value of 0.1, following the guidelines of the authors. Each algorithm produces a list, M, of groups of features α ∈ {1, . . . , d}. This list is to be compared with the one of K = 80 'true' subsets M. The performance of each algorithm is measured in terms of two criteria: the number of 'true' subsets α ∈ M that appear in M (third column of Tables 4.1 and 4.2); the number of 'errors', that is, the subsets α ∈ M that do not belong to M. These can be understood as 'false positives'. Among these errors, we make the distinction between those which are respectively proper subsets (fourth column of Tables 4.1 and 4.2) or proper supersets (fifth column) of some true β ∈ M, and the other errors (sixth column). CLEF-Peng obtains the best overall scores for both values of δ, but as explained above, a special treatment is reserved for the case ρ α ≤ 0.05, and this threshold constitutes an arbitrary tuning parameter, which can impact the performance significantly. On the other hand, CLEF-Hill does not require any other adjustment than for the special case ρ α = 0 and performs nearly as well as CLEF-Peng with δ = 0.0001 and k/n = 0.005. In addition, CLEF-Hill outperforms all the other methods. In particular, CLEF-asymptotic is globally less accurate than CLEF-Peng and CLEF-Hill. This reflects the fact that the null hypothesis in this algorithm involves an arbitrary κ min > 0 fixed by the user. Our own choice κ min = 0.08 was fixed by trial and error, which is straightforward with synthetic data and could also be achieved by cross-validation in a real use case. Finally, as expected, DAMEX obtains very low scores, because it is not designed to handle the addition of noisy features, as explained earlier.

CONCLUSION

In this work, we propose three variants of the CLEF algorithm [START_REF] Chiapino | Feature clustering for extreme events analysis, with application to extreme stream-flow data[END_REF], replacing the heuristic criterion in the original version with a formal test for asymptotic dependence, and this for all possible subsets of features among {1, . . . , d}. As in the original CLEF implementation, only a small proportion of all 2 d -1 subsets has to be examined, while the computational complexity for each such subset is low. Experimental results indicate that the CLEF algorithm is most effective when based on a test constructed from an extension of the Hill estimator [START_REF] Draisma | Bivariate tail estimation: dependence in asymptotic independence[END_REF] of the multivariate coefficient of tail dependence.

The procedure we propose is nonparametric and rank-based. Parametric approaches, based for instance on the nested asymmetric logistic distribution [START_REF] Tawn | Modelling multivariate extreme value distributions[END_REF], could have a greater sensitivity, at the cost of increased model risk and greater computational complexity. We have also assumed that the observations are serially independent; in the contrary case, the asymptotic variances of the various estimator need to be estimated by some form of bootstrap, which, in high dimensions, poses important theoretical and computational challenges; see [START_REF] Bücher | Multiplier bootstrap of tail copulas with applications[END_REF] for the bivariate and serially independent case.

PROOFS

Proof of Proposition 4.2. For ∅ = α ⊂ {1, . . . , d} and x ∈ [0, ∞) α , put

L α (x) = {y ∈ [0, ∞] d | ∃j ∈ α : y j < x j }, R α (x) = {y ∈ [0, ∞] d | ∀j ∈ α : y j < x j }.
If α = {1, . . . , d}, then just write L rather than L {1,...,d} . Note that L α (x) = L(xe α ) with e α = (1 α (j)) d j=1 and that L {j} (x j ) = R {j} (x j ) and thus W (L {j} (x j )) = W {j} (x j ). Einmahl et al. (2012, Theorem 4.6) show that, in the space ∞ ([0, T ] d ) and under Conditions 1, 2 and 3, we have weak convergence

√ k{ (x) -(x)} W (L(x)) - d j=1 ∂ j (x)W {j} (x j )
as n → ∞. Here, we have taken a version of the Gaussian process W such that the trajectories x → W (L(x)) are continuous almost surely.

As in (4.7), we have, for ∅ = α ⊂ {1, . . . , d} and x ∈ [0, ∞) α , the identity

r α (x) = ∅ =β⊂α (-1) |β|+1 (x β e β )
where x β = (x j ) j∈β . Hence, we can view the vector ( √ k( r α -r α )) ∅ =α⊂{1,...,d} as the result of the application to √ k( -) of a bounded linear map from the space 4.9. PROOFS ∞ ([0, T ] d ) to the product space ∅ =α∈{1,...,d} ∞ ([0, T ] α ). By the continuous mapping theorem, we obtain, in the latter space, the weak convergence

√ k { r α (x) -r α (x)} ∅ =β⊂α (-1) |β|+1 W (L β (x β )) -d j=1 ∂ j β (x β )W {j} (x j 1 β (j)) .
Here we used (x β e β ) = β (x β ). The set-indexed process W satisfies the remarkable property that W (A ∪ B) = W (A) + W (B) almost surely whenever A and B are disjoint Borel sets of [0, ∞] d \ {∞} that are bounded away from ∞: indeed, (4.10) implies E[{W (A∪B)-W (A)-W (B)} 2 ] = 0. It follows that the trajectories of W obey the inclusion-exclusion formula, so that, for ∅ = α ⊂ {1, . . . , d} and x ∈ [0, ∞) α , we have, almost surely,

∅ =β⊂α (-1) |β|+1 W (L β (x β )) = ∅ =β⊂α (-1) |β|+1 W j∈β R {j} (x j ) = W j∈α R {j} (x j ) = W (R α (x)) = W α (x).
We can make this hold true almost surely jointly for all such α and x: first, consider points x with rational coordinates only and then consider a version of W by extending W α to points x with general coordinates via continuity. Similarly, since W {j} (0) = W (∅) = 0 almost surely, we have

∅ =β⊂α (-1) |β|+1 d j=1 ∂ j β (x β )W {j} (x j 1 β (j)) = j∈α β:j∈β⊂α ∂ j β (x β )W {j} (x j ) = j∈α ∂r j (x)W {j} (x j ).
We have thus shown weak convergence as stated in (4.11).

Proof of Corollary 4.3. The weak convergence statement (4.14) is a special case of (4.11): set x = 1 α . The covariance formula (4.15) follows from the fact that

E[W α (1 α )W α (1 α )] = Λ({y ∈ [0, ∞] d | ∀i ∈ α ∪ α : y i < 1}) = µ({u ∈ [0, ∞) d | ∀i ∈ α ∪ α : u i > 1}) = ρ α∪α ;
the first equality follows from (4.10) and the last one from (4.3). We obtain (4.15) by expanding G α = Z α (1 α ) using (4.11) and working out E[G α G α ] with the above identity.

Proof of Proposition 4.4. Let α = {α 1 , . . . , α S } ⊂ {1, . . . , d} with S = |α| ≥ 2 and such that µ(∆ α ) > 0. In view of (4.13), we have κ α = g α (θ α ) and κ α = g α ( θ α ) where θ α = (ρ α , ρ α\α 1 , . . . , ρ α\α S ), θ α = ( ρ α , ρ α\α 1 , . . . , ρ α\α S ), and

g α (x 0 , x 1 , . . . , x S ) = x 0 S j=1 x j -(S -1)x 0 , x ∈ [0, ∞) 1+S . (4.30)
Let ∇g α (x) denote the gradient vector of g α evaluated x and let • , • denote the scalar product in Euclidean space. Proposition 4.3 combined with the delta method as in (van der Vaart, 1998, Theorem 3.1) gives, as

n → ∞, √ k( κ α -κ α ) = √ k{g α ( θ α ) -g α (θ α )} = ∇g α (θ α ), √ k( θ α -θ α ) + o P (1) ∇g α (θ α ), (G α , G α\α 1 , . . . , G α\α S ) ,
the weak convergence holding jointly in α by Slutsky's lemma and Proposition 4.3.

The partial derivatives of g α are

∂g ∂x 0 (x) = S j=1 x j { S j=1 x j -(S -1)x 0 } 2 , ∂g ∂x j (x) = -x 0 { S j=1 x j -(S -1)x 0 } 2 , j = 1, . . . , S.
Evaluating these at x = θ α and using j∈α ρ α\j -(S -1)ρ α = µ(∆ α ) as in (4.12) and (4.13), we find that

∇g α (θ α ), (G α , G α\α 1 , . . . , G α\α S ) = µ(∆ α ) -2 j∈α ρ α\j G α -ρ α j∈α G α\j ,
in accordance to the right-hand side in (4.17).

To calculate the asymptotic variance σ 2 κ,α , we introduce a few abbreviations: we write

R β = R β (1 β ) and W ∩ β = W β (1 β ) = W (R β ) for ∅ = β ∈ {1, .
. . , d} and we put W j = W {j} (1) for j = 1, . . . , d, so that G α = W ∩ α -j∈α ρj,α W j . We find

H α = i∈α ρ α\i G α -ρ α i∈α G α\i = i∈α ρ α\i W ∩ α - j∈α ρj,α W j -ρ α i∈α W ∩ α\i - j∈α\i ρj,α\i W j .
From the proof of Proposition 4.2, recall that W (A ∪ B) = W (A) + W (B) almost surely for disjoint Borel sets A and B of [0, ∞] d \ {∞} bounded away from ∞; moreover, for such A and B, the variables W (A) and W (B) are uncorrelated. Since R α\i is the disjoint union of R α and R α\i \ R α , we have therefore

W ∩ α\i = W ∩ α + W (R α\i \ R α ) almost surely.
In addition, i∈α ρ α\i = µ(∆ α ) + (S -1)ρ α by (4.12) applied to ν = µ. As a consequence,

H α = {µ(∆ α ) -ρ α }W ∩ α -ρ α j∈α W (R α\j \ R α ) + j∈α K α,j W j where K α,j = ρ α i∈α\j ρj,α\i - i∈α ρ α\i ρj,α , j ∈ α.
The S + 1 variables W ∩ α = W (R α ) and W (R α\j \ R α ), j ∈ α, are all uncorrelated, since they involve evaluating W at disjoint sets; W j = W (R {j} ) is uncorrelated with

4.9. PROOFS W (R α\j \ R α ), for the same reason. Moreover, E[W ∩ α W j ] = Λ(R α ∩ R {j} ) = Λ(R α ) = ρ α and similarly E[W (R α\i \ R α )W j ] = Λ(R α\i \ R α ) = ρ α\i -ρ α if i, j ∈ α and i = j. Hence Var(H α ) = {µ(∆ α ) -ρ α } 2 ρ α + ρ 2 α j∈α (ρ α\j -ρ α ) + i,j∈α K α,i K α,j ρ {i,j} + {µ(∆ α ) -ρ α }ρ α j∈α K α,j -ρ α j∈α K α,j i∈α\j (ρ α\i -ρ α ). As j∈α (ρ α\j -ρ α ) = µ(∆ α ) -ρ α and i∈α\j (ρ α\i -ρ α ) = µ(∆ α ) -ρ α,j , we get Var(H α ) = {µ(∆ α ) -ρ α }ρ α µ(∆ α ) + j∈α K α,j + i,j∈α K α,i K α,j ρ {i,j} -ρ α j∈α K α,j {µ(∆ α ) -ρ α\j }. (4.31)
Recall κ α (x) in (4.16). We have

∂ ∂x j 1 κ α (x) x=1α = ∂ ∂x j i∈α r α\i (x α\i ) r α (x) = ρ -2 α ρ α i∈α\j ρj,α\i -ρj,α i∈α ρ α\i = ρ -2 α K α,j .
It follows that κj,α = -ρ -2 α K α,j /(1/κ α ) 2 = -K α,j /µ(∆ α ) 2 . By (4.31), we find that σ 2 κ,α = µ(∆ α ) -4 Var(H α ) is equal to the right-hand side of (4.18). Proof of Proposition 4.5. We only need to prove that σ 2 κ,α = σ 2 κ,α + o P (1) as n → ∞. In view of the expressions (4.18) and (4.20) for σ 2 κ,α and σ κ,α , it is enough to show that κj,α,n = κα,j + o P (1), with κj,α,n in (4.19); indeed, Corollary 4.3 already gives consistency of µ(∆ α ) and ρ β . Now since 2 -1 k 1/4 {κ α (1 α + k -1/4 e j ) -κ α (1 αk -1/4 e j )} → κα,j as n → ∞, a sufficient condition is that for some > 0, sup

[1-,2+ ] α k 1/4 κ α (x) -κ α (x) = o P (1), n → ∞. (4.32)
In turn, (4.32) follows from weak convergence of k 1/2 ( κ α -κ α ) as n → ∞ in the space ∞ ([1-ε, 1+ε] α ). In light of the expressions of κ α and κ α in terms of the (empirical) joint tail dependence functions r β and r β , respectively, weak convergence of k 1/2 ( κ ακ α ) follows from Proposition 4.2 and the functional delta method (van der Vaart, 1998, Theorem 20.8). The calculations are similar to the ones for the Euclidean case in the proof of Proposition 4.4; an extra point to be noted is that if α is such that µ(∆ α ) > 0, then the denominator in the definition of κ α (x) in (4.16) is positive for all x in a neighbourhood of 1 α .

Proof of Proposition 4.6. Proposition 4.2 implies, as n → ∞, the weak convergence

√ k{ r α (2 α ) -r α (2 α }, √ k{ r α (1 α ) -r α (1 α } Z α (2 α ), Z α (1 α ) .
Now η P α = g( r α (2 α ), r α (1 α )) and η α = 1 = g(r α (2 α ), r α (1 α )) = g(2ρ α , ρ α ), with g(x, y) = log(2)/ log(x/y); note that the function r α is homogeneous. Since the gradient of g is ∇g(x, y) = log(2)(log(x/y)) -2 (-x -1 , y -1 ), the delta method gives √ k( η P -1)

∇g(2ρ α , ρ α ), Z α (2 α ), Z α (1 α ) = 1 ρ α log 2 (-1/2, 1), Z α (2 α ), Z α (1 α ) = -1 2ρ α log 2 {Z α (2 α ) -2Z α (1 α )}.
The first part of the assertion follows. As for the variance,

Var(Z α (2 α ) -2Z α (1 α )) = Var(Z α (2 α )) + 4 Var(Z α (1 α )) -4 Cov(Z α (2 α ), Z α (1 α )),
The function r α is homogeneous of order 1, so that ∂ j r α is constant along rays, that is, the function 0 < t → ∂ j r α (t x) is constant. Moreover, the measure Λ is homogeneous of order 1 too. In view of (4.10) and (4.11), it follows that Var(Z α (tx)) = t Var(Z α (x)) for t > 0; in particular Var(Z α (2 α ) = 2 Var(Z α (1 α ). Further, ρ α = (dr α (t, . . . , t)/dt) t=1 = j∈α ρj,α and thus

Var(Z α (1 α )) = ρ α -2 j∈α ρj,α ρ α + j∈α j ∈α ρj,α ρj ,α ρ {j,j } = ρ α -2ρ 2 α + j∈α j ∈α ρj,α ρj ,α ρ {j,j } .
The covariance term is

Cov(Z α (2 α ), Z α (1 α )) = ρ α - j∈α ρj,α ρ α - j∈α ρj,α r α (2 α ∧ ι j ) + j∈α j ∈α ρj,α ρj ,α r {j,j } (2, 1),
with 2 α ∧ ι j as explained in the statement of the proposition. Since j∈α ρj,α = ρ α , we can simplify and find

Var(Z α (2 α ) -2Z α (1 α )) = 6 Var(Z α (1 α )) -4 Cov(Z α (2 α ), Z α (1 α )) = 2ρ α -8ρ 2 α + 4 j∈α ρj,α r α (2 α ∧ ι j )
+ j∈α j ∈α ρj,α ρj ,α 6ρ {j,j } -4r {j,j } (2, 1) .

Divide the right-hand side by (2ρ α log 2) 2 to obtain (4.23).

PROOFS

Proof of Proposition 4.8. To alleviate notations, ∅ = α ⊂ {1, . . . , d} is fixed and the subscript α is omitted throughout the proof. Introduce the tail empirical process Q n (t) = T (n-kt ) for 0 < t < n/k. The key is to represent the Hill estimator as a statistical tail functional (Drees, 1998a, Example 3.1)

of Q n , i.e., η H = Θ(Q n ),
where Θ is the map defined for any measurable function z : (0, 1] → R as Θ(z) = 1 0 log + {z(t)/z(1)}dt when the integral is finite and Θ(z) = 0 otherwise. Let z η : t ∈ (0, 1] → t -η denote the quantile function of a standard Pareto distribution with index 1/η; it holds that Θ(z η ) = η. The map Θ is scale invariant, i.e., Θ(tz) = Θ(z), t > 0.

The proof consists of three steps:

1. Introduce a function space D η,h allowing to control Q n (t) and z η (t) as t → 0.

In this space and up to rescaling, Q n -z η converges weakly to a Gaussian process.

2. Show that the map Θ is Hadamard differentiable at z η tangentially to some well chosen subspace of D η,h .

3. Apply the functional delta method to show that η H = Θ(Q n ) is asymptotically normal and compute its asymptotic variance via the Hadamard derivative of Θ.

Step 1. : Let > 0 and h(t) = t 1/2+ , t ∈ [0, 1]. Then h ∈ H, where

H = {z : [0, 1] → R | z continuous, lim t→0 z(t)t -1/2 (log log(1/t)) 1/2 = 0}.
Introduce the function space

D η,h = {z : [0, 1] → R | lim t→0 t η h(t)z(t) = 0 ; t → t η h(t)z(t) ∈ D[0, 1]},
where D[0, 1] is the space of càdlàg functions. Notice that z η ∈ D η,h . Equip D η,h with the seminorm z η,h = sup t∈(0,1] |t η h(t)z(t)|. Let m = nq ← (k/n) , with • the ceil function, so that k/m → ρ; for self-consistency of the present paper, the roles of k and m are reversed compared to the notation in [START_REF] Draisma | Bivariate tail estimation: dependence in asymptotic independence[END_REF].

From (Draisma et al., 2004, Lemma 6.2), we have, for all t 0 > 0, in the space D η,h , the weak convergence

√ k m n Q n -z η ηt -(η+1) W (t) t∈[0,t 0 ] (4.33)
where W (t) = W (t α ), and W is defined as in the statement of Proposition 4.8. Indeed, the process W in the statement from (Draisma et al., 2004, Lemmata 6.1 and 6.2) has same distribution as W 1 (t α ) in the case ρ = 0; recall that our ρ is denoted by l in [START_REF] Draisma | Bivariate tail estimation: dependence in asymptotic independence[END_REF]. Put U i,j = 1 -F j (X i,j ), and let U (1),j ≤ . . . ≤ U (d),j be the order statistics of U 1,j , . . . , U n,j . In the case ρ > 0, W equals in distribution W dra (t α ) where W dra appears in Lemma 6.1 in the cited reference as the limit in distribution (for α = {1, 2}), for x ∈ E α , of

∆ n,k,m (x) = √ k   1 k n i=1 1{∀j ∈ α : U i,j ≤ U ( mx j ),j } -c(x)   = m k →ρ -1/2 √ m   1 m n i=1 1{∀j ∈ α : U i,j ≤ U ( mx j ),j } rn(x) with k replaced by m -r(x) k mρ →1   .
From Proposition 4.2 and Slutsky's Lemma, we have

∆ n,k,m ρ -1/2 Z α in ∞ ([0, 1] α ). Therefore, W dra = ρ -1/2 Z α , as claimed.
Step 2. : The right-hand side of (4.33) belongs to C h,η = {z ∈ D η,h | z is continuous}. To apply the functional delta-method (van der Vaart, 1998, Theorem 20.8), we must verify that the restriction of Θ to Dη,h is Hadamard-differentiable tangentially to C η,h , with derivative Θ , where Dη,h is a subspace of D η,h such that P(Q n ∈ Dη,h ) → 1 as n → ∞; see the remark following Condition 3 in Drees (1998a). Then it will follow from the scale invariance of Θ, the identities Θ(Q n ) = η H and Θ(z η ) = η, and the weak convergence in (4.33) that

√ k η H -η = √ k Θ( m n Q n ) -Θ(z η ) Θ ηt -(η+1) W (t) t∈[0,1] (4.34)
as n → ∞. From (Drees, 1998a, Example 3.1), the restriction of Θ to Dη,h , the subset of functions on D η,h which are positive and non increasing, is indeed Hadamard differentiable; letting ν denote the measure dν(t) = t η dt + d 1 (t), with 1 a point mass at 1, the derivative is

Θ (z) = 1 0 t η z(t)dt -y(1) = [0,1] z(t)dν(t).
Step 3. : The weak limit in (4.34) is thus equal to [0,1] ηt -(η+1) W (t)dν(t). From (Shorack and Wellner, 2009, Proposition 2.2.1), the latter random variable is centered Gaussian with variance

σ 2 = [0,1] 2 η 2 (st) -(η+1) Cov( W (s), W (t))dν(s)dν(t).
By definition of ν and by symmetry of the covariance,

σ 2 /η 2 = 2 1 s=0 s t=0 (st) -1 Cov( W (s), W (t))dtds A -2 1 s=0
Cov( W (s), W (1))s -1 ds B + Var( W (1)).

CLEF ALGORITHM AND VARIANTS

For any s ∈ (0, 1),

s t=0 Cov( W (s), W (t))(st) -1 dt = 1 u=0 Cov( W (s), W (us))(su) -1 du = 1 u=0 Cov( W (1), W (u))(u) -1 du = B.
The penultimate equality follows from Cov( W (λs), W (λt)) = λ Cov( W (s), W (t)) for λ > 0 and s, t ∈ (0, 1]. Therefore A = B and σ 2 = η 2 Var( W (1)), as required.

CLEF ALGORITHM AND VARIANTS

The CLEF algorithm is described at length in [START_REF] Chiapino | Feature clustering for extreme events analysis, with application to extreme stream-flow data[END_REF]. For completeness, its pseudo-code is provided below. The underlying idea is to iteratively construct pairs, triplets, quadruplets. . . of features that are declared 'dependent' whenever κ α ≥ C for some user-defined tolerance level C > 0. Varying this criterion produces three variants of the original algorithm, namely CLEF-Asymptotic, CLEF-Peng, and CLEF-Hill. The pruning stage of the algorithm is the same for all three variants.

Algorithm 4 CLEF (CLustering Extreme Features)

Input: Tolerance parameter κ min > 0. STAGE 1: constructing the collection M of tail-dependent groups.

Step 1: Put Â1 = {{1}, . . . , {d}} and S = 1.

Step s = 2, . . . , d: If Âs-1 = ∅, end STAGE 1. Otherwise:

• Generate candidates of size s:

A s = {α ⊂ {1, . . . , d} : |α| = s and α \ j ∈ Âs-1 for all j ∈ α}.

• Put Âs = α ∈ A s : κα > κ min .

• If Âs = ∅, put S = s. Abstract In a wide variety of situations, anomalies in the behaviour of a complex system, whose health is monitored through the observation of a random vector X = (X 1 , . . . , X d ) valued in R d , correspond to the simultaneous occurence of extreme values for certain subgroups α ⊂ {1, . . . , d} of variables X j . Under the heavy-tail assumption, which is precisely appropriate for modeling these phenomena, a statistical method for identifying such events/subgroups has been recently developed in Goix et al. (2016a), relying on the concept of angular measure in multivariate extreme value theory, which characterizes the dependence structure of the X j 's in the extremes. It is the purpose of this paper to exploit this approach further, by means of a mixture model that permits to describe the distribution of extremal observations and where the anomaly type α is viewed as a latent variable. In particular, the model enables to assign to any such point X a posterior probability for each anomaly type α, defining implicitely a similarity measure between anomalies. A procedure based on the EM algorithm is also proposed here to infer the parameters of the mixture model from a (truncated) training dataset and it is explained at length how the corresponding posterior similarity measure estimates permit to obtain an informative planar representation of anomalies using standard graph-mining tools. The relevance and usefulness of the 2-d visual display thus designed is illustrated on real datasets, in the aeronautics application domain.

Output: M = ∅ if S = 1 and M = S s=2 Âs if S ≥ 2.

INTRODUCTION

Motivated by a wide variety of applications ranging from fraud detection to aviation safety management through the health monitoring of complex networks, data center infrastructure management or food risk analysis, unsupervised anomaly detection is now the subject of much attention in the data science literature, see e.g. D. Gorinevsky (2012); T. [START_REF] Fawcett | Adaptive fraud detection[END_REF]; [START_REF] Viswanathan | Ranking anomalies in data centers[END_REF]. In frequently encountered practical situations and from the viewpoint embraced in this paper, anomalies coincide with rare measurements that are extremes, i.e. located far from central statistics such as the sample mean. In the 1-d setting, numerous statistical techniques for anomaly detection are based on a parametric representation of the tail of the observed univariate probability distribution, relying on extreme value theory

INTRODUCTION

(EVT) (see e.g. [START_REF] Clifton | Novelty detection with multivariate extreme value statistics[END_REF]; Lee and Roberts (2008a); [START_REF] Roberts | Extreme value statistics for novelty detection in biomedical signal processing[END_REF]; [START_REF] Tressou | Bayesian nonparametrics for heavy tailed distribution. application to food risk assessment[END_REF] among others). In (even moderately) large dimensional situations, the modelling task becomes much harder and many nonparametric heuristics for supervised classification have been thus adapted, substituting rarity for labeling, see e.g. [START_REF] Schölkopf | Estimating the support of a high-dimensional distribution[END_REF], [START_REF] Steinwart | A classification framework for anomaly detection[END_REF] or [START_REF] Liu | Isolation Forest[END_REF]. In the unsupervised setting, whereas many dimensionality reduction and visualization techniques, extending the basic linear PCA methodology, accounting for non linearities or increasing robustness for instance (cf [START_REF] Gorban | Principal Manifolds for Data Visualisation and Dimension Reduction[END_REF] and [START_REF] Kriegel | A general framework for increasing the robustness of pca-based correlation clustering algorithms[END_REF]), have been proposed in the statistics and data-mining literature to describe parsimoniously the 'center' of a massive data distribution (see e.g. [START_REF] Naik | Advances in Principal Component Analysis[END_REF] and the references therein), the issue of clustering extremes or outliers is only recently receiving attention, at the instigation of industrial applications such as those mentioned above and because of the increasing availability of extreme observations in databases: generally out-of-sample in the past, extreme values are becoming observable in the Big Data era.

It is the goal of the present article to propose a novel mixture model-based approach for clustering extremes in the multivariate setup, i.e. when the observed random vector X = (X 1 , . . . , X d ) takes its values in the positive orthant of the space R d with d > 1 equipped with the sum-norm (x 1 , . . . , x d ) = 1≤j≤d |x j |: 'extremes' coinciding then with values x such that P[ X > x ] is 'extremely small'. Precisely, it relies on a dimensionality reduction technique of the tail distribution recently introduced in Goix et al. (2017) (see also Goix et al. (2016a)) and referred to as the DAMEX algorithm. Based on multivariate extreme value theory (MEV theory), the latter method may provide a hopefully sparse representation of the support of the angular measure related to the supposedly heavy-tailed distribution of the r.v. X. As the angular measure asymptotically describes the dependence structure of the variables X j in the extremal domain (and, roughly speaking, permit to assign limit probabilities to directions x/ x in the unit sphere along which extreme observations may occur), this statistical procedure identifies the groups α ⊂ {1, . . . , d} of feature indices such that the collection of variables {X j : j ∈ α} may be simultaneously very large, while the others, the X j 's for j / ∈ α, remain small. Groups of this type being in 1-to-1 correspondence with the faces Ω α = {x ∈ R d : x = 1, x j = 0 if j / ∈ α and x j > 0 if j ∈ α} of the unit sphere composing the support of the angular measure. In practice, a sparse representation of the extremal dependence structure is obtained when only a few such groups of variables can be exhibited (compared to 2 d -1) and/or when these groups involve a small number of variables (with respect to d). Here we develop this framework further, in order to propose a (soft) clustering technique in the region of extremes and derive effective 2-d visual displays, sheding light on the structure of anomalies/extremes in sparse situations. By modelling the distribution of extremes as a specific mixture model, where each component generates a different type α of extremes, the Expectation-Maximization algorithm (EM in abbreviated form) permits to partition/cluster the set of extremal data through the statistical recovery of latent observations, as well posterior probability distributions (inducing a soft clustering of the data in a straighforward manner) and, as a by-product, a similarity measure on the set of extremes: the higher the probability that their latent variables are equal, the more similar two extreme observations X and X are considered. The similarity matrix thus obtained naturally defines a weighted graph, whose vertices are the anomalies/extremes observed, paving the way for the use of powerful graph-mining techniques for community detection and visualization, see e.g. [START_REF] Schaeffer | Graph clustering[END_REF] and [START_REF] Hu | Visualizing large graphs[END_REF] as well as the references therein. Beyond its detailed description, the methodology proposed is applied to a real fleet monitoring dataset in the aeronautics domain and shown to provide useful tools for analyzing and interpreting abnormal data.

The paper is structured as follows. Basic concepts of MEV theory are briefly recalled in Section 5.2, together with the DAMEX technique proposed in Goix et al. (2016a[START_REF] Goix | Sparse representation of multivariate extremes with applications to anomaly detection[END_REF] for estimating the (hopefully sparse) support of a heavy-tailed distribution. Section 5.3 introduces the concept of angular measure and details the mixture model we propose to describe the distribution of extreme data, based on DAMEX output, together with the EM algorithm variant we introduce in order to estimate its parameters. It is next explained in Section 5.4 how to exploit the results of this inference method to define a similarity matrix of the extremal data, reflecting a weighted graph structure of the observed anomalies, and apply dedicated community detection and visualization techniques so as to extract meaningful information from the set of extreme observations. The relevance of the approach we promote is finally illustrated by numerical experiments, on synthetic and real data in Section 5.5.

BACKGROUND AND PRELIMINARIES

As a first go, we start with recalling key notions of MEVT, as well a the inference method investigated in Goix et al. (2016a[START_REF] Goix | Sparse representation of multivariate extremes with applications to anomaly detection[END_REF] to estimate its support. Here and throughout, the Dirac mass at any point x is denoed by δ x , the indicator function of any event A by 1{A}. Capital letters generally refer to random quantities whereas lower case ones denote deterministic values. Finally, boldface letters denote vectors as opposed to Roman letters denoting real numbers.

Multivariate extreme value theory

It is the goal of Extreme Value Theory (EVT) to describe phenomena that are not governed by an 'averaging effect' but can be instead significantly impacted by very large movements. By focusing on large quantiles rather than central statistics such as the median or the sample mean, EVT provides models for the unusual rather than the usual and permits to assess the probability of occurence of rare (extreme) events. Application domains are numerous and diverse, including any field related to risk management as finance, insurance, environmental sciences or aeronautics. Risk

BACKGROUND AND PRELIMINARIES

monitoring is a typical use case of EVT. In the univariate setting, typical quantities of interest are high quantiles of a random variable X, i.e. 1 -p quantiles for p → 0. When p is of the order of magnitude or smaller than 1/N , empirical estimates become meaningless. Another issue is the estimation of the probability of an excess over a high threshold u, p u = P(X > u) when few (or none) observations are available above u. In such contexts, EVT essentially consists in using a parametric model (the generalized Pareto distributions) for the tail distribution, which is theoretically justified asymptotically, i.e. when p → 0 or u → ∞. The required assumption is the existence of two sequences {a n , n ≥ 1} and {b n , n ≥ 1}, with a n > 0 and a non-degenerate cumulative distribution function (c.d.f.) G such that

n P X -b n a n ≥ x ---→ n→∞ -log G(x) (5.1)
For all x in the continuity set of G. Notice that this assumption is satified by most textbook distributions, e.g. the normal, exponential, Cauchy, beta, gamma distributions. The reader is referred to [START_REF] Coles | An introduction to statistical modeling of extreme values[END_REF] and the references therein for an introduction to EVT and its applications.

In the multivariate setting, EVT is concerned about the tail behaviour of a ddimensional random variable X = (X 1 , . . . , X d ). The goal is to infer quantities of the kind P[X 1 > x 1 , . . . , X d > x d ] for large x 1 , . . . , x d . A natural first step is to standardize each component so as to work with identically distributed component and focus on the dependence structure. One convenient choice is to use the probability integral transform: For x = (x 1 , . . . , x d ), let F j (x j ) = P[X j ≤ x j ]. Assuming that F j is continuous, the transformed variable

V j = (1 -F j (X j )) -1 follows a Pareto distribution, P[V j > v] = v -1 , v ≥ 1. Consider the the Pareto-tranformed variable V = (V 1 , . . . , V d ). A multivariate analogue of Assumption (5.1) is nP[ V 1 n > v 1 , or . . . , or V d n > v d ] ---→ n→∞ -log G(v) (5.2)
where v = (v 1 , . . . , v d ), v j > 0 and G is a multivariate c.d.f.. Notice that the choice of Pareto margins implies normalizing sequences a n = n, b n = 0 for each component V j and that G has unit Fréchet margins, G j (v) = e -1/v , v > 0. Other standardizations are possible which lead to alternative normalizing sequences and limits. Exponent measure : To understand the right-hand-side of (5.2), the following result (see e.g. [START_REF] Resnick | Extreme Values, Regular Variation, and Point Processes[END_REF][START_REF] Resnick | Heavy-tail phenomena: probabilistic and statistical modeling[END_REF]) is key: there exists a measure µ on E = R d + \ {0} which is finite on any set A such that 0 does not belong to the closure of A,

such that -log G(v) = µ[0, v] c . µ is called the exponent measure. It is homogeneous of order -1, that is µ(tA) = t -1 µ(A), where tA = {tv, v ∈ A}, A ⊂ R + d . Another consequence of (5.2) is that for all A ⊂ R d + such that 0 / ∈ ∂A, tP[V ∈ tA] ---→ t→∞ µ(A). (5.3)
This convergence property applies immediatly to the problem of estimating the probability of reaching a set tA which is far form 0 (i.e. t is large): one may write P(V ∈ tA) ≈ 1 t µ(A), so that estimates of µ automatically provide estimates for such quantities.

In a word, µ may be used to characterize the distributional tail of V.

Support estimation

The goal of this section is to expose the connection between the support of µ and the subsets of components which may assume large values simultaneously.

Sparse support : Consider α ⊂ {1, . . . , d} a nonempty subset of features and the associated truncated cone

C α = x ≥ 0 : x ∞ ≥ 1, x i > 0 for i ∈ α, x i = 0 for i / ∈ α .
(5.4) as illustrated in Fig. 5.1. The family {C α , α ⊂ {1, . . . , d}, α = ∅} defines a partition of R d + \ [0, 1] d which is of particular interest for ou purposes: indeed, µ(C α ) > 0 if and only if the limiting rescaled probability that all feature in α are large while the others are small is non zero, see Remark 5.1.

Remark 5.1. Consider the -thickened rectangles

R α = v ≥ 0, v ∞ ≥ 1, v i > for i ∈ α,
(5.5)

v i ≤ for i / ∈ α , which defines again a partition of R d + \ [0, 1] d for each fixed ≥ 0. Also C α = ∩ >0, ∈Q R α . Thus by upper continuity of µ, µ(C α ) = lim →0 µ(R α ) with µ(R α ) = lim t→∞ tP( V ∞ > t ∀j ∈ α : V j > t , ∀j / ∈ α : V j < t ).
Now the right-hand side of the latter display corresponds to the event that all features in α are large while the other are small.

In the sequel we denote

µ α = µ(C α ), M = α ⊂ {1, . . . , d}, α = ∅, µ α > 0 .
In theory, every µ α may be positive . However a reasonable assumption in a many high dimensional contexts is that µ α = 0 for the vast majority of the 2 d -1 cones C α . In other words, not all combinations of coordinates of V can be large together, so that the support of µ is sparse. DAMEX algorithm : Earlier works (Goix et al. (2016a)) have proposed an algorithm named DAMEX which produces the list of α's such that the empirical 2017) which scale roughly as k -1/2 , where k is the order of magnitude of the number of largest observations used to learn M and the µ α 's. Given a dataset (X i ) i≤n of independent data identically distributed as X, estimation proceeds as follow: first, replace the unkown marginal distributions F j with their empirical counterpart Fj

A MIXTURE MODEL FOR MULTIVARIATE EXTREME VALUES

(x) = 1 n 1{X i,j < x}. Define then Vi,j = (1 -Fj (X i,j )) -1 and Vi = ( Vi,1 , . . . , Vi,d ). Then choose some k n large enough (typically k = O( √ n)) and define μα as the empirical counterpart of µ(R α ) with t replaced with n/k in (5.3), that is μα = 1 k n i=1 1{ Vi ∈ n k R α }.
Notice that the above definition is a variant of the orignal algorithm in Goix et al. (2016a) which uses thickened cones C α instead of R α . However finite sample guarantees in [START_REF] Goix | Sparse representation of multivariate extremes with applications to anomaly detection[END_REF] are obtained using the latter rather than the original C α 's, which is why we prefer using the R α 's.

A MIXTURE MODEL FOR MULTIVARIATE EXTREME VALUES

The idea behind this section is to consider a mixture model for µ (the distribution of the largest instances of the dataset) indexed by α ∈ M, where M is the output of the DAMEX algorithm. Thus each component α ∈ M of the mixture generated instances V such that V j is likely to be large for j ∈ α In this paper we adopt a plug-in approach and identify M (the true support of µ) with M (the output of DAMEX).

For modeling purposes, the homogeneity property µ(t•) = t -1 µ(•) suggests a preliminary decomposition of µ within a (pseudo)-polar coordinates system, as detailed next.

Angular measure

Let us consider the sum-norm v := v 1 +. . .+v d and S

d := {w ∈ R d + : w = 1} the d-dimensional simplex. Introduce the polar transformation T : v → T (v) = (r, w) defined on R d + \ {0}
, where where r = v is the radial component and w = r -1 v is the angular one. Now define an angular measure Φ on S d (see e.g. [START_REF] Resnick | Heavy-tail phenomena: probabilistic and statistical modeling[END_REF] or [START_REF] Beirlant | Statistics of Extremes: Theory and Applications[END_REF] and the references therein):

Φ(A) := µ v : v > 1, v -1 v ∈ A} ,
with A ⊂ S d . Notice that Φ(S d ) < ∞ and that by homogeneity, it may be shown that µ • T -1 (dr, dw) = r -2 drΦ(dw).

(5.6)

In other words the exponent measure µ factorizes into a radial component and an angular component. Setting R = V and W = R -1 V, a consequence of (5.3) is that

P[W ∈ A, R > tr R > t] ---→ t→∞ r -1 Φ(S d ) -1 Φ(A) (5.7)
for all measurable set A ⊂ S d and r > 1. In other words, given that the radius R is large, the radius R and the angle W are approximately independent, the distribution of W is approximately the angular measure -up to a normalizing constant Φ(S d )and R follows approximately a Pareto distribution. The transformation to unit Pareto margins and the choice of the sum-norm yield the following moment constraint on Φ:

S d w i Φ(dw) = 1, for i = 1, . . . , d.
(5.8)

In addition, the normalizing constant is explicit:

Φ(S d ) = S d Φ(dw) = S d (w 1 + . . . + w d )Φ(dw) = d.
(5.9)

Remark 5.2. The choice of the sum-norm here is somewhat arbitrary. Any other norm on R d for the pseudo-polar transformation is equally possible, leading to alternative moment constraints and normalizing constants. The advantage of the sumnorm is that it allows convenient probabilistic modeling of the angular component w on the unit simplex.

A mixture model

The partition of R d + into cones C α introduced in Section 5.2.2 induces a partition of

S d into 2 d -1 sub-simplices S α , ∅ = α ⊂ {1, . . . , d}, S α = v ∈ R d + : v = 1, v i > 0 for i ∈ α, v i = 0 for i / ∈ α, .
Also by homogeneity, the following equivalence holds:

Φ(S α ) > 0 ⇔ µ(C α ) > 0.
(5.10)

A MIXTURE MODEL FOR MULTIVARIATE EXTREME VALUES

Recall that our key assumption in this work is that the support of µ is sparse, namely we assume that |M| 2 d , where M = {α : µ(C α ) > 0} = {α : Φ(S α ) > 0}. In view of (5.10) and (5.9), the angular measure admits the decomposition

Φ( • ) = d α∈M π α Φ α ( • ) (5.11)
where Φ α is a probability measure on S α and α∈M π α = 1. We make the simplifying assumption that the sets α ∈ M are not nested, i.e. there does not exist two subsets α, β ∈ M such that α ⊂ β. Notice that this assumption could be omitted at the price of additional notational complexity. Introduce the set of coordinates which are singletons in M , E = {j ∈ {1, . . . , d} : {j} ∈ M}, as opposed to M 2 = {α ∈ M : |α| ≥ 2}. Up to relabeling we may assume that E = {1, . . . , d 1 } for some 1 ≤ d 1 ≤ d} or that E = ∅, in which case d 1 = 0. Then α∈M 2 α = {d 1 + 1, . . . , d}. For convenience let us write M 2 = {α 1 , . . . , α K } with K = |M 2 | and let us relabel the weights as π k = π α k for k ≤ K, π K+j = π {j} for j ≤ d 1 . Equipped with these notations, (5.11) becomes

d -1 Φ( • ) = K k=1 π k Φ α k ( • ) + j≤d 1 π K+j δ e j ( • ) (5.12)
where δ a is the Dirac mass at point a and e j = (0, . . . , 1, . . . , 0) the j th canonical basis vector of R d . The singletons weights derive immediately from the moment constraint (5.8): for

i ≤ d 1 , d -1 = K k=1 S {i} w i π k Φ α k (dw) =0 + j≤d 1 S {i} w i π K+j δ e j (dw) = π K+i We obtain Φ( • ) = d K k=1 π k Φ α k ( • ) + j≤d 1 δ e j ( • )
where the vector π ∈ [0, 1] K+d 1 must satisfy

K k=1 π k = 1 -d 1 /d. (5.13)
As is usual for mixture modeling purposes, we introduce a latent variable Z S {1,...,d} , the interior of S d (This is also true using the empirical versions Vi defined in Section 5.2.2).

= (Z 1 , . . . , Z K , Z K+1 , . . . , Z K+d 1 ) such that for k ≤ K (resp. k > K) , Z k = 1 if W has been generated by the mixture component Φ α k (resp. δ e k-K ) and Z k = 0 otherwise. Then for k ≤ K, P[Z k = 1] = π k while for k > K, P[Z k = 1] = d -1 .

A MIXTURE MODEL FOR MULTIVARIATE EXTREME VALUES

To account for the non-asymptotic nature of the data, we model the deviation from the asymptotic support of V, which is α∈M C α , as a noise ε with light tailed distribution, namely an exponential distribution. We denote by Ṽ = V + ε the resulting noisy vector and we assume that only Ṽ is observed (not V). This subasymptotic model may be described as follows Model 2 (Sub-asymptotic mixture model).

1. Let Ṽ be an observed random vector which marginal distributions are approximately Pareto (typically Ṽj = (1 -Fj (X j )) for Fj an estimate of the marginal distribution F j of X j ) 2. Let Z ∈ {0, 1} K+d 1 be a hidden variable as in Model 1 and let R = Ṽ . Then for 1 ≤ k ≤ K + d 1 , conditionally to { R > r 0 , Z k = 1}, the observed vector Ṽ decomposes as

Ṽ = V k + ε k = R k W k + ε k ,
(5.17)

where V k ∈ C α k , ε k ∈ C ⊥ α k are independent from each other and where

R k = V k , W k = R -1 k V k ∈ S α k are as in Model 1, i.e. R k is Pareto distributed, W k ∼ Φ k and R k , W k are independent.
3. The noise's components are independent and identically distributed as a translated exponential distribution with rate λ k : for j ∈ {1, . . . , d} \ α k ,

ε j ∼ 1 + Exp(λ k )
The unknown parameters are those inherited from Model 1, with the addition of the exponential rates λ = (λ k , 1 ≤ k ≤ K + d 1 ) where λ k > 0. The next paragraph describes an EM-algorithm for maximum-likelihood estimation of Model 2.

An EM algorithm for model inference

The likelihood for Model 2, p(ṽ|m, ν, π, λ), for one observation ṽ ∈ (1, ∞) d , ṽ ≥ r 0 , follows directly from the model specification, where f ε ( • |λ k ) denotes the marginal density for the noise ε k given the noise parameter λ k . As specified in Model 2, in this paper we set f ε (x|λ k ) = λ k e -λ k (x-1) ,

x > 1 (a translated exponential density), but any other light tailed distribution could be used instead. Notice that the term r

-|α k |-1 k = r -2 k r -|α k |+1 k
is the product of the radial Pareto density and the Jacobian term for the change of variables

T k : V k → (R k , W k ).
Recall that the constraints are

ν k > 0 (1 ≤ k ≤ K) , λ k > 0 (1 ≤ k ≤ K + d 1 ),
(5.19)

A MIXTURE MODEL FOR MULTIVARIATE EXTREME VALUES

and that π = (π 1 , . . . , π K ) and m = (m 1 , . . . , m K ) satisfy (5.13) and (5.15). The latter linear constraint on (π, m) implies that m and π cannot be optimized independently, which would complicate the M-step of an EM-algorithm. Thus we begin with a re-parametrization of the model ensuring that the moment constraint (5.8) is automatically satisfied.

Re-parametrization of the moment constraint

The main idea behind the re-parametrization is to work with the parameter ρ k,j = π k m k,j instead of (π k , m k,j ). Namely, define a K × (d -d 1 ) matrix ρ = (ρ 1 , . . . , ρ K ) where ρ k,j > 0 for j ∈ α k and ρ k,j = 0 otherwise. Then for all k ∈ {1, . . . , K}, set

π k := j∈α k ρ k,j m k,j := ρ k,j π k , ∀j ∈ α k .
(5.20)

Then (5.13) and (5.

EM algorithm

We summarize below the EM algorithm in our framework. Let n 0 ≤ n be the number of observations Ṽi such that Ṽi > r 0 . To alleviate notations, we may relabel the indices i so that these observations are Ṽ1:n 0 = ( Ṽ1 , . . . , Ṽn 0 ). Let Z i = (Z i,1 , . . . , Z i,K+d 1 ), i ≤ n 0 be the hidden variables associated with Ṽ1:n 0 . In the sequel, θ denotes the set of parameters of the re-parameterized version of Model 2, that is θ = (ρ, ν, λ), and let Θ be the parameter space, that is the set of θ's such that constraints (5.19) and (5.21) hold. Also let p(ṽ|θ, z k = 1) denote the conditional density of Ṽ given (Z k = 1, θ). In view of the likelihood (5.18), it is given for for k ≤ K by • Input Extreme standardized data Ṽ1:n 0 ,

• Initialization Choose a starting value for θ (See Remark 5.3)

• Repeat until convergence:

-E-step: compute for 1 ≤ i ≤ n 0 and k ≤ K + d 1 ,

γ i,k = P[Z i,k = 1 Ṽi , θ]
according to (5.25). Set γ = (γ i,k ) i≤n 0 ,k≤K+d 1 .

-M-step: Solve the optimization problem max θ∈Θ Q(θ, γ)where Q is a lower bound for the likelihood, namely

Q(θ, γ) = n 0 i=1 K+d 1 k=1 γ i,k log π k + log p( Ṽi |θ, z i,k = 1) ,
with π k = P(Z i,k = 1|θ), i.e. (5.24) and where p( Ṽi |θ, z i,k = 1) is given by (5.22) and (5.23).

π k =    ∈α k ρ k,l for 1 ≤ k ≤ K , d -1 for K < k ≤ K + d 1 ,
Denote by θ the solution and set θ = θ .

Remark 5.3. In this work, the output of DAMEX is used for choosing the initial value for ρ. Namely, given M 2 we compute the empirical means m k,j := 1 n 0 n 0 i=1 Ṽi,j for all j in α k and k in {1, . . . , K} and we set π 0 = . . . = π K = 1 K so that we get the corresponding ρ by ρ k,j = π k m k,j . Although it is not likely to verify (5.21), we can easily project ρ on Θ :

ρ init k,j = ρ k,j d K h=1 ρ h,j
.

We now describe at length the E-step and the M-step. E-step. : The γ i,k 's are obtained using the Bayes formula, for 1

≤ k ≤ K + d 1 , γ i,k = p(Z i,k = 1| Ṽi θ) = π k p( Ṽi |z i,k = 1, θ)
1≤ ≤K+d 1 =k π p( Ṽi |z i, = 1, θ) , (5.25)

where π k is defined in (5.24) and p( Ṽi |Z i,k = 1, θ) is given by (5.22) and (5.23). M-step. : Here optimization of Q(θ, γ) with respect ot θ = (ρ, ν, λ) is performed under constraints (5.19), (5.21). Since Q decomposes into a function of (ρ, ν) and a

GRAPH-BASED VISUALIZATION TOOLS

function of λ, and since the constraints on ρ, ν and λ are independent, maximization can be performed separtely over the two blocks. Indeed, gathering terms not depending on θ into a constant C,

Q(θ, γ) = n i=1 K k=1 γ i,k log π k • • • + log ϕ k (W i,k |ρ k , ν k ) + l∈α c k log f ε ( Ṽi,l |λ k ) • • • + K+d 1 k=K+1 γ i,k =k log f ε ( Ṽi,l |λ k ) + C = Q 1 (ρ, ν) + Q 2 (λ) + C,
where

Q 1 (ρ, ν) = n i=1 K k=1 γ i,k log l∈α k ρ k l + log ϕ k (W i,k |ρ k , ν k ) Q 2 (λ) = n i=1 K+d 1 k=1 γ i,k l∈α c k log f ε ( Ṽi,l |λ k ) .
Here we have denoted α k = {k -K} for K < k ≤ K + d 1 , in accordance with the notations from Section 5.3.2. Notice that the dependence of Q 1 and Q 2 on γ is omitted for the sake of concision.

With The function Q 1 being non-concave we use the python package mystic [START_REF] Mckerns | Building a framework for predictive science[END_REF]) to maximize it. For our choice of translated exponential noise, f ε (v|λ k ) = λ k e -λ k (v-1) , v ≥ 1, the maximizer of Q 2 has an explicit expression,

λ * k = |α c k | n i=1 γ i,k n i=1 γ i,k l∈α c k ( Ṽi, -1) , k ≤ K + d 1 .
Remark 5.4. Let γ t and θ t be the results of the t-th iteration of the algorithm then we conclude the iterative process if Q(θ t , γ t ) < Q(θ t-1 , γ t-1 ) + , with a small threshold.

GRAPH-BASED VISUALIZATION TOOLS

In this section, we explain that, beyond the hard clustering that may be straightforwardly deduced from the computation of the likeliest values z 1 , . . . , z n 0 for the hidden variables given the ṽi 's and the parameter estimates produced by the algorithm described in subsection 5.3.3, the statistical model previously introduced defines a natural structure of undirected weighted graph on the set of observed extremes, which interpretable layouts (graph drawing) can be directly derived from, using classical solutions. Indeed, a partition (hard clustering) of the set of (standardized) anomalies/extremes ṽ1 , . . . , ṽn 0 is obtained by assigning membership of each ṽi in a cluster (or cone/sub-simplex ) determined by the component of the estimated mixture model from which it arises with highest probability: precisely, one then considers that the abnormal observation ṽi is in the cluster indexed by

k i = arg max k∈{1, ..., K} γ i,k
and is of type α k i . However, our model-based approach brings much more information and the vector of posterior probabilities (γ i,1 , . . . , γ i,k ) output by the algorithm actually defines soft membership and represent the uncertainty in whether anomaly ṽi is in a certain cluster. It additionally induces a similarity measure between the anomalies: the higher the probability that two extreme values arise from the same component of the mixture model, the more similar they are considered. Hence, consider the undirected graph whose vertices, indexed by i = 1, . . . , n 0 , correspond to the extremal observations ṽ1 , . . . , ṽn 0 and whose edgeweights are w θ (ṽ i , ṽj ), 1 ≤ i = j ≤ n 0 , where w θ (ṽ i , ṽj ) = P Z i = Z j | Ṽi = ṽi , Ṽj = ṽj , θ = K k=1 γ i,k γ j,k .

Graph visualization techniques (see e.g. [START_REF] Hu | Visualizing large graphs[END_REF]), possibly combined with (spectral) graph clustering methods (see e.g. [START_REF] Schaeffer | Graph clustering[END_REF]) when the number n 0 of anomalies to be analyzed is large, can then be used to produce informative layouts. Discussing the merits and limitations of the wide variety of approaches documented in the literature in this purpose is beyond the scope of this paper. It is the goal of the next section to simply illustrate the usefulness of the weighted graph representation of the set of anomalies proposed above, when applying to it state-of-the-art graph-mining tools.

ILLUSTRATIVE EXPERIMENTS

Experiments on simulated data

To assess the performance of the proposed estimator of the dependence structure and of the EM algorithm, we generate synthetic data according to Model 2. The dimensionality is fixed to d = 100 and the mixture components, that is the elements of M = {α 1 , . . . , α K } ∪ E are randomly chosen in the power Set of {1, . . . , d}, with K = 50. The coefficients of the matrix ρ which determines the weights and centers through (5.20) is also randomly chosen, then its columns are normalized so that the 5.5. ILLUSTRATIVE EXPERIMENTS moment constraint (5.21) is satisfied. Finally we fix ν k = 20, 1 ≤ k ≤ K and λ k , 1 ≤ k ≤ K + d 1 are respectively set to (1., 0.75, 0.5, 0.25, 0.1) in the different experiments to illustrate different levels of noise. Then each point Ṽi = R i W i + ε i , i ≤ n is generated with probability π k , k ∈ {1, . . . , K} according to the mixture component k (k ≤ K), that is The threshold r 0 above which points are considered as extreme is fixed to 100. On this toy example, the pre-processing step consisting in applying DAMEX for recovering M produces an exact estimate, so that M = M. Then the procedure described in Algorithm 1 is applied.

R i ∼ P areto(1)|{R i > r 0 } W i ∼ Φ k ε i,j ∼ 1 + Exp(λ k ),
Tables 5.1 and 5.2 show the average absolute errors for the estimates ρ, ν and λ on 50 datasets of the n 0 generated extreme points, for n 0 = 1e + 3, 2e + 3, namely The performance in terms of cluster identification is measured as follows: for each point ṽi , we compare the the true label y i ∈ {1, . . . , K + d 1 } with the label obtained via assignment to the highest probable component, that is ŷi = arg max k∈{1, ..., K+d 1 } γ i,k . Table 5.3 shows the average number of labeling errors for different values of n 0 and λ k . Figure 5.4 illustrates the potential of the proposed approach in terms of visual display of anomalies. A testing set of size 100 is simulated as above, and the corresponding matrix γ is computed according to (5.25) with θ taken as the output of the training step (i.e. Algorithm 1 run with the training dataset of n 0 = 2e3 points). Finally an adjacency matrix w θ (ṽ i , ṽj ) is obtained as detailed in Section 5.4, on which we apply the spectral clustering in order to group the points according to the similarities given by w. Graph visualization of w is performed using the python package 'Networkx [START_REF] Hagber | Exploring network structure, dynamics, and function using NetworkX[END_REF], that provides a spring layout of the graph according to the Fruchtermen-Reingold algorithm [START_REF] Fruchterman | Graph drawing by force-directed placement[END_REF]. A hard thresholding on the edges in w is applied in order to improve readability: edges (i, j) such that w θ (ṽ i , ṽj ) < with a small threshold are removed. Each cluster obtained by spectral clustering is idendified with a specific color for the corresponding nodes.

Flights clustering and visualization

At Airbus, this algorithm is currently being tested to assist the building of health indicators in the context of condition based maintenance. Health indicators are then used for assessing the current state of some system and also for forecasting the future states and future degradation (ex : bleed, power systems, engine, APU, . . . ). Airlines are then informed that some systems should be maintained so as to avoid any operational event in a given time horizon (ex : such as delays, operational interruptions etc . . . ). The building of an health indicator can be basically summarized as follow :

ILLUSTRATIVE EXPERIMENTS

1. Collect health and usage data from various aircrafts (generally one has to consider some similar ones).

2. Collect some operational events happening on these aircrafts due to some aircraft systems errors (ex : operational interruption, delays, etc ...).

3. Identify anomalies in the the health and usage data.

4. Identify some dependencies between health and usage data anomalies and operational events (thanks to statistical tests but also thanks to human expertise).

5. If some dependencies are well identified, then one can quite easily build a health indicator.

One of the tricky part is the identification and the understanding of the anomalies. Indeed different operational events are often recorded corresponding to the degradation of different systems. Usually, a first step of anomalies is identified followed by a clustering of these anomalies for helping in the interpretation. One benefit of the approach proposed in this paper is that it directly provides some similarity measures associated to the anomalies. This strategy is illustrated by Figure 5.5. The proposed method was applied on a dataset of 18553 flights, each of which is characterized by 82 parameters. In order to differentiate between anomalies corresponding to unusual large and small values, each feature is duplicated and each copy of a given feature is defined as the positive (resp. negative) value of the parameter above (resp. below) its mean value.

CONCLUSION

Because extreme values (viewed as anomalies here) cannot be summarized by simple meaningful summary statistics such as local means or modes/centroids, clustering and dimensionality reduction techniques for such abnormal observations must be of very different nature than those developed for analyzing data lying in high probability regions. This paper is a first attempt to design a methodology fully dedicated to the clustering and visualization of anomalies, by means of a statistical mixture model for multivariate extremes that can be interpreted as a noisy version of the angular measure, which distribution on the unit sphere exhaustively describes the limit dependence structure of the extremes (up to a standardization). Localization of the mixture components is understood here as closeness of the data arising from them to a specific sub-simplex forming the support of the angular measure. Considering synthetic and real datasets, we also provided empirical evidence of the usefulness of (graph-based) techniques that can be straightforwardly implemented from the framework we developed.

2. 1 v

 1 Cones for d = 3, x = max d j=1 x j . . . . . . . . . . . . . . . . . . . . 2.2 Distribution of the empirical mass of µ on the cones C α for θ α ∈ {0.1, 0.5, 0.7, 0.9} (top-down, left-right). . . . . . . . . . . . . . . . . . 2.3 Hasse diagram with 4 components. . . . . . . . . . . . . . . . . . . . 2.4 Value of γα on the elements of M F and M 0 . . . . . . . . . . . . . . . 2.5 Value of κα on the elements of M F and M 0 . . . . . . . . . . . . . . . 2.6 Representation of M for the algorithm DAMEX with µ min = 0.002. . 2.7 Representation of M for the criterion {γ α > γ min } with γ min = 0.2. . 2.8 Representation of M for the algorithm CLEF with κ min = 0.2. . . . . 2.9 Values of ηα -(1 -q 1-δ σα √ k ) on the elements of M F and M 0 . . . . . . 2.10 Pseudo-angle and residual. . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Output of CLEF for the stream-flow dataset: Maximal groups of stations α ∈ M that are likely to be jointly impacted by an extreme event. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Stability region for k (number of extreme points) on the stream-flow data. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Stability region for k (number of extreme points) on simulated data. . 5.1 Truncated cones C α in 3D . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Bivariate illustration of Model 1. . . . . . . . . . . . . . . . . . . . . 5.3 Trivariate illustation of the sub-asymptotic model 2 . . . . . . . . . . 5.4 Spectral clustering visualization of a synthetic anomaly test data of size 100 with d = 20 and |M| = 12 . . . . . . . . . . . . . . . . . . . . 5.5 Spectral clustering visualization of flights anomalies with agglomerated Nodes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.6 Spectral clustering visualization of flights anomalies . . . . . . . . . . List of Tables 3.1 Output of Goix et al. (2016b)'s DAMEX algorithm with the hydrological dataset. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Average number of errors (non recovered and falsely discovered clusters) of CLEF, Apriori and DAMEX with simulated, noisy data. . . 4.1 Average number of recovered clusters and errors of CLEF-asymptotic (κ min = 0.08), CLEF-Peng, CLEF-Hill, CLEF and DAMEX on 50 datasets. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Same setting as Table 4.1 with δ = 0.0001 . . . . . . . . . . . . . . . 5.1 Average error on the model parameters, n 0 = 1e3 . . . . . . . . . . . 5.2 Average error on the model parameters, n 0 = 2e3 . . . . . . . . . . . 5.3 Average number of labeling errors . . . . . . . . . . . . . . . . . . . . vii

  La contrainte (1.13) est donc une condition nécessaire pour que Φ soit une mesure angulaire, et c'est aussi une condition suffisante (voir la proposition 5.11 dans[START_REF] Resnick | Extreme Values, Regular Variation, and Point Processes[END_REF]). En conséquence directe de la contrainte des moments, dans le cas où • est la norme L 1 , i.e. x = |x 1 | + . . . + |x d |, nous avons Φ(S d ) = d. En effet:

Figure 1

 1 Figure 1.1: Cônes pour d = 3, x = max d j=1 x j

  Notons que pour ∈ (0, 1) l'ensemble des R α forme toujours une partition de [0, 1] c (pour la norme infinie). De plus, supposons µ(C α ) > 0, en vertu de la remarque 3 le sous-espace C ⊂ C α , minimal pour l'inclusion et tel que µ(C) > 0, est un cône. Or tout cône inclu dans C α est nécessairement d'intersection non vide avec R α . De ce fait, les cônes élargis R α vérifient µ(C α ) > 0 ⇒ µ(R α ) ≥ µ(C) > 0 et sont de mesure de Lebesgue non nulle.

Figure 1

 1 Figure1.2: Répartition de la masse empirique de µ sur les cônes C α pour θ α ∈ {0.1, 0.5, 0.7, 0.9} (de haut en bas, de gauche à droite).

Figure 1 . 3 :

 13 Figure 1.3: Diagramme de Hasse pour 4 composantes.

  Figure 1.4: Valeurs de γα pour les sous-groupes de composantes de M F et de M 0 .

Figure 1

 1 Figure 1.5: Valeurs de κα pour les sous-groupes de composantes de M F et de M 0 .

Figure 1 . 6 :

 16 Figure 1.6: Représentation de M pour l'algorithme DAMEX avec µ min = 0.002.

Figure 1 . 7 :

 17 Figure 1.7: Représentation de M pour le critère {γ α > γ min } avec γ min = 0.2.

Figure 1

 1 Figure 1.9: Valeurs de ηα -(1 -q 1-δ σα √ k ) pour les sous-groupes de composantes de M F et de M 0 .

  sur le modèle de mélange est alors opérée pour les pseudo-angles W α = Vα Vα et les résidus ε α c via un algorithme espérance-maximisation (EM) adapté.

Figure 1 .

 1 Figure 1.10: Pseudo-angle et résidu.

Figure 2

 2 Figure 2.1: Cones for d = 3, x = max d j=1 x j

Figure 2

 2 Figure 2.2: Distribution of the empirical mass of µ on the cones C α for θ α ∈ {0.1, 0.5, 0.7, 0.9} (top-down, left-right).

Figure 2

 2 Figure 2.3: Hasse diagram with 4 components.

  6), (2.7) and (2.8) respectively show the results of the DAMEX algorithm, and the CLEF algorithm with the empirical mass (2.31) and with κα (2.32), for the respective thresholds µ min = 0.002, γ min = 0.2 and κ min = 0.2.Figure (2.6) displays M and the figures (2.7) and (2.8) display M. Each subgroup of dependent coordinates is represented by the convex hull of the corresponding stations. Each color corresponds to a specific size of subgroup. We see that unlike the empirical mass (2.31) and the DAMEX algorithm, the criterion {κ α ≥ κ min } allows to reach large subgroups of stations without being disturbed by small subgroups of stations.

Figure 2

 2 Figure 2.4: Value of γα on the elements of M F and M 0 .

Figure 2

 2 Figure 2.5: Value of κα on the elements of M F and M 0 .

Figure 2 . 6 :

 26 Figure 2.6: Representation of M for the algorithm DAMEX with µ min = 0.002.

Figure 2 . 7 :

 27 Figure 2.7: Representation of M for the criterion {γ α > γ min } with γ min = 0.2.

Figure 2 . 8 :

 28 Figure 2.8: Representation of M for the algorithm CLEF with κ min = 0.2.

  Figure 2.9: Values of ηα -(1 -q 1-δ σα √ k ) on the elements of M F and M 0 .

;

  [START_REF] Cooley | The pairwise beta distribution: A flexible parametric multivariate model for extremes[END_REF];[START_REF] Sabourin | Bayesian model averaging for multivariate extremes[END_REF] for parametric modeling,Boldi and Davison (2007a);[START_REF] Guillotte | Non-parametric bayesian inference on bivariate extremes[END_REF];[START_REF] Sabourin | Bayesian dirichlet mixture model for multivariate extremes: A re-parametrization[END_REF];Fougeres et al. (

Figure 3

 3 Figure 3.2: Stability region for k (number of extreme points) on the stream-flow data. Upper panel: number of detected clusters, lower panel: average cluster size. Vertical red lines (k ∈ {1000, . . . , 1200} / t ∈ [320, 400]): stability region.

Figure 3

 3 Figure 3.3: Stability region for k (number of extreme points) on simulated data.

  STAGE 2: pruning, keeping maximal groups α only.If S = 1, then M max = ∅. Otherwise: Initialization: M max ← ÂS . for s = (S -1) : 2, for α ∈ Âs , If there is no β ∈ M max such that α ⊂ β, then M max ← M max ∪ {α}.

Figure 5 . 1 :

 51 Figure 5.1: Truncated cones C α in 3D

Figure 5

 5 Figure 5.2: Bivariate illustration of Model 1. Here, V is generated from component α k = {1, 2} and the model has two components: M = {{2, }, {1, 2}}. The red line is the Dirichlet density ϕ k on the unit simplex S {1,2} . The red point represent the other component δ e 2 .

Figure 5 .

 5 Figure 5.3 illustrates Model 2 in dimension d = 3.

Figure 5 . 3 :

 53 Figure 5.3: Trivariate illustation of the sub-asymptotic model 2 Here the observed point Ṽ has been generated by component α k = {1, 2}. The grey triangle is the unit simplex, the shaded red area represents the Dirichlet density ϕ k .

  p(ṽ|m, ν, π, λ) =r 0 K k=1 π k r -|α k |-1 k ϕ k (w k |m k , ν k ) j∈α c k f ε (ṽ j |λ k )

  15) together are equivalent to {k:j∈α k } ρ k,j = 1 d , ∀j ∈ {d 1 + 1, . . . , d}. (5.21) In the sequel we denote respectively by p(ṽ|ρ, ν, λ) := p(ṽ|π, m, ν, λ) and ϕ k (w|ρ k , ν k ) := ϕ k (w|m k , ν k ) the likelihood and the Dirichlet densities in the reparameterized model, where (m, π) are obtained from ρ via (5.20).

  p(ṽ|z k = 1, θ) = r -|α k |-1 k ϕ k (w k |ρ k , ν k ) j∈α c k f ε (ṽ j |λ k ) (5.22) and for K < k ≤ K + d 1 , p(ṽ|z k = 1, θ) = ṽ-2 k j∈{1,...,d}\k f ε (ṽ j |λ k ) (5.23)Algorithm 1 (EM algorithm for Model 2).

|

  λ k -λ k |On this toy example, estimation of the means and weights, as well as the noises parameters are almost exact. The estimator for the ν k 's is not so precise, but as shown next, this drawback does not jeopardize clusters identification. Table 5.1: Average error on the model parameters, n 0 = 1e3 λ k = 1. λ k = 0.75 λ k = 0.5 λ k = 0.25 λ k = 0Average error on the model parameters, n 0 = 2e3 λ k = 1. λ k = 0.75 λ k = 0.5 λ k = 0.25 λ k = 0

  Table 5.3: Average number of labeling errorsλ k = 1. λ k = 0.75 λ k = 0.5 λ k = 0.25 λ k = 0

Figure 5 . 4 :

 54 Figure 5.4: Spectral clustering visualization of a synthetic anomaly test data of size 100 with d = 20 and |M| = 12 Each point is represented as a numbered node. The number is the true label, while the colors indicate the clusters obtained by spectral clustering. The spatial arrangement of the nodes is obtained by the Fruchtermen-Reingold algorithm.

  Figure 5.5: Spectral clustering visualization of flights anomalies with agglomerated Nodes The agglomerated visualization is obtained via spectral clustering: each node represents a cluster. Levels of blue show the intern connectivity between the original nodes so that darker clusters have strongly connected elements. The size of each node is proportional to the number of points forming the cluster.

Figure 5 .Figure 5 . 6 :

 556 Figure 5.5 and figure 5.6 display the clustering of 300 'extremal' flights into 18 groups, showing on the one hand the output of the spectral clustering applied to the similarity graph w θ and on the other hand the underlying graph obtained with the same procedure as in Figure 5.4.

ESTIMATION DU SUPPORT DE LA MESURE ANGULAIRE faut

  parcourir peut être beaucoup plus grand, i.e. O(2 d 0 -1) pour d 0 ≤ d. C'est une conséquence du fait que l'union ∪ α⊂{1,...,d} Γ α ne forme plus une partition de l'espace considéré. Plus généralement, cette problématique exploratoire se traduit de la façon suivante: nous voulons trouver tous les sous-groupes I d'un certain ensemble d'indices T tel que I vérifie une condition particulière C(I). La propriété essentielle étant que pour tout I 1 ⊂ I 2 alors C(I 2 ) implique C(I 1 ). Ce problème est directement lié à l'algorithme 'Apriori' introduit dans

	1.4.2.1 Algorithme 'Apriori'

Cette approche a une contrepartie: elle implique une explosion combinatoire. En effet, à la différence du nombre de cônes C α à estimer, qui est au maximum le nombre de points extrêmes du jeu de données, le nombre de rectangles imbriqués Γ α qu'il 1 1.4.

. ESTIMATION DU SUPPORT DE LA MESURE ANGULAIRE

	1	1.4Algorithm 1 CLEF (CLustering Extreme Features) Entrée: Seuil κ min > 0.
			grands | Tous sont grands sauf au plus un ]
							(1.33)
		L'estimateur κα de la quantité (1.32) est donc:
		κα =	μ(	γα	Γ β )	.	(1.34)
				β⊂α
			|β|=|α|-1
		L'algorithme permettant l'estimation de M est décrit par la procédure (1).

PHASE 1: Construire l'ensemble M 0 de tous les sous-groupes asymp- totiquement dépendant. Étape 1: M

  2 0 = α ⊂ {1, . . . , d} : |α| = 2, κα > κ min , S = 2.

Étape s = 3, . . . , d: Si M s-1 0 = ∅, finir PHASE 1. Sinon:

  (en bleue). La ligne bleu représente max α∈M F γα (resp. max α∈M F κα ) et la ligne rouge représente min α∈M 0 γα (resp. min α∈M 0 κα ). Les sous-groupes rouges sont classés par taille croissante, et l'on remarque que la masse empirique (1.31) décroit avec la taille des sous-groupes. De plus, comme max α∈M F γα > min α∈M 0 γα , il est impossible de trouver un seuil γ min pour lequel on peut retrouver M.

		.36)
	Si (1.36) est vérifié, il suffit de choisir un seuil κ min appartenant à l'intervalle
	( max α∈M F κα , min α∈M 0	κα ) pour retrouver M. Les figures (1.4) et (1.5) montrent la valeur
	de la masse empirique γα (1.31) et de κα (1.32) sur les ensembles M 0 (en rouge)
	et M F Nous comparons à présent l'algorithme DAMEX ainsi que l'algorithme CLEF
	pour les différents critères C(α) = {γ α > γ min } et C(α) = {κ α > κ min } sur des
	données journalières de débit fluvial. Le jeu de données provient de d = 92 stations
	réparties sur les différents fleuves de France, enregistré entre le 1er janvier 1969 et
	le 31 décembre 2008, et comprend n = 14610 vecteurs (X 1 , . . . , X n ) à valeur dans
	R d + . L'objectif est de trouver les sous-groupes de stations susceptibles de subir une

, 1 1.5. MODÈLE PARAMÉTRIQUE POUR LA MESURE ANGULAIRE

														13			nb stations : 2 nb stations : 3 nb stations : 4 nb stations : 5 nb stations : 6
	17 18 1920	21	15 22 23 16 24 26 25 48 14	47 49	45 46	42 44	43 34 41 64	33 63 39 40	11 35 36 37 38 61 62	10 32 60 59	8 9 30 31 57 58	27 29 28	5 6 7 12	2 80 68 67 74 75 84 85 82 83 81 86 69 4 76 70 71 72 77 78 79 73	3	0 1
											55 53	56 52			87 88	89
											54
			66	65		51								91
							50			90

). 1 1.5. MODÈLE PARAMÉTRIQUE POUR LA MESURE ANGULAIRE

	18 1920 17	21	16 22 23 24 25 26 15 48 14	49 47	45 46	42 44	34 41 64 43	33 63 39 40	11 35 36 37 38 61 62	10 32 60 59	8 9 30 31 57 58	27 29 28 13	5 6 7 12	2 80 82 68 83 81 74 84 85 86 69 76 77 78 75 70 71 72 79 73 67 4	3	0 1 nb stations : 2 nb stations : 3 nb stations : 4 nb stations : 5 nb stations : 6 nb stations : 7 nb stations : 8 nb stations : 9
														55 53	56 52	87 88	89
														54
				66	65		51								91
								50			90
	Figure 1.8: Représentation de M pour l'algorithme CLEF avec κ min = 0.2.
	Nous pouvons donc opérer la décomposition:
			1 Φ(S d )	Φ(dw) =	α∈M	π α φ α (w α )dw α	(1.39)
	avec les poids π α ∈ (0, 1) tels que α π α = 1 et φ α une densité non nulle sur S α . Le
	modèle de mélange, afin de correspondre à une mesure angulaire valide, doit vérifier
	la contrainte des moments (1.13). Les densités φ α ainsi que les poids π α doivent
	donc vérifier la contrainte:
	α∈M	π α	Sα	w j φ α (w α )dw α =	1 d		, ∀j ∈ {1, . . . , d}.	(1.40)
	j∈α
	La loi de Dirichlet étant une distribution naturelle sur le simplexe, nous proposons
	comme modèle paramétrique un mélange de Dirichlet. Cette loi peut être paramétrée
	par sa moyenne m α ∈ S α et un paramètre de concentration ν α > 0. Nous avons,
	pour tout w ∈ S α :

  , up to 220 gauging stations presenting missing or censored records, which have been removed from our analysis, which results in n = 14610 vectors X 1 , .., X n in R d , with d = 92 the number of stations. The reader is referred to[START_REF] Giuntoli | Low flows in france and their relationship to large-scale climate indices[END_REF] for more details. Related work. Dimensionality reduction for extreme value analysis has emerged very recently in the literature. As far as we know, the seminal contribution is[START_REF] Chautru | Dimension reduction in multivariate extreme value analysis[END_REF] and is restricted to moderate dimensional settings (d ≤ 20, see Section 3.3.1 for more details). The methodology proposed by[START_REF] Chautru | Dimension reduction in multivariate extreme value analysis[END_REF] allows to recover groups of components (features) which may take large values simultaneously, while the other features stay small. For the purpose of anomaly detection,

Table 3 .

 3 1: Output of[START_REF] Goix | Sparse representation of multivariate extremes with applications to anomaly ranking[END_REF]'s DAMEX algorithm with the hydrological dataset. indicate respectively the number of thickened cones C α, with non zero empirical mass, and the percentage of cones (among those such that µ n (C α, ) > 0) containing less than 1% of the 'extreme data', that is of #{i : Vi ∞ > t}.

		# {α : µ n (C α, ) > 0} % α :	#{i:t -1 V i ∈Cα, } #{i: V i ≥t}	< 1%
	0.01	740	100%
	0.05	688	98%
	0.1	639	94%
	0.2	559	88%
	Columns 1 and 2

Table 3 .

 3 2: Average number of errors (non recovered and falsely discovered clusters) of CLEF, Apriori and DAMEX with simulated, noisy data.

	p # errors CLEF # errors Apriori # errors DAMEX
	40	1.2	6.4	72.2
	50	3.5	10.9	91.0
	60	6.3	14.6	112.4
	70	10.1	25.8	134.0

Table 4 .

 4 1: Average number of recovered clusters and errors of CLEF-asymptotic (κ min = 0.08), CLEF-Peng, CLEF-Hill, CLEF and DAMEX on 50 datasets.

	n = 5e4	k/n	recovered	subset errors superset errors other errors
	CLEF-asymptotic 0.003	71.1 (3.0)	7.4 (4.7)	5.1 (2.1)	28.0 (13.3)
		0.005	73.0 (3.7)	8.0 (6.3)	2.4 (1.7)	14.6 (8.9)
	CLEF-Peng	0.003 79.70 (0.7) 1.00 (2.5)	0. (0.)	3.9 (2.7)
		0.005 79.98 (0.1) 0.06 (0.4)	0. (0.)	0.9 (0.9)
	CLEF-Hill	0.003	79.0 (1.4)	2.4 (3.5)	0.04 (0.2)	17.9 (7.0)
		0.005	75.7 (2.4)	9.2 (6.8)	0. (0.)	0. (0.)
	CLEF	0.003	69.9 (4.4)	16.2 (8.1)	0.5 (0.6)	2.3 (2.2)
		0.005	75.0 (3.6)	8.1 (6.4)	0.2 (0.5)	0.9 (1.2)
	DAMEX	0.003	0.6 (0.2)	1.7 (1.4)	32.9 (5.6)	45.4 (5.9)
		0.005	0.1 (0.4)	2.4 (1.5)	18.3 (5.5)	59.1 (5.9)
	n = 1e5
	CLEF-asymptotic 0.003	73.2 (3.7)	9.5 (6.7)	0.9 (0.8)	4.7 (2.7)
		0.005	72.6 (4.4)	11.7 (7.6)	0.1 (0.4)	0.5 (0.9)
	CLEF-Peng	0.003 79.9 (0.2)	0.2 (1.0)	0. (0.)	0.1 (0.4)
		0.005 80.0 (0.)	0. (0.)	0. (0.)	0. (0.)
	CLEF-Hill	0.003	77.0 (2.0)	6.1 (4.6)	0. (0.)	0. (0.)
		0.005	67.2 (4.8)	22.8 (10.4)	0. (0.)	0. (0.)
	CLEF	0.003	75.2 (3.2)	7.5 (5.9)	0.0 (0.2)	0.2 (0.5)
		0.005	77.9 (2.3)	3.2 (3.9)	0.02 (0.1)	0.02 (0.1)
	DAMEX	0.003	0.04 (0.2)	1.3 (1.0)	24.4 (6.7)	54.2 (7.0)
		0.005	0.1 (0.3)	1.9 (1.6)	10.3 (3.7)	67.6 (4.7)

Table 4 .

 4 2: Same setting as Table4.1 with δ = 0.0001

	n = 5e4	k/n	recovered subset errors superset errors other errors
	CLEF-asymptotic 0.003 71.8 (2.4)	2.3 (2.5)	7.8 (2.8)	41.9 (19.3)
		0.005 73.5 (2.8)	3.7 (3.8)	4.8 (2.5)	25.8 (12.2)
	CLEF-Peng	0.003 79.7 (0.7)	1.0 (2.5)	0. (0.)	3.9 (2.7)
		0.005 80.0 (0.1)	0.1 (0.4)	0. (0.)	0.9 (0.9)
	CLEF-Hill	0.003 79.5 (0.8)	0.3 (1.1)	0.5 (0.8)	142.2 (33.2)
		0.005 79.2 (1.0)	1.6 (2.3)	0. (0.)	0.2 (0.5)
	CLEF	0.003 69.9 (4.4)	16.2 (8.1)	0.5 (0.6)	2.3 (2.2)
		0.005 75.0 (3.6)	8.1 (6.4)	0.2 (0.5)	0.9 (1.2)
	DAMEX	0.003	0.6 (0.2)	1.7 (1.4)	32.9 (5.6)	45.4 (5.9)
		0.005	0.1 (0.4)	2.4 (1.5)	18.3 (5.5)	59.1 (5.9)
	n = 1e5
	CLEF-asymptotic 0.003 75.7 (2.8)	3.7 (3.8)	2.0 (1.4)	11.0 (5.5)
		0.005 76.0 (2.9)	5.6 (4.5)	0.4 (0.7)	1.9 (1.9)
	CLEF-Peng	0.003 79.9 (0.2)	0.2 (1.0)	0. (0.)	0.1 (0.4)
		0.005 80. (0.)	0. (0.)	0. (0.)	0. (0.)
	CLEF-Hill	0.003 79.5 (1.0)	1.2 (2.3)	0. (0.)	0.1 (0.2)
		0.005 75.4 (2.8)	8.7 (5.2)	0. (0.)	0. (0.)
	CLEF	0.003 75.2 (3.2)	7.5 (5.9)	0.0 (0.2)	0.2 (0.5)
		0.005 77.9 (2.3)	3.2 (3.9)	0.02 (0.1)	0.02 (0.1)
	DAMEX	0.003 0.04 (0.2)	1.3 (1.0)	24.4 (6.7)	54.2 (7.0)
		0.005 0.1 (0.3)	1.9 (1.6)	10.3 (3.7)	67.6 (4.7)

  j ∈ {1, . . . , d} \ α k , and with probability 1 d according to component k ∈ {K, . . . , K + d 1 } according to:R i ∼ P areto(1)|{R i > r 0 } W i = 1 ε i,j ∼ 1 + Exp(λ k ), j ∈ {1, . . . , d} \ {k}.

2.2. EXTREME VALUE THEORY

2.4. ESTIMATION OF THE SUPPORT OF THE ANGULAR MEASURE

3.4. EMPIRICAL CRITERION AND IMPLEMENTATION

ACKNOWLEDGMENTS

Part of this work has been funded by the the 'LabEx Mathématiques Hadamard' (LMH) project, by the 'AGREED' project from the PEPS JCJC program (INS2I, CNRS) and by the chair 'Machine Learning for Big Data' from Télécom ParisTech. The authors would like to thank Benjamin Renard for interesting discussions about the hydrological use case and for sharing the data.

Dirichlet model : One natural model for probability distributions on a simplex S α is the Dirichlet family. Such distributions admit a density ϕ α with respect to the (|α| -1) Lebesgue measure which we denote dw for simplicity. It can be parameterized by a mean vector m α ∈ S α and a concentration parameter ν α > 0, so that for w ∈ S α , .14) In this paper we model Φ α by a Dirichlet distribution with unknown parameters m α , ν α .

Using the standard fact that for such a distribution, Sα wϕ α (w|m α , ν α )dw = m α , the moment constraint (5.8) becomes:

(5.15)

where we have set

The Dirichlet mixture model may be summarized as follows:

Model 1 (Dirichlet mixture model).

1. Consider a standardized random vector V such that V j has standard Pareto distribution (see Section 5.2.1)

3. Fix some high radial threshold r 0 , typically a large quantile of the observed radii.

4. Let Z be a hidden variable indicating the mixture component responsible for W, and k such that

The unknown parameters of the model are

, where ν k > 0 and where π, m must satisfy the constraints (5.13) and (5.15)

Figure 5.2 illustrates Model 1 in the bivariate case. Sub-asymptotic model: incorporating a noise : Recall from (5.7) that Φ is the limiting distribution of W for large R's. In practice the observed angles corresponding to radii R > r 0 follow a sub-asymptotic version of Φ. In particular, the marginal variables V j have continuous Pareto distributions, so that P(V j > 1) = 1, j ∈ {1, . . . , d}. As a consequence with probability 1, all the V i = (V i,1 , . . . , V i,d ), 1 ≤ i ≤ n, lie in the central cone C {1,...,d} and the angles W i = (W i,1 , . . . W i,d ) lie in