Aurore A Lavigne
email: aurore.lavigne@univ-lille3.fr

Nicolas Eckert

Liliane Bel

Éric Michaël Deschâtres

E Parent

Modelling the spatio-temporal repartition of right-truncated data: an application to avalanche

Keywords: Bayesian hierarchical model, geostatistic, snow avalanche, truncated data

Introduction

Hazard maps are useful tools to localize the danger engendered by extreme avalanches in mountains [START_REF] Jamieson | Application and limitations of dynamic models for snow avalanche hazard mapping[END_REF][START_REF] Barbolini | A new method for avalanche hazard mapping using a combination of statistical and deterministic models[END_REF].

According to a reference hazard related to a set of suitable variables, at risk zones are dened by combining information about avalanche frequency and intensity. Whereas frequency is often considered stationary and rather easily estimated, mapping intensity is the crux of the issue [START_REF] Ancey | Computing extreme avalanches[END_REF][START_REF] Eckert | Bayesian stochastic modelling for avalanche predetermination: from a general system framework to return period computations[END_REF]. Two main approaches are encountered in the literature to evaluate avalanche intensity, namely, statistical methods and dynamic methods.

The former are based on direct observations of avalanche intensity, commonly measured by runout distance, i.e. the distance covered by the avalanche. In dynamical approaches, runout distributions are derived from factors at the origin of avalanches, such as the accumulated snow depth in the release area, which are propagated according to dynamical models [START_REF] Barbolini | Application of statistical and hydraulic-continuum dense-snow avalanche models to ve real European sites[END_REF].

However, these methods suer from the small amount of data avalaible, and from many uncertainties in physical laws driving the motion of snow [START_REF] Jamieson | Application and limitations of dynamic models for snow avalanche hazard mapping[END_REF].

Focusing on a region rather than on a path has a major advantage: observations coming from multiple paths can be gathered. However, runout distances cannot be directly compared between paths, since they depend on local climatic and topographical factors. As a consequence, data standardisation is needed. Within the common alpha-beta [START_REF] Lied | Empirical calculations of snow-avalanche run-out distance based on topographic parameters[END_REF] and runout ratio [START_REF] Mcclung | Statistical and geometrical denition of snow avalanche runout[END_REF] methods, each runout distance is transformed into a transferred variable, as a function of typical slope values from the paths where it occurred. All transferred variables are then assumed to be independent, identically distributed with their sample size large enough for robust inference. While these methods benet from all regional information, they may be ineective when the local path topography is irregular [START_REF] Smith | Characteristics and prediction of high-frequency avalanche runout[END_REF]. In addition, they do not make use of any notion of proximity.

Even if, to our knowledge, the spatial dependance of avalanche intensity has never been demonstrated from a statistical point of view, some hints may suggest its reality. Indeed, avalanches arise from topographic and weather factors that are both spatially distributed. For instance, studying the spatial repartition of extreme snowfalls [START_REF] Gaume | Mapping extreme snowfalls in the French Alps using max-stable processes[END_REF] and snow depths [START_REF] Blanchet | Spatial modeling of extreme snow depth[END_REF], two factors impacting directly avalanche activity, have shown that they are spatially dependent. Moreover, the existence of signicant spatial dependence in avalanche frequencies in the northern French Alps has been shown in [START_REF] Eckert | Hierarchical Bayesian modelling for spatial analysis of the number of avalanche occurrences at the scale of the township[END_REF][START_REF] Eckert | A spatio-temporal modelling framework for assessing the uctuations of avalanche occurrence resulting from climate change: application to 60 years of data in the northern French Alps[END_REF].

Bayesian hierarchical spatio-temporal models (BHM) [START_REF] Cressie | Statistics for spatio-temporal data[END_REF] are popular tools to analysis complex spatio-temporal structure on many kinds of data; they are encountered in various domains such as epidemiology [START_REF] Best | A comparison of Bayesian spatial models for disease mapping[END_REF], ecology [START_REF] Munoz | Estimation and prediction of the spatial occurrence of sh species using Bayesian latent Gaussian models[END_REF] and environment [START_REF] Vanem | A Bayesian hierarchical spatio-temporal model for signicant wave height in the North Atlantic[END_REF]. In the BHM framework, observations are modelled conditionally to the spatio-temporal process and the conditional denition allows for the splitting of the model into three layers: observation -process -prior. The observa-tion model takes into account data specicities, whereas the spatio-temporal process is generally modelled with the convenient Gaussian process. Under the Bayesian framework, inference is generally straightforward thanks to the Gibbs algorithm (Gilks et al, 1995b).

The main issue in applying the framework to our case arises from data truncation. Indeed, in the French avalanche database we work with, avalanches are recorded only if their runout altitude is lower than a given threshold. This special case diers from the more common spatial censored data (Stein, 1992), because several repetitions of truncated data are supplied for each site (for us, the avalanche path), and the proportion of non-observed data is unknown and depends on the site. In such a case, it is known that ignoring the truncation will produced biased estimators (A 'Hearn, 2004), and so a conditional distribution should be assigned to observations. However, parameter estimation under this distribution is uneasy, because truncation distorts the likelihood surface to create large at zones of instabilities. As a consequence, the variance of maximum likelihood estimators is high, and modied estimators obtained with penalised likelihood [START_REF] Cope | Penalized likelihood estimators for truncated data[END_REF], or restricted likelihood (A 'Hearn, 2004) have been developed. These two methods aim at introducing a small bias in order to decrease the estimator's variance. Under the Bayesian framework, [START_REF] Zhou | Bayesian estimation of truncated data with applications to operational risk measurement[END_REF] report a good behaviour of the Bayesian estimates, even if non-informative Jerey's priors are used.

In this paper, we develop a spatio-temporal BHM to integrate within a unied framework, all the information about avalanche intensity in Haute-Savoie, the northern most French Alps department. On the observation layer the truncated observations are modelled with a conditional distribution. We reduce the instabilities linked to this distribution by restricting the parameter space to a reasonable domain, with constraints on the mean inter-annual runout altitudes, given by physical laws or expert knowledge. On the process layer, the spatio-temporal covariance is assumed to be separable, and we model the spatial dependence between mean inter-annual runout altitudes with a Matérn covariance function. As a result, we show that spatial dependence exists in runout altitudes. We then take advantage of the identied spatial pattern to predict avalanches on non-documented paths.

The data, their origins and their specicities are described in Section 2. Our specically designed BHM with a Gaussian truncated distribution as model for observations is described in Section 3, as well as the constraints on parameters. This model is inferred under the Bayesian paradigm. In Section 4, we discuss two points concerning Bayesian inference: the prior choice and the specicities of the Gibbs sampler we propose. In Section 5, we display results regarding the application on avalanche runout altitude. Finally, in Section 6 we explore the predictive abilities of our approach.

Avalanche data in Haute-Savoie

A systematic survey of avalanche events began in France, in the early 20 th cen- tury under the scope of the EPA (Enquête Permanente des Avalanches ). The aim of the EPA survey was to better understand avalanche activity in order to limit damage caused to forests, roads and inhabited areas. Consequently, a set of paths was chosen and the runout altitude of avalanches occurring on these paths has constantly been recorded since then. EPA data are accessible from the website http://www.avalanches.fr/epa_presentation-donnees/.

In 2002, an observation threshold was further dened for each path. It corresponds to a level above the bottom of the path, such that only the avalanches running down below it are recorded. The thresholds are chosen to be easily detectable by observers such as roads and rivers, but also such that a substantial part of the total avalanche activity, and at least the avalanches arriving close to potential elements at risk, are recorded. Table 1 displays the distribution of the altitudes between valley altitudes and thresholds. On average, the thresholds are located 300 m above the valley altitude, but this distance varies from one path to another within a range of ∼1000 m. Before 2002, the threshold was implicitly used by most of the observers and we do not know the proportion of recorded and non-recorded avalanches above the threshold.

Consequently, here, over the entire of this study period, we consider only the avalanches with a runout altitude below the threshold and treat the data as right-truncated.

We have at our disposal a database of n = 5331 avalanches that have been recorded in Haute-Savoie during the 1925-2012 period, indexed by year and path. They occurred on 389 paths, among the 553 paths currently surveyed by the EPA in Haute-Savoie. These 389 paths are located in the eastern part of the department where altitudes are on average high, and most of them are aligned on both sides of valleys, as shown in Fig. 1. In particular, more than 90 paths are concentrated along the Arve valley, where the well-known township of Chamonix-Mont-Blanc is situated, and where avalanche danger is known to be very high [START_REF] Naaim | Return period calculation and passive structure design at the Taconnaz avalanche path, France[END_REF][START_REF] Naaim | Dense avalanche friction coecients : inuence of physical properties of snow[END_REF].

Table 1 provides some indicators of the path features, emphasizing the heterogeneity of the paths. Concerning the avalanche record frequency, fewer than 5 avalanches are recorded for 41 % of the paths, but for 6 % of them, the number of records exceeds 50. The within path avalanche activity is also heterogeneous: the standard deviation (sd) of recorded runout altitude is high for some paths, but it generally does not exceed 50 m (for 61 % of the paths).

Specically, for 42 paths, the standard deviation is null, meaning that all the avalanches have exactly the same runout altitude. This result is not surprising, given that the precision of measurement is at least 10 m and surely larger in the past. More generally, this heterogeneity in local runout variability expresses the large range of topographies in the sample of studied paths. Indeed, regular or moderately varying path proles lead to low variability of runout altitudes, whereas for paths with more complex topography the variability can ]0, 100 [ [100, 200[ [200, 400[ [400, 600[ [600, 1020[ ]0, 15 [ [15, 30[ [30, 50[ [50, 75[ [75, 100[ [100, 250[ Nb. 65 42 34 59 65 74 26 24 be much higher, as demonstrated by long range runout simulations [START_REF] Meunier | Towards a conceptual approach to predetermining long-return-period avalanche run-out distances[END_REF][START_REF] Eckert | Bayesian optimal design of an avalanche dam using a multivariate numerical avalanche model[END_REF].

All in all, the avalanche altitude data at our disposal is rather impressive, but with many particular features, which have to be taken into account in an ecient modelling.

3 Hierarchical framework for kriging truncated data 3.1 A spatio-temporal model for right-truncated data

We make use of a hierarchical specication to take into account all the data features at the observation level, and to easily model space and time uctuations of runout altitudes on a Gaussian latent layer [START_REF] Cressie | Statistics for spatio-temporal data[END_REF].

We opt for Gaussian truncated distributions as an observation model to take into account that the data are observed below a threshold; this choice is discussed in Section 5.5. In addition to the truncation, we chose an heteroscedastic model, i.e. path-specic variances in order to tackle the large between path variability of empirical standard deviations.

We denote Y cti the observed runout altitude for avalanche i

∈ {1 • • • n ct }, for the path c ∈ {1, • • • n c }, during the year t ∈ {1 • • • n t },
s c the observation threshold for path c, µ ct the mean of the complete (that is not truncated) Gaussian distribution and σ 2 c its variance. We assume that conditionally to µ ct and σ 2 c , the variables Y cti are independent and their probability density function (pdf ) is expressed as

f (y cti |µ ct , σ 2 c ) = 1 Φ sc-µct σc 2πσ 2 c exp -1 2σ 2 c (y cti -µ ct ) 2 1 ycti<sc , (1)
where Φ() is the cumulative density function (cdf ) of the standard normal distribution.

On the latent layer, we aim at modelling the spatio-temporal uctuations of the runout altitude. Our choice, discussed in Section 5.5 and based on an empirical analysis of covariance matrices (Appendix A), is made on a simple spatio-temporal factorial decomposition. The annual mean by path µ ct is supposed to be the sum of an annual term B t and a spatial term C c . For symmetry purposes between space and time terms, we assume that c C c = t B t = 0 and introduce a constant term α, and so µ ct can be then written as

µ ct = α + B t + C c . (2)
The temporal term B t is supposed to be the sum of a second order random walk g t of variance δ 1 and an independent Gaussian noise t of variance δ 0 , re-written as B t |g t , δ 0 ∼ N (g t , δ 0 ).

(3) Under the Bayesian paradigm, the posterior mean of g t approximates a smoothing spline [START_REF] Wahba | Improper priors, spline smoothing and the problem of guarding against model errors in regression[END_REF]. This non-parametric smooth term has been proved to be convenient to model large-scale inter-annual variability, for instance in avalanche data [START_REF] Lavigne | A model for spatio-temporal clustering using multinomial probit regression: application to avalanche counts[END_REF], or related proxies [START_REF] Rabatel | Changes in glacier equilibrium-line altitude (ELA) in the western Alps over the 1984-2010 period: evaluation by remote sensing and modeling of the morpho-topographic and climate controls[END_REF].

The spatial term aims at reproducing the local variations and the spatial dependence between paths of the mean runout altitude. For the sake of simplicity, we denote as well c the centroid of the contour line of the path c. We assume that C is a Gaussian random eld that can be written

C c = µ c + A c , (4)
where A c is a zero mean Gaussian random eld with a Matérn isotropic covariance function [START_REF] Matérn | Spatial variation[END_REF]. Thus, the covariance between A c and A c is given by

cov(A c , A c ) = τ 2 1 2 ν-1 Γ (ν) h cc φ ν K ν h cc φ if c = c ρ 2 + τ 2 otherwise, (5)
where h cc is the Euclidean distance between the centroids c and c , K ν is the modied Bessel function of the second kind of order ν > 0 and Γ is the gamma function. The Matérn class of covariance functions is recommended for its exibility [START_REF] Banerjee | Hierarchical modeling and analysis for spatial data[END_REF]. Indeed, in addition to the sill τ 2 , the nugget ρ 2 and the range φ, the Matérn covariance function is also parametrised by the smoothness parameter ν. This parameter is related to the smoothness of the spatial process at small distances: the larger it is, the smoother the process. In particular, when the smoothness parameter is 0.5 the Matérn covariance function is the exponential covariance function, and in the extreme case where ν → ∞, the Matérn covariance is the Gaussian covariance. As it is well-established that this parameter is dicult to infer, we consider that it belongs to the set {0.5, 1, 1.5}. We do not consider higher values, because in practice, data do not carry enough information to distinguish ν = 2 from ν > 2 [START_REF] Banerjee | Hierarchical modeling and analysis for spatial data[END_REF].

Two covariates are used to model the mean spatial term µ c , namely, the altitude valley h c and the release area exposure o c , so that the mean µ c is expressed as

µ c = a + bh c + o c . (6)
o c ∈ {o N , o S } is a categorical variable with two levels, S for the south, southeast and southwestern directions and N for the other directions. It is intended to capture the fact that south facing slopes behave dierently because of the inuence of sun exposure on snow (melting, metamorphism, etc.). For identiability purposes, we set o N = 0.

Constraints on the mean runout altitude

In contrast to the ususal case, maximum likelihood estimators for the mean and variance of a Gaussian truncated sample are not independent. Thus, the likelihood surface for both parameters is distorted and presents a large at zone shown in Fig. 2. A couple of estimates close to the maximum likelihood estimate is as likely as a couple of estimates far from the likelihood maximum. As a consequence, estimation algorithms may fail to converge. We thus add constraints on the mean runout altitude to restrict the parameter space to reasonable values. More specically, we have at our disposal the following information:

at least 50% of the data are observed for each path, in accordance with the threshold selection;

the valley bottom is a lower bound for all avalanches, hence the mean runout altitude is higher than (or ultimately equal to) the valley altitude.

These assertions can be expressed as constraints on the sum α + C c . The rst can be written q Yc 0.5 < s c , with q Yc 0.5 the 50% quantile of the Gaussian distribution N (α + C c , σ 2 c ). It means that the runout altitude mean by path α + C c must not be higher than the observation threshold s c . The second is simply written h c ≤ α + C c . These constraints are not very restrictive, but they are essential to make inference tractable.

Predictions

The spatial dependence in the mean runout altitude helps in making predictions for a set of undocumented paths c 0 . Since a dierent variance is set by path, we assume that the variance distribution for the predicted paths σ 2 c0 is the mixture of the posterior distributions of the variances of the mod- elled path set. We also chose to make predictions for the entire study period, and thus we integrate over time. Basically, we krige the spatially dependent Prior distribution assignment to the parameters is a sensitive task since they must reect the a priori knowledge about the parameters [START_REF] Asa Carlin | Bayes and empirical Bayes methods for data analysis[END_REF][START_REF] Gelman | Bayesian Data Analysis[END_REF]. Here, we use both non-informative and informative prior distributions elicited from expert knowledge. We assume that parameters are independent a priori, and most of the time we make recourse to the convenient conjugate priors.

term A = (A 1 , • • • , A nc )
Non-informative vague priors are assigned to the regression parameters, α prior is taken improper, and a zero mean Gaussian distribution with variance 10 7 is chosen for a, b and o S . For the variance components of the temporal part, the improper density 1/ √ δ 0 δ 1 is selected, as advised by [START_REF] Speckman | Bayesian nonparametric regression and autoregression priors Stein ML (1992) Prediction and inference for truncated spatial data[END_REF]. Independent inverse gamma distributions with 0.0001 for rate and shape parameters are assigned to τ 2 and ρ 2 . A uniform distribution [0, φ max ] is assigned to the range parameter φ. Since, φ is highly dependent on the smoothness parameter ν, φ max is selected according to the value of ν considered, such that the eective range, i.e. the distance beyond which the covariance function does not exceed 5% of the sill , follows a uniform [0, 100 km] distribution.

Since 100 km is about the distance between the two furthest sites of the study region, this prior can be considered as non-informative.

The n c variances σ 2 c , c ∈ {1, • • • n c } are assumed to be independent and inverse gamma distributed. To introduce extra data information and help ensure the stability of the algorithm, informative priors are elicited. The hyperparameters of each gamma distribution are set such that the standard deviation mean and variance equal 25 meters, a value considered reasonable by experts. This distribution is then right and left truncated to ensure that no unrealistic variance can be sampled; the standard deviation σ c must not be smaller than 30 m, representing the observation error, and not larger than 316 m.

A Gibbs sampler to draw from the posterior

Relying on the conditional independence relations in the BHM, for each value of the smoothness parameter ν, the joint posterior distribution can be written as

[ α, C, B, σ 2 c , g, a, b, o S , A, δ 0 , δ 1 , φ, τ 2 , ρ 2 |Y, ν] ∝ nc c=1 nt t=1 nct i=1 [Y cti |α, B t , C c , σ 2 c ] × nc c=1 [C c |a, b, o S , A c , ρ 2 ] × ( nt t=1 [B t |g t , δ 0 ]) ×[A|τ 2 , φ, ν] × [g|δ 1 ] ×[α, σ 2 c , δ 0 , δ 1 , τ 2 , ρ 2 , φ, a, b, o S ], (7)
which is a decomposition of the conditional equations dening the observation model (Eq. 1), the process model (Eqs. 2, 3, 4) and the prior model. A Gibbs sampler is implemented in C++ to draw parameters and latent variables from their joint posterior distribution. In case of a standard linear regression, the algorithm would be straightforward. Here, the truncation, the constraints on the parameters and the spatial covariance function complicate the inference.

For the variance components of the spatial term, we rewrite the algorithm supplied by Diggle and Ribeiro ( 2007) and use a Metropolis-Hastings step to sample the range parameter.

For dealing with the truncation, we follow Griths ( 2004)' algorithm which proposes to add latent variables Z cti sampled in the non-truncated distribution N (µ ct , σ 2 c ). These latent variables can be interpreted as pseudo-data in the case of non-truncation, and conditionally to µ ct and σ 2 c they are linked to Y cti , by the relationship P (Z cti < z cti ) = P (Y cti < y cti ) leading to

z cti = µ ct + σ c Φ -1 Φ( ycti-µct σc ) Φ( sc-µct σc ) . ( 8
)
It is then straightforward to sample µ ct and σ 2 c conditionally to Z cti . However, µ ct and σ 2 c are correlated to Z cti , and to decrease the dependence between µ ct and σ 2 c into the MCMC sample, we sample directly σ 2 c from Y cti and not from Z cti . The posterior conditional distribution is then given by

[σ 2 c |Y, α, C c , B] ∝ nt t=1 nct i=1 1 Φ( Y cti -µ ct σc ) 1 σ 2 c a0+ 1 2 mc+1 × exp -1 σ 2 c b 0 + nt t=1 nct i=1 (y cti -µ ct ) 2 , ( 9
)
where m c is the number of avalanches which occurred on path c, and where a 0 and b 0 are the hyperparameters of the inverse gamma prior distribution. We draw σ 2 c according to the pdf of Eq. 9 thanks to the adaptive rejection Metropolis sampling (ARMS) algorithm (Gilks et al, 1995a), an adaptation of the adaptive rejection sampling [START_REF] Gilks | Adaptive rejection sampling for Gibbs sampling[END_REF] suitable when the pdf is nearly log-concave.

Another diculty comes from the constraints assigned to the parameters, which are of two types. Firstly, the identiability constraints stipulate that nt t=1 B t = 0 and nc c=1 C c = 0. A deterministic relation links the components of vectors B, respectively C, and we have to draw them in subspaces of dimension n t -1, respectively n c -1. For this, we follow the recommendations of [START_REF] Gelman | Analysis of variancewhy it is more important than ever[END_REF]

and sample B 1 , • • • , B nt-1 conditionally to nt t=1 B t = 0, then we set B nt = - nt-1 t=1 B t .
Secondly, the mean runout altitude by path, i.e. α+C c , is bounded. Specically, we have h c ≤ α+C c < s c for all c ∈ {1, • • • n c }.

These constraints are linear and we follow the algorithm of [START_REF] Rodriguez-Yam | Ecient Gibbs sampling of truncated multivariate normal with application to constrained linear regression[END_REF] to sample the variables α and C in their conditional complete posterior distribution. This algorithm consists in i) removing the dependence within elements of C or α thanks to a simple linear transformation, ii) sampling the independent elements of the transformed variable in the convenient truncated univariate normal distribution. For this step, we use the method and the code proposed by [START_REF] Chopin | Fast simulation of truncated Gaussian distributions[END_REF]. Details of the implemented algorithm are given in Appendix C.

After testing the ability of the algorithm to estimate the parameters from simulated data correctly, three chains of 40,000 iterations including a burn-in period of 10,000 iterations were launched for each of the tested models (each of the ν values). In each case, the good mixing of the chains was checked with the Gelman-Rubin convergence diagnostic [START_REF] Brooks | Some Issues for Monitoring Convergence of Iterative Simulations[END_REF] and conrmed by visual inspection.

Results and Discussion

Choice of the smoothness parameter

The choice of the smoothness parameter of the Matérn covariance function is based on the overall model t, measured by the expected deviance, where deviance is dened as -2 times the log-likelihood [START_REF] Gelman | Bayesian Data Analysis[END_REF]. The results for each of the three considered smoothness parameters are given in Table 2. Surprisingly, we obtain the same expected deviance, meaning that for each of the three models, the t to the data is roughly the same. However, we notice an interesting dierence between models. Indeed, the eective range and the sill decrease with the increase of the smoothness parameter, whereas the nugget increases. This suggests that when the smoothness parameter is high, the spatial dependence is weaker and model predictions based on spatial dependence are less accurate. Consequently, in what follows, we decide to assign the value 0.5 to the parameter ν. This choice will be corroborated in section 6 when we compare prediction performances with the three values of this parameter. Table 3 gives the posterior mean, standard deviation and median of all estimated parameters when the exponential covariance (ν = 0.5) is considered.

Spatial repartition

The µ ct posterior variance is given by Var(C c )+Var(B t ). are explained for 81% by the covariates, the valley altitude, h c and the path exposition o c . The value of the valley altitude coecient (b=0.83 with standard deviation 0.03) indicates that its inuence is dominating. Its introduction into the modelling enables the comparison of the runout altitude irrespective of the valley altitude. The inuence of the path exposition is not so clear. Avalanches on south exposed paths generally stop higher (7.20 m), but the variability on the parameter is large (the standard deviation is 9.35 m), and we cannot fully conclude to a real eect of the exposition in the runout altitude, even if the obtained result is physically consistent. In fact, it is expected that for south exposed slopes, there is generally less snow and it is less cold than in north exposed slopes; the runout distance is therefore shorter and the runout altitude higher. Whereas the variability on the sum of sill and nugget is high, for the range parameter, nearly all of the realizations agree, suggesting that the spatial dependence does not originate from embedded-scale phenomena.

This eective range obtained here for runout altitudes is shorter than the eective range of 42 km found in the study of avalanche frequencies on the French Alps [START_REF] Lavigne | Adding expert contributions to the spatiotemporal modelling of avalanche activity under dierent climatic inuences[END_REF]. Two factors account for this dierence. First, the avalanche frequency is primarily sensitive to climate variations whereas the runout altitude is, in addition, strongly constrained by local topography (paths more or less channelized, land cover, etc.), leading to range reduction.

Second, the frequency data considered in [START_REF] Lavigne | Adding expert contributions to the spatiotemporal modelling of avalanche activity under dierent climatic inuences[END_REF] were aggregated by township, attenuating the eect of local topography even more, and leading to a range probably longer than if counts by paths had been considered.

Nevertheless the signicant value for the range we obtain here demonstrates the presence of spatial dependence in the runout altitude after crude correction of the topographic eects. This result is of crucial importance, since it has 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 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 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 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 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 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 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 q q q q qqq 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 -124 -60 -44 -34-21.3 -3.1515.9 35 69 228 m 10 km Fig. 3 Posterior mean of the random eld A of the spatial term. In black solid line, the Arve river. never been obtained before to our knowledge and proposes an alternative to the standard approaches based on topographic variables [START_REF] Lied | Empirical calculations of snow-avalanche run-out distance based on topographic parameters[END_REF][START_REF] Mcclung | Statistical and geometrical denition of snow avalanche runout[END_REF][START_REF] Sigurôsson | Transferring avalanches between paths[END_REF]. More precisely, it suggests that it is conceivable to transfer information within a 10 km neighbourhood, as an alternative to the introduction of additional variables required by standard methods for avalanche hazard mapping, which are not always easily available.

Temporal evolution

The temporal term represents only 1.3% of the variability of the spatio-temporal term µ ct , the largest part being explained by spatial variability. However, signicant year-to-year uctuations and a smooth trend can be extracted, as shown in Fig. 5. The year-to-year uctuations are essentially captured by the unstructured temporal term whose variance (δ 0 ) is estimated to 138 m 2 .

These random uctuations mainly result from the annual variations of extreme weather conditions which inuence avalanching.

At a broader time scale, the inter annual variations are modelled by the structured temporal term g. From 1925 to 1955, the mean runout altitude decreases slowly at a rate of about 0.6 m/year, a trend dicult to interpret since, even if climate reconstructions at this time are not very precise, they do not seem to show colder temperature or higher snow records [START_REF] Beniston | Climatic change at high elevation sites: an overview[END_REF][START_REF] Corona | Millennium-long summer temperature variations in the European Alps as reconstructed from tree rings[END_REF]. Furthermore, we do not have any other reconstruction of runout altitude uctuations at this early period, in France or other countries.

On the other hand, the 1955-1975 period is marked by a sudden increase of the mean runout altitude, increasing from 1200 to 1245 m in only 20 years. Finally, over the most recent years, this trend stops in accordance with winters which are somewhat harsher [START_REF] Eckert | Temporal trends in avalanche activity in the French Alps and subregions: from occurrences and runout altitudes to unsteady return periods[END_REF].

High variability in the between path variance

For most of the paths (91%), the runout altitude variance is small and does not exceed 5000 m 2 , but few paths can have a much higher variance reaching 25000 m 2 . Understanding the origin of this between path variability is a key step for in the slope at the bottom of the path seems to be the cause for the runout altitude dispersion, because trees are able to stop the less dynamic events [START_REF] Bartelt | The inuence of tree and branch fracture, overturning and debris entrainment on snow avalanche ow[END_REF]. Taking into account covariates in the variance model linked to topography or land cover would therefore be the next step to improve prediction. However, at this stage, it is impossible to obtain this information automatically, and spatial prediction will therefore be performed by employing regional distribution, i.e. without considering local specicities in runout distance dispersion around their mean.

Model choice

The BHM chosen may seem rather simple since it combines a truncated Gaussian distribution on the observation layer and an additive spatio-temporal process. However, we think that introducing more complexity in the modelling is not appropriate given the data quality and our objectives.

In the context of extreme event prediction, generalised extreme values (GEV) distributions are traditionally used to describe avalanche runout distances [START_REF] Mcclung | Extreme avalanche runout: a comparison of empirical models[END_REF][START_REF] Keylock | An alternative form for the statistical distribution of extreme avalanche runout distances[END_REF]. More generally, skewed distributions are popular for tting environmental data. For instance the gamma distributon is broadly used in hydrology [START_REF] Zalina | Selecting a probability distribution for extreme rainfall series in Malaysia[END_REF], and the lognormal or Birnbaum-Saunders distributions [START_REF] Leiva | Extreme value Birnbaum Saunders regression models applied to environmental data[END_REF][START_REF] Garcia-Papani | Birnbaum Saunders spatial modelling and diagnostics applied to agricultural engineering data[END_REF] often model air quality measurements. Our runout altitude records do not allow us to study the distribution tail, because for many avalanche paths, the runout zone is at, leading to a left censorship. The GEV distributions have thus been ruled out. This choice is conrmed by a preliminary analysis, in which the adjustment of the normal, lognormal, Weibull, and gamma truncated distributions are compared for paths on which more than 5 avalanches are recorded. The four distributions are adjusted for each path, and the AICs are averaged over the paths. The best distributions seem to be the normal (average AIC equals 236.0) and the lognormal (average AIC = 236.2). However the Weibull distribution (AIC = 237.1) performs quite well also. In contrast, the adjustement of the gamma distribution is worse (AIC = 282.8). Considering the ease of use of the normal distribution, we chose it rather than the lognormal, keeping in mind that the model is less accurate for tail distribution.

The additive spatio-temporal structure chosen (Eq. 2) enables us to directly apply the Matérn covariance and the second order random walk, two ecient techniques developped for spatial and temporal data, respectively. Moreover, thanks to the spatio-temporal separability, the computational cost is considerably reduced. However, this simple model does not allow interactions between time and space, a constraint which is not always fulllled in practice [START_REF] Fuentes | Testing for separability of spatialtemporal covariance functions[END_REF]. In order to investigate the presence of interactions between space and time, a preliminary study of the empirical spatio-temporal covariance matrices [START_REF] Cressie | Statistics for spatio-temporal data[END_REF] has been carried out (see Appendix A). It shows that assuming space-time interactions is not relevant, essentially because spatial dependance is much higher than temporal dependance. This point is conceptually corroborated from a physical point of view, as soon as we recognise that in this relatively small region the climate is constant. We expect that the same variations of the climate will have the same eects on runout altitudes for all the paths of the study region.

Prediction performances

We use the prediction method detailed in Section 3 to assess the ability of the model to predict runout altitude distributions on new undocumented sites.

We consider a cross-validation framework, in which the 51 paths with more than 30 observations are used as a test sample. They are split in 10 groups and, alternatively, each group is removed from the study. For each path of the group, draws are generated from the predictive distribution. We then compare this predictive distribution to the empirical distribution originating directly from the data, with numerical and graphical tools.

First, the Kullback-Leibler (KL) divergence of the predictive from the empirical distribution is examined. This non-symmetric measure of dierence between two probabilities can be interpreted as the loss in information by using the predictive distribution rather than the empirical one. Thus, the KL divergence is always positive, and null when the empirical and predictive distributions are identical. Here, we compute an empirical KL divergence thanks to the formula provided by Pérez-Cruz ( 2008), which has good asymptotic properties. Figure 6 maps the predicted paths, in colors, and their associated KL divergence. The best predicted paths, with the smallest KL divergences, are located along the Arve valley, a well documented region with paths of rather regular topography. For one path highlighted in Fig. 6, predictions are particularly bad: its KL divergence reaches 3.17, whereas the divergence of the other paths does not exceed 2.1. By studying this site in detail, we note that it is located at a higher altitude, with a very high observation threshold and avalanches runout at about 2200 m a.s.l (above sea level). Hence, this elevated path does not look like the standard paths of our database, which partially explains its bad prediction.

As predictions are carried out on the basis of spatial dependence, predictive paths with numerous well documented paths nearby should be better predicted than isolated paths. To check this assertion, we investigate the correlation coecient between the KL divergence and the number of observations within a given radius around the predicted path (Tab. 4). The path mentioned above is excluded from this study. When the radius is small (2 km) no correlation is observed, but when it is greater (5-15 km), there is a signicant negative correlation meaning that the loss in information by using the model rather than the empirical distribution decreases when the amount of information around a predicted path increases. This consistent result is encouraging, since it means that our approach relying on the spatial dependence between avalanche runout altitudes is an ecient means to predict the distribution on a new path.

Figure 7 shows the probability-probability plot of the predictive and empirical distributions. It illustrates the variability of the prediction ability of our method. Indeed, for some paths, the curve is close to the rst bissector, but in some other cases, it can be very far from it. Two hypotheses could be at the origin of these relatively bad predictions (i) the choice of the covariance model, or (ii) the lack of topographic variables in the model.

In order to check the rst hypothesis, we compare prediction abilities of three Matérn covariance models with smoothness parameters equal to 0.5, 1 and 1.5. The prediction ability is measured by summing the KL divergence of all of the 51 paths. It equals 59.13 when the exponential covariance model is considered, whereas it reaches 59.78 (resp. 60.33) when ν equals 1 (resp.

1.5). This small dierence in favour of the exponential model, suggests that, considering the data at our disposal, the covariance model does not seem to be the cause of the rather bad predictions, and that sticking on the exponential choice is sensible.

In order to identify the relation between local topography and predictions, four paths are investigated more deeply. Their predictive and empirical cumulative probability distributions are given in Fig. 8, while topographic details are supplied in Tab. 5. In the three cases of bad prediction (2, 3 and 4), a topographic particularity exists in the predicted path or in its surroundings (the slope is irregular with at or steep parts). On the contrary, for the case of rather good prediction (1), the path is regular and ends up far from the valley bottom. This suggests that taking into account complementary topographic features would improve the prediction for bad predicted paths. So far, only dynamical models based on a detailed topography of the path are able to deal with these irregularities, at the price of a computational burden and only if a reliable, detailed record of local activity is available. However, this leaves 4 Correlation coecient between the Kullback-Leibler divergence and the number of observation at 2, 5, 10, and 15 km around the predicted path, pvalue of the test of the null assumption of no correlation and condence interval.

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • ••• • • • • ••• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • 10 km • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • ••• • • •
• 0 0.3 0.6 0.9 1.2 1.5 1.8 3 m KL=3.17 q 1 q 2 q 3 q 4 the door open for further developments of our regional approach integrating dynamical simulations to cope for the specic topography of each path.

Conclusion and outlooks

This Table 5 Details of four predicted paths: Kullback-Leibler divergence between the predicted and empirical distribution (KL div.), prediction quality and descriptive features of the path, i.e. topographic elements which can explain the good or bad prediction. For more details, the topographic maps where paths are located can be download from ftp://avalanchesftp. grenoble.cemagref.fr/epaclpa/EPA_feuilles_carte/alpes/, the EPA map gives the name of the map to select, and the paths can be localized with their names (the rst ve numbers correspond to the postal code of the township and the three last to the path number).

strate the presence of a spatial dependence, with an eective range of about 10 km. In comparison to other quantities such as avalanche frequencies, this eective range is clearly shorter. The dierence may be explained by the factor at the origin of the spatial structure. In our case, the topography has been identied as the major control variable, whereas in the frequency studies, climate may play a larger role. The temporal term extracted from the model presents high year-to-year uctuations and a global increase of the runout altitude over recent decades.

The results roughly reect with those obtained in the northern Alps by [START_REF] Eckert | Temporal trends in avalanche activity in the French Alps and subregions: from occurrences and runout altitudes to unsteady return periods[END_REF], but a 1955-1975 period of high increase has appeared here, presumably caused by an intense construction of avalanche defense structures in Haute Savoie at this time. The moderate increase over more recent decades is coherent with the temperature increase reported in the Northern Alps [START_REF] Durand | Reanalysis of 44 Yr of Climate in the French Alps (1958-2002): Methodology, Model Validation, Climatology, and Trends for Air Temperature and Precipitation[END_REF], and may be related to global warming.

Relying on spatial dependence, a consistent framework for prediction based on kriging equations has been proposed. It is a major advance with regards to existing approaches in the avalanche literature, in which the site eect is only addressed with topographic variables. However, since this work is a rst-ever approach, certain limitations remain.

The variability between paths is high, and ground factors such as the presence of forest or the path shape have been identied as good candidates to explain it. Unfortunately, it is not possible to automatically extract these variables easily and quickly at a regional scale. Consequently, at this stage, the path-specic variances in the model are an obstacle to a conclusive predictive framework, and is partially responsible for some poor predictions. Introducing dynamical calculation in the hierarchical framework would be an attractive possibility for further improvements. Indeed, it could solve the two major diculties identied: distribution tail and local topography.

Finally, to achieve the ultimate goal, i.e. hazard mapping, spatialized return levels are necessary. The results obtained for the runout altitude in a spatial context should now be coupled with frequency estimates to supply them as function of space and, possibly, time.

Appendices

A Empirical study of the spatio-temporal covariance structure

In order to investigate the presence of spatio-temporal dependence in runout altitudes, empirical covariance matrices are estimated in the non-separable, separable and additive cases (Fig. 9). Let Y ict be the runout altitude of avalanche i on path c the year t. The empirical spatio-temporal covariance at spatial lag [h k-1 , h k ] and time lag [τ l-1 , τ l ] is given by: with N (h k , τ l ) the set of avalanche pairs (ict, i c t ) such that the distance between sites c and c is in the interval ]h k-1 , h k ] and the lag time t -t belongs to the interval ]τ l-1 , τ l ]. |N (h k , τ l )| is the number of pairs. The reference mean μY (c) is calculated by linear regression on the valley altitude, the best covariate: μY (c) = 6.3 + 7.410 -4 hc. From the empirical non-separable covariance matrix, the additive and separable covariance matrices have been estimated by least squares. In details, for each spatial and time lag, estimates of γ T l and γ S k are found such that:

ĈY (h k , τ l ) = 1 |N (h k , τ l )|
ĈY (h k , τ l ) = γ T l + γ S k + kl
for the additive covariance matrix, and:

ĈY (h k , τ l ) = γ T l γ S k + kl
for the separable covariance matrix, where kl is the error.

Figure 9 shows that the non-separable covariance matrix is similar to the separable and additive ones, essentially because the spatial dependence is much higher than the temporal dependence. This pleads further in favor of our choice of a simple model with additive eects.

B Decomposition of the predictive distribution

First we integrate the annual predictive distribution [Yc 0 t|Y] for time t and paths c 0 over the time to make prediction for the entire study period:

[Yc 0 |Y] = n t t=1 1 nt [Yc 0 t|Y].
By using the hierarchical structure of the model we can express this distribution conditionally to the process level:

(B.3)

The joint Gaussian distribution for the vector [A, Ac 0 |φ, τ 2 ] is:

A Ac 0 ∼ N 0 0 , τ 2 Σ .. φ Σ .0 φ Σ 0. φ Σ 00 φ ,
where Σ .. φ and Σ 00 φ are the variance-covariance matrices of A and Ac 0 , respectively, and where Σ .0 φ and Σ 0. φ are covariance matrices between elements of A and Ac 0 . Then the formulas for conditional Gaussian distributions give:

E(Ac 0 |A, φ, τ 2 ) = Σ 0. φ Σ ..-1 φ A V (Ac 0 |A, φ, τ 2 ) = Σ 00 φ -Σ 0. φ Σ ..

C Algorithm details

We detail the algorithm to draw C under the constraint nc c=1 Cc = 0, the same method is used to draw B under the constraint n t t=1 Bt = 0.

Let X C the n × nc matrix such that X C ij = 1 if the avalanche i occurs on path j, 0 either, D the diagonal variance-covariance matrix of Z, D = diag(σ 2 c i ) i=1,n , with c i the path label for avalanche i, and 1c the nc length column vector of ones.

The complete posterior distribution of C is given by,

C|Z, B, α, D, E[C], ρ 2 ∼ N (m, Σ) with m = (X C D -1 X C + 1 ρ 2 Ic) -1 (X C D -1 (Z -α -X B B) + 1 ρ 2 E(C)) Σ = (X C D -1 X C + 1 ρ 2 Ic) -1

  Fig. 1 Localization of the EPA avalanche paths in Haute-Savoie, French Alps. The paths are gured out by a contour line, the gray scale in the background represents the altitude.

  Fig. 2 log-likelihood for a sample of 100 independent observations of a standard Gaussian distribution right-truncated at the value 0.1.

  and add the variability corresponding to the several sources of uncertainties in order to obtain the predictive distribution [Y c0 |Y], where Y = (Y 111 , • • • , Y ncntnn c n t ) . Considering the hierarchical structure of our model, it is straightforward to simulate predictive runout altitudes for new paths from the MCMC outputs. Details of the method are given in

Figure 3

 3 Figure 3 displays the purely spatially structured random term A c . We distinguish clearly a spatial pattern. For instance in the Arve valley (highlighted with a black solid line) avalanches seem to stop lower on the west side than on the east side. The estimates of the variance parameters conrm the presence of spatial dependence. The sill is estimated to 4669 m 2 , whereas the nugget is only 2479 m 2 and the eective range is about 9.6 km. The uncertainty regard- ing the three parameters (sill, range and nugget) are shown in Fig. 4, which represents the posterior mean semivariogram and two hundred realizations.

Fig. 4

 4 Fig. 4 Red: mean posterior variogram, obtained with sill, range and nugget posterior mean; black: posterior variogram realizations. (m 2 )

Fig. 5

 5 Fig.5Posterior mean of the temporal term (solid line), posterior mean of the smooth term (dotted line), and 95% credible interval of the smooth term (grey area).

  10 -4 [-0.66; -0.22] Table

Fig. 6

 6 Fig. 6 Kullback-Leibler divergence for the 51 predictions paths (colored). The paths in black are in the learning sample.

Fig. 7

 7 Fig. 7 Probability-probability plot of the predictive and the empirical distributions. Each curve represents one of the 51 paths where the runout altitude is predicted (validation sample). Number labels identify the four paths described in Tab. 5.

  scheme.Important quantities for the avalanche hazard assessment eld have been inferred. The most relevant are the spatial covariance parameters, which demon-

Fig. 8

 8 Fig. 8 Predictive (red line) and empirical (green line) cumulative probability distributions for the four paths detailed in Tab. 5.

Fig. 9

 9 Fig.9Empirical non-separable (left), separable (center) and additive (right) spatiotemporal covariance matrices. The mean is only spatially indexed, and is computed from the valley altitude.

(

  ict,i c t )∈N (h k ,τ l ) (Y ict -μY (c))(Y i c t -μY (c ))

[

  Yc 0 t|Y] = [Yc 0 t|Cc 0 , Bt, α, σ 2 c 0 ][Cc 0 , Bt, α, σ 2 c 0 |Y]d(Cc 0 , Bt, α, σ 2 c 0 ). (B.1) Let θ 1 = (Bt, α, σ 2 c 0 ), we decompose the distribution [Cc 0 , θ 1 |Y] in the same way to make the spatial term Ac 0 appear,[Cc 0 , θ 1 |Y] = [Cc 0 , θ 1 |Ac 0 , a, b, o S , ρ 2 ][Ac 0 , a, b, o S , ρ 2 |Y ]d(Ac 0 , a, b, o S , ρ 2 ).

  express [Ac 0 , θ 2 |Y], with θ 2 = (a, b, o S , ρ 2 ) , conditionally to A rather than Y: [Ac 0 , θ 2 |Y] = [Ac 0 , θ 2 |A, φ, τ 2 ][A, φ, τ 2 |Y]d(A, φ, τ 2 ).

  . (B.2) and (B.3) within Eq. (B.1), and using the Bayes' rule we obtain[Yc 0 ,t|Y] = [Yc 0 t|Cc 0 , θ 1 ] [Cc 0 |Ac 0 , θ 2 ][Ac 0 |θ 3 ] [θ 1 , θ 2 , θ 3 |Y]d(θ 1 , θ 2 , θ 3 ), (B.4)with θ 3 the vector (A, φ, τ 2 ) .

  Table 2 Comparison of the parameters posterior mean for the three Matérn variogram models.

	ν	0.5		1		1.5
	φ	3209 m	1652 m	1262 m
	Eective range	9460 m	6591 m	5982 m
	sill τ 2	4669 m	2	3968 m	2	3738 m	2
	nugget ρ 2	2478 m	2	2842 m	2	3019 m	2
	Expected Deviance	4.718.10 4 4.718.10 4 4.718.10 4
						The ratio	Var(Cc) Var(Cc)+Var(Bt)
	equals 98.7% indicating that spatio-temporal variations of runout-altitudes are
	mainly due to spatial variations. The large variations of the spatial component

Table 3

 3 Parameter estimates: posterior mean, posterior standard deviation and posterior median (standard deviation 162.6 m) and corresponds to quite a well documented path with 10 observations, but for which two avalanches stop on a gentler slope before the valley. Thus, the specic topography of the site (two distinct runout zones) is responsible for the high variance observed. In other cases of large variance, the paths are still well documented, and the presence of forest

	parameters	mean				standard deviation	median
	φ	3209.40 m	2102.98 m	2719.85 m
	τ 2	4669.10 m	2	1230.73 m	2	4540.29 m	2
	ρ2	2478.59 m	2	651.17 m	2	2485.95 m	2
	δ 0	137.69 m	2	33.47 m	2	133.76 m	2
	δ 1	0.85 m	2			1.06 m	2	0.49 m	2
	α	1237.86 m	1.22 m	1237.84 m
	a	-935.56 m	31.03 m	-936.02 m
	b	0.83				0.03		0.83
	o S	7.20 m			9.35 m	7.13 m
	making predictions on new paths. Here, the largest variance reaches 26425 m	2

Acknowledgements The authors thank the ANR research program MOPERA (Modélisation probabiliste pour l'Étude du Risque d'Avalanche -www.avalanches.fr/mopera-projet/) for funding this work.

At a second step we write the joint distribution of C|Z, B, α, D, E[C], ρ 2 , noted C|γ and its constraint 1cC = 0,

The conditional distribution is then:

Only the nc -1 rst components of C are sampled, Cn c is given by Cn c = -nc-1 c=1 Cc.

Then we constrain the mean runout altitude by path, α + Cc. We have hc < α + Cc < sc for all c ∈ {1, • • • nc}. In a Gibbs sampler, C and α are drawn successively in their complete conditional distributions. We describe the method to take into account the constraint for C, α is sampled in a similar way. The nc -1 rst components of C must fullll the constraints:

that we can write matrix-wise VC < v with V the k × nc -1 matrix of constraints. Here k = 2nc, since there are two constraints for each line of the system Eq. (C.2), one bears on sc, the other one on hc. The posterior distribution of C is thus a multivariate truncated distribution whose mean and variance for the complete version are given by Eq. (C.1) and whose support S = {x ∈ R nc-1 : Vx < v}. We follow the two-step algorithm of Rodriguez-Yam et al ( 2004). In the rst step, we remove the dependence between the components of C by considering the vector w = LC, with L such that LΣ 0 L = I nc-1 and I nc-1 the identity matrix. It has the following truncated distribution:

with support: S = {w ∈ R nc-1 : VL -1 w < v}.

In the second step, each component of w is drawn successively conditionally to the others in a Gibbs sampling. In order to avoid the computation of the entire vector of constraints VL -1 w at each draw, we may notice that for each constraint i ∈ {1, • • • k}, and each component of

The constraint i for the j 0 component of w is written

. To update w j 0 , one only has to compute u j 0 i for each constraint. This step demands k(nc -1) elementary operations, instead of k(nc -1) 2 operations when VL -1 is computed naively.

Finally we implement the following algorithm for sampling the nc -1 components of C from their truncated Gaussian distributions:

Compute the vector u = VC; Initiate the matrix L, its inverse L -1 , and the vector w = LC;

For each

, and u -j = u -u j w j ; Find the interval [a j , b j ] in which w j is drawn to satisfy the k constraints; Sample w j in the normal truncated distribution N ((Lm 0 ) j , 1) with support [a j , b j ].

Here, we use the method and the code proposed by [START_REF] Chopin | Fast simulation of truncated Gaussian distributions[END_REF]; Update u as u = u -j + u j w j ; Set C as C = L -1 w.