Anne Sabourin
email: anne.sabourin@telecom-paristech.fr

SEMI-PARAMETRIC MODELING OF EXCESSES ABOVE HIGH MULTIVARIATE THRESHOLDS WITH CENSORED DATA

Keywords: Multivariate extremes, censored data, data augmentation, semi-parametric Bayesian inference, MCMC algorithms

How to include censored data in a statistical analysis is a recurrent issue in statistics. In multivariate extremes, the dependence structure of large observations can be characterized in terms of a non parametric angular measure, while marginal excesses above asymptotically large thresholds have a parametric distribution. In this work, a flexible semi-parametric Dirichlet mixture model for angular measures is adapted to the context of censored data and missing components. One major issue is to take into account censoring intervals overlapping the extremal threshold, without knowing whether the corresponding hidden data is actually extreme. Further, the censored likelihood needed for Bayesian inference has no analytic expression. The first issue is tackled using a Poisson process model for extremes, whereas a data augmentation scheme avoids multivariate integration of the Poisson process intensity over both the censored intervals and the failure region above threshold. The implemented MCMC algorithm allows simultaneous estimation of marginal and dependence parameters, so that all sources of uncertainty other than model bias are captured by posterior credible intervals. The method is illustrated on simulated and real data.

Introduction

Data censoring is a commonly encountered problem in multivariate statistical analysis of extreme values. A 'censored likelihood' approach makes it possible to take into account partially extreme data (non concomitant extremes): coordinates that do not exceed some large fixed threshold are simply considered as left-censored. Thus, the possibly misleading information carried by non-extreme coordinates is ignored, only the fact that they are not extreme is considered [START_REF] Smith | Multivariate threshold methods[END_REF]; [START_REF] Ledford | Statistics for near independence in multivariate extreme values[END_REF]; [START_REF] Smith | Markov chain models for threshold exceedances[END_REF], see also [START_REF] Thibaud | Efficient inference and simulation for elliptical pareto processes[END_REF] or [START_REF] Huser | A comparative study of likelihood estimators for multivariate extremes[END_REF]). However, there are other situations where the original data is incomplete. For example, one popular way to obtain large sample sizes in environmental sciences in general and in hydrology in particular, is to take into account data reconstructed from archives, which results in a certain amount of left-and right-censored data, and missing data. As an example, what originally motivated this work is a hydrological data set consisting of daily water discharge recorded at four neighboring stations in the region of the Gardons, in the south of France. The extent of systematic recent records is short (a few decades) and varies from one station to another, so that standard inference using only 'clean' data is unfeasible (only 3 uncensored multivariate excesses of 'large' thresholdsfixed after preliminary uni-variate analysis-are recorded). Historical information is available, starting from the 17 th century, a large part of it being censored: only major floods are recorded, sometimes as an interval data (e.g. 'the water level exceeded the parapet but the Mr. X's house was spared'). These events are followed by long 'blank' periods during which the previous record was not exceeded. Uni-variate analysis for this data set has been carried on by [START_REF] Neppel | Flood frequency analysis using historical data: accounting for random and systematic errors[END_REF] but a multivariate analysis of extremes has never been accomplished, largely due to the complexity of the data set, with multiple censoring.

While modeling multivariate extremes is a relatively well marked out path when 'exact' (non censored) data are at stake, many fewer options are currently available for the statistician working with censored data. The aim of the present paper is to provide a flexible framework allowing multivariate inference in this context.

Here, the focus is on the methodology and the inferential framework is mainly tested on simulated data with a censoring pattern that resembles that of the real data. A detailed analysis of the hydrological data raises other issues, such as, among others, temporal dependence and added value of the most ancient data. These questions are addressed in a separate paper, intended for the hydrological community [START_REF] Sabourin | Combining regional estimation and historical floods: a multivariate semi-parametric peaks-over-threshold model with censored data[END_REF] 1 .

Under a standard assumption of multivariate regular variation (see Section 2), the distribution of excesses above large thresholds is characterized by parametric marginal distributions and a non-parametric dependence structure that is independent from threshold. Since the family of admissible dependence structures is, by nature, too large to be fully described by any parametric model, non-parametric estimation has received a great deal of attention in the past few years [START_REF] Einmahl | Nonparametric estimation of the spectral measure of an extreme value distribution[END_REF][START_REF] Einmahl | Maximum empirical likelihood estimation of the spectral measure of an extreme-value distribution[END_REF][START_REF] Guillotte | Non-parametric bayesian inference on bivariate extremes[END_REF]. To the best of my knowledge, the non parametric estimators of the so-called angular measure (which characterizes the dependence structure among extremes) are only defined with exact data and their adaptation to censored data is far from straightforward.

For applied purposes, it is common practice to use a parametric dependence model. A widely used one is the Logistic model and its asymmetric and nested extensions [START_REF] Gumbel | Distributions des valeurs extrêmes en plusieurs dimensions[END_REF][START_REF] Coles | Modeling extreme multivariate events[END_REF][START_REF] Stephenson | High-dimensional parametric modelling of multivariate extreme events[END_REF][START_REF] Stephenson | Simulating multivariate extreme value distributions of logistic type[END_REF][START_REF] Fougères | Models for dependent extremes using stable mixtures[END_REF]. In the logistic family, censored versions of the likelihood are readily available, but parameters are subject to non linear constraints and structural modeling choices have to be made a priori, e.g., by allowing only bi-variate or tri-variate dependence between closest neighbors.

One semi-parametric compromise consists in using mixture models, built from a potentially infinite number of parametric components, such as the Dirichlet mixture model (DM), first introduced by [START_REF] Boldi | A mixture model for multivariate extremes[END_REF]. They have shown that it can approach arbitrarily well any valid angular measure for extremes. A re-parametrized version of the DM model [START_REF] Sabourin | Bayesian dirichlet mixture model for multivariate extremes: A re-parametrization[END_REF], allows for consistent Bayesian inference -thus, a straightforward uncertainty assessment using posterior credible sets -with a varying number of mixture components via a reversible-jumps algorithm. The approach is appropriate for data sets of moderate dimension (typically, d ≈ 5).

The purpose of the present work is to adapt the DM model to the case of censored data. The difficulties are two-fold: First, from a modeling perspective, when the censoring intervals overlap the extremal thresholds (determined by preliminary analysis), one cannot tell whether the event must be treated as extreme. The proposed approach here consists in reformulating the Peaks-over-threshold (POT) model originally proposed by [START_REF] Boldi | A mixture model for multivariate extremes[END_REF] and [START_REF] Sabourin | Bayesian dirichlet mixture model for multivariate extremes: A re-parametrization[END_REF], in terms of a Poisson model, in which the censored regions overlapping the threshold have a well-defined likelihood. The second challenge is numerical and algorithmic: for right-censored data above the extremal threshold (not overlapping it), the likelihood expression involves integrals of a density over rectangular regions, which have no analytic expression. The latter issue is tackled within a data augmentation framework, which is implemented as an extension of [START_REF] Sabourin | Bayesian dirichlet mixture model for multivariate extremes: A re-parametrization[END_REF]'s algorithm for Dirichlet mixtures.

An additional issue addressed in this paper concerns the separation between marginal parameters estimation and estimation of the dependence structure. Performing the two steps separately is a widely used approach, but it boils down to neglecting marginal uncertainty, which confuses uncertainty assessment about joint events such as probabilities of failure regions. It also goes against the principle of using regional information together with the dependence structure to improve marginal estimation, which is the main idea of the popular regional frequency analysis in hydrology. In this paper, simultaneous inference of marginal and dependence parameters in the DM model is performed, which amounts in practice to specifying additional steps for the marginal parameters in the MCMC sampler. The rest of this paper is organized as follows: Section 2 recalls the necessary background for extreme values modeling. The main features of the Dirichlet mixture model are sketched. This POT model is then reformulated as a Poisson model, which addresses the issue of variable threshold induced by the fluctuating marginal parameters. Censoring is introduced in Section 3. In this context, the Poisson model has the additional advantage that censored data overlapping threshold have a well defined likelihood.The lack of analytic expression for the latter is addressed by a data augmentation scheme described in Section 4. The method is illustrated by a simulation study in Section 5: marginal performance in the DM model and in an independent one (without dependence structure) are compared, and the predictive performance of the joint model in terms of conditional probabilities of joint excesses is investigated. The model is also fitted to the hydrological data. Section 6 concludes. Most of the technicalities needed for practical implementation, such as computation of conditional distributions, or details concerning the data augmentation scheme and its consistency are relegated to the appendix.

X t = (-1/ log(F 1 (Y 1,t )), . . . , -1/ log(F d (Y d,t ))) .
Then the X j,t 's have unit-Fréchet distribution, P(X j,t ≤ x) = e -1/x , x > 0. It is mathematically convenient to switch to pseudo-polar coordinates,

R = d j=1 X j (radial component), W = 1 R X ∈ S d (angular component) ,
where S d = {x : x j ≥ 0, d j=1 x j = 1} is the unit simplex. The radial variable R corresponds to the 'amplitude' of the data whereas the angular component W characterizes their 'direction'. Asymptotic theory [START_REF] Resnick | Extreme values, regular variation, and point processes, volume 4 of Applied Probability[END_REF][START_REF] Beirlant | Statistics of extremes: Theory and applications[END_REF][START_REF] Coles | An introduction to statistical modeling of extreme values[END_REF] tells us that, under mild assumptions on F (namely, belonging to a multivariate maximum domain of attraction), an appropriate model, commonly referred to as a multivariate Peaks-over-threshold (POT) model, for (R, W) over high radial thresholds r 0 , is

P(R > r, W ∈ A | R > r 0 ) = r r 0 H(A) , r 0 > r, A ⊂ S d , (2.1)
where H is the so-called 'angular probability measure' (called the 'angular measure' in the remainder of this paper). The angular measure is thus the limiting distribution of the angle, given that the radius is large. Concentration of H's mass in the middle of the simplex indicates strong dependence at extreme levels, whereas Dirac masses only on the vertices characterizes asymptotic independence. This paper focuses on the case where H is concentrated on the interior of the simplex, so that all the variables are asymptotically dependent.

Because of the standardization to unit Fréchet, a probability measure H on S d is a valid angular measure if and only if

S d w j dH(w) = 1 d (1 ≤ j ≤ d)
. This moments constraint is the only condition on H, so that the angular measure has no reason to be part of any particular parametric family.

2.2. Dirichlet mixture angular measures. In this paper, the angular measure H is modeled by a Dirichlet mixture distribution [START_REF] Boldi | A mixture model for multivariate extremes[END_REF][START_REF] Sabourin | Bayesian dirichlet mixture model for multivariate extremes: A re-parametrization[END_REF]. A Dirichlet distribution can be characterized by a shape ν ∈ R + and a center of mass µ ∈ S d , so that its density with respect to the d -

1 dimensional Lebesgue measure dw = dw 1 • • • dw d-1 , is diri ν,µ (w) = Γ(ν) d j=1 Γ(νµ j ) d j=1 w νµ j -1 j (w ∈ S d ).
(2.2)

A parameter for a k-mixture is of the form

ψ = ((p 1 , . . . , p k ), (µ 1 , . . . , µ k ), (ν 1 , . . . , ν k )) ,
with weights p m > 0, such that k m=1 p m = 1. This is summarized by writing ψ = (p 1:k , µ 1:k , ν 1:k ). The corresponding mixture density is

h ψ (w) = k m=1 p m diri ν,µ m (w) .
(2.3) the moments constraint is satisfied if and only if

k m=1 p m µ m = (1/d, . . . , 1/d) ,
which, in geometric terms, means that the center of mass of the µ 1:k 's, with weights p 1:k , must lie at the center of the simplex. As established by [START_REF] Boldi | A mixture model for multivariate extremes[END_REF] and mentioned in the introduction, the family of Dirichlet mixture densities satisfying the moments constraint is weakly dense in the space of admissible angular measure. In addition, in a Bayesian framework, [START_REF] Sabourin | Bayesian dirichlet mixture model for multivariate extremes: A re-parametrization[END_REF] have

shown that the posterior is weakly consistent under mild conditions. These two features put together make the Dirichlet mixture model an adequate candidate for modeling the angular components of extremes.

2.3. Model for margins. The above model for excesses concerns standardized versions X t of the data Y t involving marginal cumulative distribution function

F j (1 ≤ j ≤ d)
, which have to be estimated. As a consequence of uni-variate extreme value theory [START_REF] Pickands | Statistical inference using extreme order statistics[END_REF], uni-variate excesses above large thresholds v j (1 ≤ j ≤ d) are approximately distributed according to a Generalized Pareto distribution with parameters ξ j (shape) and σ j (scale parameter),

P (Y j > y | Y j > v j ) ≈ v j →∞ (1 + ξ j y -v j σ j ) -1/ξ j .
A widely used method to model the largest excesses is a follows: Define a high multivariate threshold v = (v 1 , . . . , v d ) and call 'marginal excess' any Y j,t > v j . Then, marginal excesses above v j are modeled as generalized Pareto random variables with parameters ξ j and σ j . The marginal parameters are gathered into a (2d)-dimensional vector χ = (log(σ 1 ), . . . , log(σ d ), ξ 1 , . . . , ξ d ) ∈ R 2d .

Let F v j denote the j th marginal distribution conditionally on Y j not exceeding v j , and let ζ j = P(Y j > v j ) denote the probability of excursion above v j . The j th marginal model (1

≤ j ≤ d) is thus F (χ) j (y) = P(Y j,t ≤ y | ξ j , σ j ) =    1 -ζ j 1 + ξ j y-v j σ j -1/ξ j (y ≥ v j ), (1 -ζ j )F v j (y) (y < v j ).
(2.4)

It is common practice [START_REF] Coles | Modeling extreme multivariate events[END_REF][START_REF] Davison | Models for exceedances over high thresholds[END_REF] to use an empirical estimate ζ = ( ζ1 , . . . , ζd ) for the vector of probabilities of marginal excursion, and to ignore any estimation error, so that ζ is identified to ζ is the sequel.

2.4. Joint inference in a Poisson model. When it comes to simultaneous estimation of the margins and of the angular measure, the angular model (2.1) for radial excesses is difficult to handle, because a radial failure region r > r 0 on the Fréchet scale (i.e. , in terms of X's) corresponds to a complicated shaped failure region on the original scale, which depends on the marginal parameters and, accordingly, potentially contains a varying number of data points. It seems more reasonable to use a failure region which is fixed on the original scale (in terms of Y's). Further, a common criticism towards radial failure regions [START_REF] Ledford | Statistics for near independence in multivariate extreme values[END_REF] is that the marginal Pareto model is not valid near the axes of the positive orthant. Last but not least, censoring occurs along the directions of the Cartesian coordinate system, which prevents using the polar model (2.1) as it is.

To address these issues, the statistical model for threshold excesses developed in this paper uses a 'rectangular' threshold. Also, it will be very convenient (see Section 3.3) to adopt a Poisson process representation of extremes (see e.g. [START_REF] Coles | Modeling extreme multivariate events[END_REF] as an alternative to the POT model (2.1), with a 'censored likelihood' near the axes.

Poisson model. Under the same condition of domain of attraction as above, the point process formed by time-marked, standardized and suitable re-scaled data converges in distribution to a Poisson process (see e.g. [START_REF] Resnick | Extreme values, regular variation, and point processes, volume 4 of Applied Probability[END_REF][START_REF] Resnick | Heavy-Tail Phenomena: Probabilistic and Statistical Modeling[END_REF][START_REF] Coles | Modeling extreme multivariate events[END_REF],

n t=1 1 ( t n , X t n ) w -→ PP( ⊗ λ) , (2.5)
in the space of point measures on

([0, 1]×E), where E = [0, ∞] d \{0}.
The temporal component of the limiting intensity measure denotes the Lebesgue measure on R and λ, the so-called exponent measure, is homogeneous of order -1, and is related to the angular measure H via dλ(r, w) = d r 2 dr dH(w) .

(2.6)

From a statistical perspective, consider a failure region

A v = E \ [0, v],
where v is the high multivariate threshold introduced in section 2.3 and

[0, v] = [0, v 1 ] × • • • × [0, v d ] \ {0}. Call 'excess above v' any point Y t in A v , as opposed to marginal excesses Y j,t > v j . The Fréchet re-scaled multivariate threshold is u = T(v) = -1/ log(1 -ζ)
and does not depend on χ. Consider the re-scaled region on the Fréchet scale

A u,n = 1 n T(A v ) = [0, ∞] d \ [0, u 1 n ] × • • • × [0, u d n ] ,
and denote A u = A u,1 . Applying the marginal transformations

T χ j (y) = -1/ log F (χ)
j (y) , the marginal variables X j,t = T χ j (Y j,t ) have unit Fréchet distribution, as required in (2.5). The point process N = n t=1 1 ( t n , X t n ) composed of the excesses X t ∈ A u (i.e. Y t ∈ A u ) is modeled according to the limit in (2.5),

n t=1 1 ( t n , X t n ) ∼ PP( ⊗ λ) on [0, 1] × A u,n ,
where λ is of the form (2.6), with angular component H written as a Dirichlet mixture of the form (2.3).

Joint likelihood of uncensored data. Let θ = (χ, ψ) be the parameter for the joint model. As explained at the beginning of this section, the model likelihood needs to be expressed in Cartesian coordinates. The density of an exponent measure λ with respect to the d-dimensional Lebesgue measure dx = dx 1 • • • dx d , is [START_REF] Coles | Modeling extreme multivariate events[END_REF], Theorem 1)

dλ dx (x) = d . r -(d+1) h(w) .
(2.7)

Denote by λ ψ the exponent measure corresponding to the Dirichlet mixture h ψ . Then, in the simplified case where the Y j,t 's are exactly observed and where the marginals F j 's below threshold are known, the likelihood in the Poisson model over

A v is L v ({y t } 1≤t≤n , θ) ∝ e -n λ ψ (Au) nv i=1 dλ ψ dx (x t i ) j:y j,t i >v j J χ j (y j,t i ) , (2.8)
where t 1 , . . . , t nv are the occurrence times of excesses y t i ∈ A v , x j,t i = T χ j (y j,t i ), and the J χ j are Jacobian terms resulting from marginal transformations T χ j (see Appendix A for details).

Censored model

3.1. Causes of censoring. The presence of censored observations is the result of two distinct causes: first, data are partially observed, which results in interval-or right-censoring, which we call natural censoring. In addition, observed data points that exceed at least one threshold in one direction do not necessarily exceed all thresholds, so that the marginal extreme value model does not apply. Following [START_REF] Ledford | Statistics for near independence in multivariate extreme values[END_REF], those components are also considered as left-censored.

This second censoring process is thus a consequence of an inferential framework which is designed for analyzing extreme values only, and we call it inferential censoring.

The total censoring process C , which results from the juxtaposition of natural and inferential censoring, is assumed to be non informative. This means that (see also [START_REF] Gómez | Frequentist and bayesian approaches for interval-censored data[END_REF], if F is the marginal c.d.f. for Y j and f is the marginal density, then Y j 's distribution conditional upon having observed only the left and right censoring bounds

(L, R) is f ( • )/[F (R) -F (L)].
This definition is easily extended to the multivariate case by replacing F (R) -F (L) by the integral of the density over the censored directions (see e.g. [START_REF] Schnedler | Likelihood estimation for censored random vectors[END_REF], for a proof of consistency of maximum censored likelihood estimators).

3.2. Natural censoring. Call 'Natural censoring' the one which occurs independently from the choice of an extreme threshold v by the statistician. The observed process is denoted O = (O t ) t , with O t = (O 1,t , . . . , O d,t ). One marginal observation O j,t consists in a label κ j,t indicating presence or absence of censoring, together with the exact data Y j,t (if observed) or the censoring bounds (L j,t , R j,t ) (where L j,t may be set to 0 in the case of left-censoring or missing data and L j,t = +∞ in the case of right-censored or missing data). In the sequel, κ j,t = 0 (resp. 1, 2, 3) refer respectively to missing, exact, right-and left-censored data.

In this context, the 'position' of a marginal data O j,t with respect to the threshold is not necessarily well defined for censored data. Recall that we consider a case where the process of interest (Y t ) t≥0 is stationary, whereas the censoring bounds vary with time, as a result of external factors on the observation process. When the censoring interval overlaps the threshold, i.e. , v j ∈ (L j,t , R j,t ), κ j ∈ {2, 3}), the statistician does not know if an excess occurred or not. This situation is described here as O j,t marginally overlapping the threshold. The different positions of O j,t with respect to v j encountered in the data set of interest in this paper are summarized in Figure 1. For a multivariate observation O t , if at least one coordinate marginally overlaps the threshold or is missing, and if the others are below the threshold, then the position of O t with respect to the multivariate threshold v is undetermined. Indeed, O t is below threshold if all its marginals are below the corresponding marginal threshold, and above threshold (in the failure region) otherwise. In the undetermined case, O t is qualified as globally overlapping the threshold.

Inferential censoring below threshold. Since the marginal distributions

F v j 's, conditional upon not exceeding v j , are unknown, the X j,t 's such that Y j,t < v j Above Below Undetermined (1) (3) (2) Type Position (1) (3) (3) (2) v Discharge Figure 1.
Position of marginal data points with respect to a marginal threshold v (horizontal line). Black dots: marginal data points of type κ = 1; vertical arrows: data of type κ ∈ {2, 3}.

are not available. Instead of attempting to estimate the F v j 's, one option is to censor the Fréchet-transformed components below threshold. More precisely, for a raw observation O j,t = (κ j,t , Y j,t , L j,t , R j,t ), let us denote by C χ j,t = (κ j,t , X j,t , Lj,t , Rj,t ) the corresponding 'Fréchet transformed' and censored one, and

C χ t the multivari- ate observation (C χ 1,t , . . . , C χ d,t ). The transformation O t → C χ t is illustrated in
Figure 2 in the bi-variate case.

In a nutshell, inferential censoring occurs when the censoring intervals are below threshold or marginally overlapping it, or when an exact data component Y j,t < v j is recorded. A formal definition of C χ j is as follows:

• (κ j,t , X j,t ) =      (3, NA) if κ j,t = 1 and Y j,t < v j , (0, NA) if κ j,t = 2 and L j,t < v j , κ j,t , T χ j (Y j,t ) otherwise. • Lj,t = 0 if L j,t < v j , T χ j (L j,t ) otherwise. • Rj,t = u j if R j,t < v j , T χ j (R j,t ) otherwise.
In the above definition, NA stands for a missing value and it is understood that T χ j (NA) = NA.

In the end, observations globally overlapping threshold have their marginal lower bounds Lj,t set to zero if L j,t < v j . The interest of using a Poisson model becomes clear at this point. Indeed, censored observations C χ t obtained from observations O t globally overlapping threshold correspond to events of the kind 'The observation at time t belongs to [0, Rt ]', which, by contrapositive, means 'No point is observed outside of [0, Rt ] between t and t + 1'. This is written in terms of the Poisson process N as

N t n , t + 1 n × 0, Rt n c = 0 , (3.1)
which is a measurable event with respect to N. The overlapping observations thus have a well defined likelihood in a Poisson model, as detailed in the next section, whereas they could not be taken into account in Sabourin and Naveau ( 2014)'s POT model.

3.4.

Poisson likelihood with censored and missing data. Due to the combination of natural and inferential censoring, the data set (from which missing days are excluded) is decomposed into data in the failure region, data overlapping threshold and data below threshold. Let n v , n v and n v be the respective number of observations in each category. The number of non missing days is thus

n obs = n v + n v + n v
, and the number of 'determined' data (i.e. not overlapping v) is

n det = n v + n v . The n v Fréchet-transformed observations {C χ t i }(i ∈ {1, . . . , n v }) correspond to events of the kind X t i ∈ [ Lt i , Rt i ] (rectangular region)
where Rj,t i = Lj,t i = X j,t i in the case of exact data.

Observations overlapping threshold correspond to events of the kind (3.1) introduced in the previous section. When a limited number I of right censoring bounds R t = (R 1,t , . . . , R d,t ) are present, it is convenient to classify these overlapping events accordingly, writing n v = I i=1 n i where n i is the number of observations with right censoring bound R t i . Since the Poisson process is temporally exchangeable, there is no loss of generality in assuming that the latter observations occur at consecutive dates (t i , . . . , t i + n i -1). The region E \ [0, Rt i ] is 'missed' by the Fréchet re-scaled process X t during this time period. With theses notations, the censored likelihood in the Poisson model may be written

L v (O, θ) = exp -n det λ ψ (A u ) - I i=1 n i λ ψ (A i ) × • • • • • • nv i=1 [ Lt i , Rt i ] dλ ψ dx d i (x) j:y j,t i >v j J χ j (y j,t i ) , (3.2)
where notation ' d i (x)' in the integral terms is a shorthand for 'the Lebesgue measure of dimension equal to that of [ Lt i , Rt i ]' when the latter is greater than one, or 'the Dirac mass at x t i = Lt i = Rt i ' for exact data. Compared with the uncensored likelihood (2.8), n has been replaced with n obs , the exponential term for the non overlapping data follows from

exp - n det n obs λ ψ A u,n obs = exp -n det λ ψ A u ,
and a similar argument yields the additional terms exp(-n i λ ψ (A i )) for overlapping data.

At this stage, the model has been entirely specified. The remaining issue concerns the treatment of the integral terms

[ Lt i , Rt i ] dλ ψ dx d i (x) (3.3)
and the exponential terms

exp -n det λ ψ (A u ) and exp -n i λ ψ (A i ) , (3.4)
which have no analytic expression, as they require integrating λ ψ over rectangular regions. First, the dimension of numerical integration can be reduced as far as 'missing coordinates' are involved, because partial integration of λ ψ over [0, ∞] in one direction has an exact expression. The model is stable under marginalization, in the sense that the obtained marginal exponent measures correspond again to Dirichlet mixtures on a lower dimensional simplex (see Appendix B for details). However, no closed form is available for the integral in the remaining censored directions, nor for the exponent measures of A u or the A i 's. This problem is tackled in the next section via a data augmentation method.

Data augmentation

4.1. Background. In a Bayesian context, one major objective is to generate parameter samples approximately distributed according to the posterior. In classical MCMC algorithms, the value of the likelihood is needed to define the transition kernel. Evaluating the integrated likelihood L v (O, θ) at each iteration of the algorithm seems unmanageable: The dimension of integration can grow up to d, for each observation, and the shape of the integrand varies from one iteration to another, which is not favorable to standard quadrature methods. In particular, large or low values of the shape parameters ν m in ψ induce concentration of the integrand around the centers µ m or unboundedness at the simplex boundaries. Instead, data augmentation methods (see e.g. [START_REF] Tanner | The calculation of posterior distributions by data augmentation[END_REF]Van Dyk and Meng, 2001) treat missing or partially observed data as additional parameters, so that the numerical integration step is traded against an increased dimension of the parameter space. In this section, [ • ] denotes the distribution of a random quantity as well as its density with respect to some appropriate reference measure. Proportionality between σ-finite measures is denoted by ∝. Thus, [θ] is the prior density and

[θ | O] ∝ [θ] L v (O, θ) is the posterior.
The main idea is to define an augmentation space Z, and a probability measure [ • |O] + , on the augmented space Θ × Z, conditional on the observations, which may be sampled using classical MCMC methods, and which is consistent with the 'objective distribution' on Θ. That is, the posterior on Θ must be obtained by marginalization of

[ • | O] + , [θ|O] = Z [z, θ | O] + dz . (4.1)
In the sequel, the invariant measure (or its density) [z, θ | O ] + is referred to as the augmented posterior. can be directly simulated as block proposals in a Metropolis-within-Gibbs algorithm [START_REF] Tierney | Markov chains for exploring posterior distributions[END_REF]. To ensure the marginalization condition (4.1), the augmented posterior

[Z, θ| O] + is not the same as [Z, θ| O]. Instead it has a density of the form [z, θ | O] + ∝ [z, O | θ ] [ θ ] ϕ(z) , (4.
2) where ϕ(z) is any weight function allowing to enforce the consistency condition (4.1). The [θ] terms cancel out and the latter condition is equivalent to

L v (O, θ) ∝ [z, O | θ] ϕ(z) dz . (4.3)
In the end, a posterior sample from [θ|O] is simply obtained by ignoring the Zcomponents from the one produced with the 'augmented' Markov chain.

In our case, the augmented data Z consist of two parts, Z = (Z above , Z ). The first one, Z above = {Z t i } i≤nv , accounts for integral terms (3.3), while the second one, Z = (Z u , {Z i } i≤I ), accounts for the exponential terms (3.4). The Z t i 's have an intuitive interpretation, which is standard in data augmenting methods: they replace the censored X j,t i 's. On the contrary, the Z i 's and the Z u are just a computational trick accounting for the exponential terms: they are the points of independent Poisson processes defined on 'nice' radial sets containing the failure regions of interest A i 's and A u , and ϕ is a smoothed version of an indicator function of the failure regions.

The augmented likelihood factorizes as

[ z, O | θ ] = nv i=1 [ z t i , O t i | θ ] [ z u | ψ] I i=1 [z i | ψ] , (4.4)
and the weight function ϕ is of the form

ϕ(z) = ϕ u (z u ) I i=1 ϕ i (z i ) .
A , where ψ is the DM parameter and W is an angular data set consisting of the angular components W t of data X t that have been normalized to Fréchet margins in a preliminary step.

In contrast, the present algorithm handles 'raw' data (not standardized to Fréchet margins), with censoring of various types as described in Section 3, so as to approach the full posterior distribution (marginal and dependence parameters) [ θ | O ] . In practice, this amounts to allowing two additional move types in the Metropolis-within-Gibbs sampler: marginal moves (updating the marginal parameters) and augmentation moves updating the augmentation data Z described above. The reversible jumps and the moves updating the DM parameters keeping the dimension constant are unchanged.

Simulations and real case example

Keeping in mind the application, the aim of this section is to verify that the algorithm provides reasonable estimates with data sets that 'resemble' the particular one motivating this work.

After a description of the experimental setting, an example of results obtained with a single random data set is given, then a more systematic study is conducted over 50 independent data sets. Finally, a brief description of the results obtained with the original hydrological data is provided. The latter part is kept short because, as mentioned in introduction, from a hydrological point of view, a full discussion of the benefits brought by historical information, as well as treatment of temporal dependence in the case of heavily censored data is needed and will be the subject of a separate paper.

5.1. Experimental setting. For this simulation study, the marginal shape parameters are constrained to be equal to each other, ξ 1 = • • • = ξ d := ξ, in accordance with the regional frequency analysis hydrological framework [START_REF] Hosking | Regional frequency analysis: an approach based on L-moments[END_REF], where the different gauging stations under study are relatively close to each other (in the same watershed). Also, the dimension is set to d = 4, as it is for the hydrological data set of interest.

Preliminary likelihood maximization (with respect to the marginal parameters) is performed on the hydrological data, again assuming independence between locations for the sake of simplicity and imposing a common shape parameter (this latter hypothesis not being rejected by a likelihood ratio test). Then, data sets are simulated according to the model for excesses (2.1), with marginal parameters and threshold excess probabilities (for daily data) approximately equal to the inferred ones, i.e.

ζ ≈ (0.021, . . . , 0.021) , ξ = 0.4 , log(σ) = (4.8, 4.6, 5.9, 5.1) , for a total number of days n = 148 401, which is the total number of days in the original data. The four-variate dependence structure is chosen as a Dirichlet mixture distribution h ψ , where A radial threshold for simulation on the Fréchet scale and the number of points simulated above the latter are respectively set to r s = -1/ log(1 -max(ζ)), n rad.exc = n * 4/r s . The remaining n -n rad.exc points are arbitrarily scattered below the radial threshold, so that the proportion of radial excesses is 4/r s , the exponent measure of the region {x ∈ R 4 : x 1 > r s }.

Afterwards, the data set is censored following the real data's censoring pattern, i.e. censoring occurs on the same days and on the same locations (here, a location is a coordinate j ∈ {1, . . . , 4}) as for the real data, with same censoring bounds (on the Fréchet scale), so that the data are observed only if the censoring threshold is exceeded. Finally, in order to account for the loss of information resulting from time dependence within the real data (whereas the simulated data are time independent), only n eff data out of n are kept for inference, where n eff = n/mean cluster size = 118 911 (see Section 6 for an explanation about clusters). The vast majority of data points (real as well as simulated) are either censored or below the multivariate threshold: in the real data set, only 125 multivariate excesses above threshold are recorded, among which only 3 have all their coordinates of type 1 (exact data). In such a context, a simplified inferential framework in which censored data would be discarded is not an option: only 3 points would be available for inference. For simulated data, the threshold is arbitrarily set to the same value as the one defined for real data, i.e. v = (300, 320, 520, 380).

The MCMC sampler described in Section 4.3 and Appendix D is run for each simulated data set, yielding a parameter sample, which approaches the posterior distribution given the censored data. To allow comparison with the default spaceindependent framework in terms of marginal estimates, such as probabilities of a marginal excess, Bayesian inference is also performed in the independent model defined as follows: the marginal models are the same as in (2.4), the Y j,t , 1 ≤ j ≤ 4 are assumed to be independent while the shape parameters ξ j , 1 ≤ j ≤ 4 are, again, equal to each other. A MCMC sampler is straightforwardly implemented for this independent model, following the same pattern as defined in the marginal moves for the full model (see Appendix D). In the sequel, the full Poisson model with Dirichlet mixture dependence structure is referred to as the DM Poisson model, or the dependent model, as opposed to the independent model. 5.2. Illustration of the method with one simulated data set. The censoring pattern described above yields a censored data set which resembles the real data set is term of average number of exact (uncensored) coordinates in each observation, as shown in Figure 3: most of the extracted data have only one exact coordinate. In addition to data above threshold, the number of threshold-overlapping blocks (made of data which position with respect to the threshold is undetermined, see Section 3) is I = 39 in this simulated data set, with block sizes varying between 1 and 39 845, and a total number of 'undetermined' days n v = 112 676. To fix ideas, the bi-variate projections onto the six pairs (i < j) of a simulated data set are displayed in Figure 4.

Let us turn to results obtained with this particular data set. A 'flat' prior was specified for the DM Poisson model parameters and MCMC proposals for the DM parameter were set in a similar way as in [START_REF] Sabourin | Bayesian dirichlet mixture model for multivariate extremes: A re-parametrization[END_REF], which resulted in satisfactory convergence diagnostics after 10 6 iterations, see Appendix D.6 for details. Figure 5 displays the posterior predictive for bi-variate marginalization's of the angular density, obtained via equation (B.2), together with the true density and posterior credible sets around the estimates. The estimated density captures reasonably well the features of the true one and the posterior quantiles are rather concentrated around the true density. This is all the more satisfying 0 1000 2000 3000 4000 S1 S2 q q q q q q q q q q q q 0 1000 2000 3000 4000 0 1000 2000 3000 4000 5000 6000 S1 S3 q q q q q q q q q q q q 0 1000 2000 3000 4000 0 1000 2000 3000 4000 S1 S4 q q q q q q q q q q q 0 500 1000 1500 2000 2500 S2 S3 q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0 500 1000 1500 2000 2500 0 1000 2000 3000 4000 S2 S4 q q q q q q q q 0 1000 2000 3000 4000 5000 6000 0 1000 2000 3000 4000 S3 S4 q q q q q q q q q q q q q q q q q q q q q Figure 4. Bi-variate projection of the simulated data set after censoring. White points correspond to pairs for which both coordinates are observed. Superimposed Gray rectangles (resp. segments) represent pairs for which both coordinates (resp. one coordinate) are (is) censored. The white striped rectangle is the region below the multivariate threshold v.

that, at first view (Figure 4), the censored data used for inference seem to convey little information about the distribution of the angular components. In risk analysis, especially in hydrology, return level plots (i.e. quantile plots) are used to summarize the marginal behavior of extremes. The return level Q for a return period T at location j, when marginal data are distributed according to F χ j and where there is no temporal dependence, is usually defined as the 1 -1/Tquantile of F χ j . Figure 6 compares the return levels obtained both in the dependent and the independent models, together with the true return levels. The posterior estimates in the dependent model are very close to the truth, relatively to to the size of the credible intervals. In contrast, estimation in the independent model under-estimates the return levels, and the true curve lies outside the posterior quantiles at two locations out of four. This single example is however not enough to conclude that the dependence structure improves significantly marginal estimation.

The absence of significant improvement (or deterioration) of marginal estimates is indeed one of the conclusions of the next subsection. q q q q q q q q q q q q q q q q q q q q q q q q q q qqqqqqq q q q q q q q S1 0 1 2 3 0 1250 2500 3750 5000 return period (years, log-scale) discharge :

q q q q q q q q q q q q q q q q q q qqqqq qqq q q q q q q S2 0 1 2 3 0 3750 7500 11250 15000 return period (years, log-scale) discharge :

q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqq qq qq q q q q q q q q S3 0 1 2 3 0 2000 4000 6000 8000 return period (years, log-scale) discharge :

q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q qqq qqq q q q q q q q S4 Figure 6. Return level plots: Quantile versus logarithm (base 10) of the return period at the four locations. Grey points: empirical return level of observed threshold excesses (simulated data); Black line: true curves; Gray dashed line and shaded area: posterior mean and 0.05 -0.95-quantiles in the dependent model with Dirichlet mixture angular structure; Black dash-dotted line and area: idem in the independent model.

5.3.

Simulation study with 50 data sets. The aim of this section is to verify that the posterior distribution of the DM Poisson model parameters is reasonably informative, even with censored data. The procedure described above is applied to generate independently 50 data sets.

Marginal predictive performance. The model's ability to estimate of the probability of a marginal excess is first investigated. Large thresholds (V 1 , . . . , V 4 ) are specified so that their true marginal probability of exceedance is P 0 = 1/(10 * 365), an approximate ten years return level. The quantities of interest are the posterior distributions of ∆ j (θ) = P θ (Y j > V j ), which we hope to be concentrated around P 0 . Each posterior samples (θ ι ) ι∈1..N issued by a MCMC algorithm is transformed into a series of exceedance probabilities, (δ j (θ ι )) ι∈1..N which empirical distribution

Fj = 1 N N ι=1 1 δ j (θι) ( • )
approximates the posterior distribution of ∆ j . The performance of the posterior may then be investigated in terms of posterior quadratic loss,

QL( Fj ) = E Fj (∆ j ) -P 0 2 + Var Fj (∆ j ) .
This loss corresponds to the the predictive model choice criterion (PMCC) [START_REF] Laud | Predictive model selection[END_REF]. In the framework of scoring rules [START_REF] Gneiting | Strictly proper scoring rules, prediction, and estimation[END_REF] the PMCC is not a 'proper score' in general, but it is so when the true probability distribution of the quantity of interest ∆ j is a Dirac mass, which is the case here (Dirac mass at P 0 ). After 1.10 6 MCMC iterations, at least one chain (out of six chains run in parallel for each simulated data set) passed the Heidelberger and Welch's stationarity tests [START_REF] Heidelberger | Simulation run length control in the presence of an initial transient[END_REF] at level 10 -4 , for each data set. Table 1 gathers, for the margins j ∈ {1, . . . , 4}, the mean and standard deviation of the QL scores, normalized by the (squared) true probability P 2 0 for readability,

QL j = QL( Fj ) P 2 0 .
QL j 1 2 3 4 mean 0.55 0.17 0.22 0.32 standard error 0.64 0.14 0.18 0.30 first quartile 0.14 0.08 0.08 0.10 third quartile 0.77 0.20 0.32 0.45

Table 1. Normalized scores QL j for the marginal probability of an excess: mean and standard deviation, first and third quartiles over the 50 data sets. Column j corresponds to an excess at location j (1 ≤ j ≤ 4).

Although the variability of the scores (standard deviations) is relatively high, normalized third quartile less than one indicate that the posterior distribution concentrates in regions where the probability of a marginal excess is of the same order of magnitude as the true probability.

In order to verify that introducing a rather complex dependence structure model does not deteriorate the marginal estimates, the same quadratic loss score is computed with samples issued from the independent model. In view of Figure 7, displaying box-plots of the scores computed in both models (dependent versus independent), there is no significant difference between the two models in terms of marginal estimation. q q q q q q q q q dependent independent 0 2 4 6 8 q q q q q q q q q dependent independent 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 q q q q q q dependent independent 0.0 0.5 1.0 1.5 2.0 2.5 q q q q q q q q q dependent independent Estimation of the angular measure. To assess the performance of the DM Poisson model in terms of estimation of the dependence structure of extremes, a similar scoring procedure as above is followed, with quantities of interest ∆ j defined as the probability of a joint excess of the V i 's, given a marginal excess of V j ,

∆j = P(Y 1 ≥ V 1 , . . . Y 4 ≥ V 4 |Y j ≥ V j ) .
These quantities do not have an explicit expression in the DM model, but are easily approached by standard Monte-Carlo sampling. Lacking a reference model in this context (the independent one is obviously unable to predict these quantities), only the scores in the dependent model are available. Table 2 summarizes the results in terms of normalized scores

QL j = QL( ∆j ) P 2 j
, where P j is the true conditional probability of a joint excess.

A comforting fact is that the dependence structure seems to be even better estimated than the marginal distributions of extremes, in view of Tables 1 and2. 5.4. Real data. Analyzing the hydrological data set presented in the introduction requires an additional declustering step: a multivariate run-declustering scheme was implemented. A cluster starts when one component exeeds the threshold. It and ends when, during τ consecutive days, all components are below there respective threshold, or (when censoring is present) with undertermined position, so that the observer can not ascertain that an excess occurred. The lag parameter τ = 3 was set by considering stability regions of the estimates, and additional physical characteristics of the hydrological catchment, see [START_REF] Sabourin | Combining regional estimation and historical floods: a multivariate semi-parametric peaks-over-threshold model with censored data[END_REF] for details. After-wise the model is fitted to the extracted four-variate QL j 1 2 3 4 mean 0.06 0.18 0.25 0.05 standard error 0.05 0.11 0.13 0.05 first quartile 0.02 0.11 0.17 0.01 third quartile 0.10 0.23 0.33 0.06

Table 2. Normalized QL scores QL j for conditional probabilities of a joint excess: mean, standard deviation, first and third quartiles over the 50 simulated datasets. Column j corresponds to conditioning upon an excess at location j.

cluster maxima. Again, 6 parallel MCMC's of length 10 6 are run, with satisfactory convergence diagnostic after a 2 10 5 burn-in period. The posterior mean estimates of marginal parameters are close to the maximum likelihood estimates mentioned in subsection 5.1. Figure 8 displays bi-variate versions of the four-variate posterior predictive angular measure, together with point-wise 90% posterior credibility sets. The inferred dependence structure is rather complex, which fosters the use of such a flexible semi-parametric model. The posterior credible bounds of the angular density are relatively narrow around the posterior predictive in the central regions of the simplex (which is the segment [0, 1] for bi-variate data) and indicate that asymptotic dependence is present. On the other hand, high levels (even unboundedness) of the predictive density near certain edges indicates a 'weakly asymptotically dependent regime', for the considered pair, i.e. a regime where one component may be large while the other one is not. Possibly, some 'true' asymptotic angular mass is concentrated on these edges, which translates in the Dirichlet mixture model (which only allows mass on the topological interior of the simplex) into high densities near the edges. The observed widening of posterior credible regions near the edges where the density is unbounded is not surprising : an unbounded density corresponds to Dirichlet parameters with components νµ i < 1, for which a small variation of the parameter value induces a large variation of the density near the edges.

Conclusion

In this work, a flexible semi-parametric Bayesian inferential scheme is implemented to estimate the joint distribution of excesses above multivariate high thresholds, when the data are censored. A simulation example is designed on the same pattern as a real case borrowed from hydrology. Although the tuning of the MCMC algorithm requires some care, taking into account all kinds of observations for various censoring bounds allows to obtain satisfactory estimates, despite the loss of information relative to the angular structure induced by the censoring process. In particular, accurate enough estimations of quantities of interest such as marginal or conditional probabilities of an excess of a large threshold can be obtained. The main methodological novelty consists in taking advantage of the conditioning and marginalizing properties of the Dirichlet distributions, in order to simulate augmentation data which 'replace' the missing ones. Also, exponential terms in the likelihood with no explicit expressions are handled by sampling well chosen functionals of augmentation Poisson processes. This new inferential framework opens the road to statistical analysis of the extremes of data sets that would otherwise have been deemed unworkable.

k m=1      p m Γ(ν m ) d j=1 Γ(ν m µ j,m ) d j=1 x νmµ j,m -1 j   d j=1 x j   -(νm+1)      . (A.1)
The likelihood expression (2.8) is simply that of a Poisson process on the region [0, 1] × A u,n . Recall that, for a Poisson process with intensity µ on a region A, the likelihood of n points (Z 1 , . . . , Z n ) observed in A is proportional to

e -µ(A) n i=1 dµ dz (Z i ).
The exponential term e -n λ ψ (Au) in (2.8) follows from the homogeneity property of λ ψ , e -([0,1]) . λ ψ (Au,n) = e -nλ ψ (Au) .

The terms J χ j are the inverse Jacobians of the marginal transformations T χ j : Y j → X j , i.e.

J χ j (y j ) = σ -1 j (ζ j ) -ξ j x 2 j e 1 x j 1 -e -1 x j 1+ξ j
, where x j = T χ j (y j ), Appendix B. Integration of the exponent measure along directions of missing components

The exact expression for the integral of the exponent measure's density in (3.2) along the axes [0, ∞] corresponding to missing coordinates is given below. Let D = {j 1 , . . . , j r } be the non missing coordinates (r < d). Integrating dλ ψ dx over (0, ∞) in the missing directions D 0 = i 1 , . . . , i d-r yields the marginal density of λ with respect to the Lebesgue measure on the vector space spanned by D, i.e.

∂λ ψ (x) ∂x j 1 ••• ∂x jr .
With a DM angular measure, the integral has an analytic expression, which is, using (2.2) and (2.3),

∂λ ψ (x) ∂x j 1 • • • ∂x jr = {z:z j =x j (j∈D),z i ∈R + (i∈D 0 )} dλ ψ dx (z) dz i 1 , . . . , dz i d-r = r k m=1    p 0 m Γ(ν 0 m ) j∈D Γ(ν 0 m µ 0 j,m ) j∈D x ν 0 m µ 0 j,m -1 j   j∈D x j   -(ν 0 m +1)    , (B.1) with ν 0 m = ν m (1 - i∈D 0 µ i,m ) , µ 0 m = (1 - i∈D 0 µ i,m ) -1 µ m , p 0 m = d r (1 - i∈D 0 µ i,m )p m . (B.2)
This is the spectral measure associated with another angular DM distribution on S r with parameter ψ 0 = (ν 0 1:k , µ 0 1:k , p 0 1:k ). The censored likelihood (3.2) can thus be re-written as

L v (O, θ) = exp -n det λ ψ (A u ) - I i=1 n i λ ψ (A i ) × • • • • • • nv i=1 [ L0 t i , R0 t i ] ∂ r(i) λ ψ ∂x j 1 (t i ) • • • ∂x j r(i) (t i ) (x) d 0 i j:Y j,t i >v j J χ j (y j,t i ) , (B.3)
where integration is performed in the censored, non-missing directions, and where 0 i is the Lebesgue measure on the corresponding subspace of R d and [ L0 t i , R0

t i ] are the original bounds [ Lt i , Rt i ], modulo the subspace of missing components.

Appendix C. Data augmentation details

Here is detailed the construction of augmentation data Z = (Z above , Z ), first introduced in Section 4.2.

C.1. Definition of augmentation the variable Z above = {Z t i } i≤nv . Consider a Fréchet-transformed, censored observation C χ t i = (C χ 1 , . . . , C χ d ), with C χ j = (κ j,t i , X j,t i , Lj,t i , Rj,t i ), as in Section 3.3. Let D c (i) = (j 1 (i), . . . , j c (i)) be the censored coordinates in observation C χ t i . The latent variables Z t i = (Z j 1 ,t i , . . . , Z j c ,t i ) are defined so as to 'replace' those coordinates: More formally, let

Xt i | ψ ∼ 1 λ ψ (A u ) 1 Au ( • ) λ ψ ( • ) , (C.1)
be an uncensored d-dimensional variable with Fréchet margins and dependence structure given by λ ψ on A u . Then, Z t i is defined through its joint distribution with the {X j,t i : j / ∈ D c (i)}, conditionally on θ,

Z t i , {X j,t i : j / ∈ D c (i)} θ distribution = Xt i |ψ.
Then, conditionally on the observation O t i ,

Z t i | O t i , θ = [( Xj 1 ,t i , . . . , Xj c ,t i ) | O t i , θ] = ( Xj 1 ,t i , . . . , Xj c ,t i ) Xt i ∈ [ Lt i , Rt i ], ψ . (C.2)
The contribution of the 'augmented' data point (Z t i , O t i ) to the augmented likelihood (4.4) is

[z t i , O t i |θ] = dλ ψ dx (x t i ) j:κ j,t i =1 J χ j (y j,t i ) , (C.3)
where xj,t i = T χ j (y j,t i ) if κj,t i = 1 , z j,t i otherwise .

Remark. With missing components D 0 (i) = {j : κ j,t i = 0} = ∅, integration in the direction D 0 (i) can be performed analytically (see Appendix B), which reduces the dimension of the augmented data. Indeed, in such a case, the corresponding Z j,t i 's need not be included in Z t i , the uncensored variable Xt i is defined on the quotient spaces E/D 0 (i) and its distribution is proportional to the exponent measure λ ψ 0 defined by equations (B.1) and (B.2), with density

∂ r(i) λ ψ ∂x j 1 (t i ) ••• ∂x j r(i) (t i ) ( • ) . as in equation (B.3), so that [z t i , O t i |θ] = ∂ r(i) λ ψ ∂x j 1 (t i ) • • • ∂x j r(i) (t i ) (x t i ) j:κ j,t i =1 J χ j (y j,t i ) , (C.4) C.2. Full conditional distribution of augmented data [ Z above | O, θ ]. The full conditionals [Z j,t i |{Z s,t i } s =j , O t i , θ
] are functions of truncated Beta distributions that can easily be sampled in a Gibbs step of the algorithm, as shown below.

In the remaining of this subsection, we omit the temporal index t i . If ψ is a mixture of k Dirichlet distributions, as in (2.3), then, for any bounded, continuous function g defined on [ Lj , Rj ], the conditional expectation of g(Z j ) is, up to a multiplicative constant, , s j = r =j x r and ρ j = r =j x νmµr,m-1 r

E [g(Z j ) |O, {Z r } r =j , θ ] = E g( Xj ) Xj ∈ [ Lj , Rj ], Xr = x r (r = j), ψ = Rj Lj g(x j )h ψ x r x r
. Changing variable with u = x j /(x j + s j ), the integration bounds are

R j = Rj s j + Rj , L j = Lj s j + Lj ,
and we have

I m = γ m ρ j R j L j g s j u 1 -u s j u 1 -u νmµ j,m -1 (s j + s j u 1 -u ) -νm-1 s j (1 -u) 2 du = γ m ρ j s -νm(1-µ j,m )-1 j R j L j g s j u 1 -u u νmµ j,m -1 (1 -u) νm(1-µ j,
D m = Γ(a m + b m ) Γ(a m )Γ(b m )IB am,bm (L j , R j ) = Γ(ν m ) Γ(ν m µ j,m )Γ(ν m (1 -µ j,m )) 1 IB am,bm (L j , R j )(1 -µ j,m )
Finally, we have

I m = C m .D m R j L j g( s j u 1 -u ) u νmµ j,m -1 (1 -u) νm(1-µ j,m ) du with C m = (1 -µ j,m ) Γ(ν m (1 -µ j,m )) s =j Γ(ν m µ s,m ) IB am,bm (L j , R j ) ρ j s -νm(1-µ j,m )-1 j , (C.6)
so that that the conditional expectation (C.5) is that of a mixture distribution,

E [g(Z j )|O, {Z r,r =j }] = j m=1 p m E g( s j U m 1 -U m ) with weights (p m ) m≤k , p m = p m C m , (C.7)
where C m is given by equation (C.6).

As a conclusion, the conditional variable

[Z j |O, Z s =j , θ] is a mixture distribution of k components p m , V j,m = s j U m 1 -U m 1≤m≤k . (C.8)
In presence of missing coordinates, (C.8) still holds, up to replacing the mixture parameters (p, µ, ν) with (p 0 , µ 0 , ν 0 ) as in equation (B.2).

Z u ∼ PP(λ ) on E u,n det , ϕ u (Z u ) = (1 -1/τ ) Z u (Au,n det ) , (C.9)
where Z u (A u,n det ) is the number of points forming Z u which hit A u,n det . Full justification and simulation details are given in the next subsection (Appendix C.4).

Let {x u,s = (r u,s , w u,s )} s≤N u be the points of Z u in E u,n det , the density of Z u over E u,n det , which contributes to the augmented likelihood (4.4), is

[z u | ψ] = 1 N u ! e -n det τ d min j≤d u j N u s=1 τ d (r u,s ) 2 h ψ (w u,s ) .
(C.10)

The processes Z i 's and the weights ϕ i 's are defined similarly, replacing n det with n i and u j with Rj,t i (1 ≤ i ≤ I ).

C.4. Consistency of the augmentation model . The augmentation process Z and the weight function ϕ have been constructed with a hint towards using the general expression of the Laplace transform of a Poisson process to prove consistency of the augmented posterior.

Proposition 1. Define the factors

[ z t i , O t i | θ ] (i ≤ n v ) and [ z u | ψ], [z i | ψ] (1 ≤ i ≤ I )
composing the augmented likelihood (4.4), according to equations (C.3) and (C.10), and let the factors ϕ u , ϕ i , i ≤ I of the weight function ϕ in (4.2) be defined by equation (C.9).

Then, the augmented posterior

[ z, θ | O ] + ∝ [ θ ][z, O | θ ] ϕ(z) is consistent, in
the sense that condition (4.3), whence (4.1), is satisfied.

Proof .

It is enough to show that, on the one hand,

[z above , O | θ ] dz above = nv i=1 [ Lt i , Rt i ] dλ ψ dx d i (x) j:y j,t i >v j J χ j (y j,t i ) , (C.11)
and, on the other hand,

E [ϕ u (Z u )] = exp (-n det λ ψ (A u )) , E [ϕ i (Z i )] = exp (-n i λ ψ (A i )) (i ≤ I ) . (C.12)
where the above expectations are taken with respect to

[Z u | ψ ] and [Z i | ψ ]. Establishing (C.11) is immediate: from definition (C.3) of the [z t i , O|θ]'s com- posing [z above , O|θ], [ Lt i , Rt i ] [z t i , O t i |θ] dz t i = [ Lt i , Rt i ]
dλ ψ dx d i (x) j:y j,t i >v j J χ j (y j,t i ) , which yields (C.11) by taking the product over indices 1 ≤ i ≤ n v . It remains to show (C.12). To wit, ϕ u is a smoothed version of the indicator function 1 {N u (Au,n det )=0} , which expectancy is P(N u (A u,n det ) = 0) = e -n det λ ψ (Au) , as soon as N u is a Poisson process with intensity measure λ ψ . For f a bounded, continuous function defined on a nice space E and N a point process on E, denote

N (f ) = E f dN = N (E) i=1 f (s i ).
Then, if N is a Poisson process PP(λ) on E, the Laplace transform Lap N (f ) =E(e -N (f ) ) is (Resnick, 1987, Chap. 3)

Lap N (f ) = exp - E (1 -e -f (s) ) dλ(s) .
Consider the region E = E u,n det as in (C.9) and take

f u (x) = -log(1 -1/τ )1 Au,n det (x) , so that 1 -e -fu = 1 τ 1 Au,n det .
With these notations,

ϕ u (Z u ) = 1 - 1 τ Z u (Au,n det ) = exp -Z u (f u ) , whence E ϕ u (Z u ) = Lap Z u (f u ) = exp - Eu,n det (1 -e -fu ) dλ = exp - Eu,n det 1 τ 1 Au,n det d(τ λ ψ ) = exp (-λ ψ (A u,n det )) .
This shows the first equality in (C.12). The second one is derived with a similar argument.

The points of Z u can easily be simulated (see Resnick, 1987, Chap.3): the number of points N u in E u,n det is a Poisson random variable with mean equal to

λ (E u,n det ) = τ d min j (u j /n det ) ,
and each point has density in polar coordinates equal to

1 λ (Eu,n det ) τ d r 2 h ψ (w).
Remark. One may be tempted to define Z u as a Poisson process with intensity λ ψ on some E ⊃ A u,n det , and ϕ(Z u ) as the indicator 1 z u (Au,n det )=0 , with a similar definition for the ϕ i 's and the Z i 's. As pointed out in the proof of Proposition 1, one would have E 1 Z u (Au,n det )=0 = exp (-n det λ ψ (A u )), as required. However, even if this construction is valid in theory, it leads to a very large rate of rejection in the Metropolis algorithm: ϕ(Z u ) has too much variability around its mean value and the proposal is systematically rejected each time a point in the augmentation process hits the failure region. C.5. Expression of the augmented posterior. Recall from Section 4.2, equation (4.2), that the augmented posterior density to be sampled by the MCMC algorithm is

[ z, θ | O ] + ∝ [ θ ][ z, O | θ ]ϕ(z) .
Combining equations (C.3) and (C.9), and integrating out missing components as in Appendix B, the developed expression is

[ z, θ | O ] + ∝ [ θ ] [z u |ψ] I i=1 [z i |ψ] • (1 -1/τ ) z u (Au,n det )+ I i=1 z i (A i,n i ) • • • • • • nv i=1 ∂ r(i) λ ψ ∂x j 1 (i) • • • x j r(i) (i) (x t i ) j:Y j,t >v j J χ j (y j,t ) , (C.13)
where [z u |ψ] and the [z i |ψ]'s are given by equation (C.10).

Appendix D. MCMC algorithm

The MCMC algorithm generates a sample (θ ι , Z ι ) ι=1,...,N which distribution converges to the invariant distribution of the chain, which is the augmented posterior distribution [Z, θ|O] + , as defined in Section 4 and Appendix C.5. The quantity of interest here is the joint parameter θ, which is the concatenation of the marginal parameters and the dependence parameter: θ = (χ, ψ). We recall that MCMC algorithms aiming at sampling a quantity ∆ ∈ E according to a density π( • ) proceed typically as follows • Start with any value δ(0) ∈ E • for ι ∈ {1, . . . , N } (1) Generate δ * according to a proposal distribution with density q(δ(ι), • )

(2) Compute the acceptance ratio

α = π(δ * ) q(δ * , δ(ι)) π(δ(ι)) q(δ(ι), δ * ) ,
and generate U , a uniform random variable on [0, 1].

(3) If U > α, 'reject' the proposal and set δ(ι + 1) = δ(ι) .

Otherwise (i.e. with probability α), set δ(ι + 1) = δ * (4) Set ι = ι + 1, go to (1).

• Return δ(ι min ), . . . , δ(N ) , where ι min is the length of the burn-in period, after which the chain is deemed to have reached a stationary behavior. In our case, the unnormalized posterior measure [Z, θ|O] + (see equation (C.13)) on the augmented parameter space plays the role of the objective density π above.

In a Metropolis-within-Gibbs MCMC, several proposal kernels q 1 , . . . , q T are defined, each of them corresponding to a type of move, which is randomly chosen among {1, . . . , T } at each iteration ι and allows to modify some subset of components in δ alone. The algorithm developed in this paper builds on the MCMC algorithm proposed by [START_REF] Sabourin | Bayesian dirichlet mixture model for multivariate extremes: A re-parametrization[END_REF], in which several types of move modifying the dependence structure (Dirichlet mixture parameter ψ) have been defined. Those are kept as they are in the present work, the novel part of which concerns the definition of marginal moves (modifying the marginal parameter χ) and augmentation moves (modifying the augmentation data Z). Additional notations distinguishing between the quantities appearing in the general MCMC algorithm above, according to the type of move, are omitted in the remainder of this section. D.1. Starting values. In a preliminary step, likelihood optimization is performed in the independent model (the likelihood for one multivariate observation is the product of d Pareto densities). This provides starting values for the marginal parameters as well as a Hessian matrix H, that may be used as the inverse of a reference covariance matrix when updating the marginal parameters. D.2. Marginal moves. The marginal parameter χ is updated as a block: The proposal is normal, with mean at χ(ι) and co-variance matrix Σ = δ H -1 , where δ is a scaling factor fixed by the user, that may typically be set around 0.5 and H is the Hessian matrix computed in the preliminary step. Since the proposal density is symmetric, and since this move does not modify the dependence structure, neither the terms involving the proposal density, nor the point processes Z , appear in acceptance ratio. The augmented variables Z above are left unchanged. If any augmented component Z j,t i is outside of the candidate censoring interval [ L * t i , R * t i ] (on the new Fréchet scale) resulting from the modification of the marginal parameters, the move is rejected, since in such a case, the candidate has augmented likelihood [z, O | θ * ] = 0, i.e. α = 0.

Otherwise, the uncensored Fréchet-transformed variables (such that κj,t i = 1) are updated to X * t i = T χ * j,t i (Y j,t i ). The acceptance ratio is

α = [χ * ] [χ(ι)] nv i=1 ∂ r(i) λ ψ ∂x j 1 (i) • • • ∂x j r(i) (i) ( X * t i ) ∂ r(i) λ ψ ∂x j 1 (i) • • • ∂x j r(i) (i)
( Xt i (ι))

-1 j:Y j,t i >v j J χ * j (Y j,t i ) J χ(ι) j (Y j,t i ) , where j 1 , . . . , j r(i) are the non-missing components in the censored observation C χ t i and Xj = X j for uncensored components, Xj,t i = Z j,t i otherwise (see Appendix B). D.3. Augmentation moves for Z above . The augmented components {Z j,t i } = { Xj,t i : κj ∈ {2, 3}} (c.f. Section 4.2) are re-sampled, one coordinate at a time, from their exact conditional distribution given the other coordinates, as derived in Appendix C.1. Since no other component of (Z, θ) is modified, the proposal density equals the objective density, and the acceptance ratio is thus set to α = 1. [(Z u ) * (Au,n det )-Z u (ι)(Au,n det )]+ i≤I (Z i ) * (A i,n i

)-Z i (ι)(A i,n i )

.

D.5. Dependence moves. These types of moves allow to update ψ(ι).The only difference between the present algorithm and what is described is [START_REF] Sabourin | Bayesian dirichlet mixture model for multivariate extremes: A re-parametrization[END_REF] is that not enough exact angular data are available to construct proposals for moving or splitting a Dirichlet mixture components µ m . Indeed, most of the observations have at least one coordinate missing or censored, so that no 'angle' is available. Consequently, the latter proposal is a simple Dirichlet distribution with mode at µ m (ι), with re-centering parameter 0 < < 0.5, q(µ m (ι), • ) = diri d ,γ * ( • ) , with γ * = (1 -) µ m + ( 1 d , . . . , 1 d ). Each dependence move (except for a shuffling move which only affects the representation of the angular distribution, see [START_REF] Sabourin | Bayesian dirichlet mixture model for multivariate extremes: A re-parametrization[END_REF]) is systematically followed by an augmentation move updating the Poisson processes Z , which improves the chain's mixing properties. This also avoids the computation of the 'costly' term involving the density [z |ψ] (see equation (C.10)). Indeed, the acceptance ratio for the two consecutive moves (dependence move followed by a augmentation move) is

α = [ψ * ] [ψ(ι)]
q(ψ * , ψ(ι)) q(ψ(ι), ψ * )

× • • • • • • (1 -1/τ ) [Z * u (Au,n det )-Z u (ι)(Au,n det )]+ i≤I Z * i (A i,n i )-Z i (ι)(A i,n i ) × • • • • • • nv i=1 ∂ r(i) λ ψ * ∂x j 1 (i) • • • ∂x j r(i) (i) ( Xt i (ι)) ∂ r(i) λ ψ(ι) ∂x j 1 (i) • • • ∂x j r(i) (i) ( Xt i (ι)) -1
. D.6. MCMC settings and convergence diagnostics in the simulation study.

For the simulation study, the prior on the Dirichlet mixture distributions is specified in a similar way as in [START_REF] Sabourin | Bayesian dirichlet mixture model for multivariate extremes: A re-parametrization[END_REF]. The number k of mixture components has truncated geometric distribution, [k] ∝ 1 -1 λ k-1 1 λ 1 [1,kmax] (k) with upper bound k max = 10 and mean parameter λ = 4. Also, for the sake of simplicity, all the marginal parameters are assumed to be a priori independent, with normal distributions (after log-transformation of the scales). The shape parameter has standard normal distribution and the logarithms of the scales have mean and standard deviation both equal to 5.

As for the augmentation Poisson process data, the multiplicative constant τ involved in the Poisson intensity is set to 50. It appeared that smaller values of τ (close to 1) considerably affected the mixing properties of the chains.

Convergence of the dependence parameters ψ(ι) can be monitored using functionals based on integration of the simulated densities against Dirichlet test functions (see Sabourin and Naveau, 2014, for details). To detect possible mixing defects, six chains of 10 6 iterations each are run in parallel. Standard convergence diagnostic tests are implemented in R [START_REF] Heidelberger | Simulation run length control in the presence of an initial transient[END_REF][START_REF] Gelman | Inference from iterative simulation using multiple sequences[END_REF], respectively testing for non-stationarity and poor mixing. For example, the stationarity test detects three non-stationary chains out of six for the simulated data set exemplified in Section 5.2. The mixing properties of the three retained ones, as measured by a variance ratio inter/intra chains, are satisfactory enough: all the potential scale reduction factors [START_REF] Gelman | Inference from iterative simulation using multiple sequences[END_REF] are below 1.1. The same is true of the marginal parameter component of the chains, (χ(ι)) ι .

Figure 2 .

 2 Figure 2. Example in the two-dimensional case of marginal transformation and censoring below threshold, for data of marginal types (1, 1), (upper panel) (3, 1) (middle panel) and (3, 3) (lower panel). In these three cases, observations are represented in black, respectively by a dot, an horizontal arrow and a rectangle. The Grey areas represent respectively the multivariate threshold v (left side) and the Fréchet transformed one u (right side). The two upper panels correspond to observations above threshold, while the lower panel shows an undetermined observation.

4. 2 .

 2 Data augmentation in the Poisson model. In our context, finding an augmentation random variable Z + with easily manageable conditional distributions, such that the augmented posterior [z, θ| O] + satisfy (4.1), is far from straightforward, mainly due to the exponential terms (3.4) in the likelihood. Instead, an intermediate variable Z is introduced, which plays the role of a proposal distribution in the MCMC algorithm. Z is defined conditionally to θ, through the augmented likelihood [ z, O | θ ], so that the full conditionals [Z | θ, O]

4. 3 .

 3 Implementation of a MCMC algorithm on the augmented space. This section describes only the main features of the algorithm, more details are provided in Appendix D. The present algorithm is an extension of Sabourin and Naveau (2014)'s one, who proposed a Metropolis-within-Gibbs algorithm to sample the posterior distribution of the angular measure in a POT framework (2.1), within the Dirichlet mixture model (2.3). The number of components in the Dirichlet mixture is not fixed in their model, and the MCMC allows reversible-jumps between parameters sub-spaces of fixed dimension, each corresponding to a fixed number of components in the Dirichlet mixture. Their algorithm approximates the posterior distribution[ψ|W]

Figure 3 .

 3 Figure 3. Average number of exactly observed coordinates above threshold in each recorded multivariate excess in the real data set (hatched bars) and in a simulated data set (Gray bars).

Figure 5 .

 5 Figure 5. Bi-variate angular predictive densities (thick black line) with posterior credible sets (Gray regions) corresponding to the posterior 0.05 -0.95 quantiles, together with the true angular density (dashed line).

Figure 7 .

 7 Figure 7. Normalized QR scores in the dependent Dirichlet-Poisson model versus the independent model, at the four different locations, for 50 simulated datasets

Figure 8 .

 8 Figure 8. Bi-variate angular measure for the hydrological data set: posterior predictive density (black line) and posterior 0.05 and 0.95 quantiles.

  (C.1) for the definition of Xj ). Each term I m (m ≤ k for a mixture of k components) isI m = γ m )-( r≤d (νm µr,m-1)) + x j ) -νm-1 dx j ,

  m ) du One recognizes in the integrand the unnormalized density of a Beta random variable U m ∼ beta(a m , b m ), with a m = ν m µ j,m , b m = ν m (1 -µ j,m ) + 1 ; Let IB a,b (x, y) denote the incomplete Beta function (i.e. the integral of the Beta density) between truncation bounds x and y. The missing normalizing constant in the integrand is

C. 3 .

 3 Augmentation Poisson process Z = {{Z i } i≤I , Z u } and weight function ϕ. Let us define a region E u,n det = {x ∈ (R + ) d : x 1 > min j ( u j n det )}, so that A u,n det ⊂ E u,n det . Choose a multiplicative constant τ > 0 and define a Poisson intensity measure λ ( • ) = τ λ ψ ( • ). The augmentation process Z u is a Poisson processes which is defined together with ϕ u by

D. 4 .

 4 Augmentation moves for Z . During this move, proposals Z * = Z * u , Z * i , i ≤ I , for the augmentation Poisson processes introduced in the end of Section 3.4, are sampled under their exact distribution, q Z (ι), Z * = [Z * |ψ(ι)] .The latter is determined by their intensity measure λ = τ λ ψ : the multiplicative constant τ and the sampling procedure have been described in Appendix C.3. The acceptance ratio is thusα = [ Z * , θ(ι) | O ] + [ Z (ι) | ψ(ι) ] [ Z * , θ(ι) | O ] + [ Z * | ψ(ι)] all the terms cancel out except the ratio ϕ(Z * )/ϕ(Z (ι)), so that α = (1 -1/τ )

  precise definition of Z above is given in Appendix C.1, together with the expression of the corresponding contributions to the augmented likelihood, [ z t i O t i |θ ] . The full conditionals [Z above |O , θ ] are derived in C.2. The augmentation Poisson processes Z , together with the weight ϕ are defined in Appendix C.3, and the compatibility condition (4.3) is proved to hold in Appendix C.4.

preprint available on https://hal.archives-ouvertes.fr/hal-01087687

Acknowledgments

Part of this work has been supported by the EU-FP7 ACQWA Project (www.acqwa.ch), by the PEPER-GIS project, by the ANR-MOPERA project, by the ANR-McSim project and by the MIRACCLE-GICC project. The author would like to thank Benjamin Renard for providing the hydrological data that motivated this work and for his useful comments, and Anne-Laure Fougères and Philippe Naveau for their advice and interesting discussions we had during the writing of this paper.