Detecting and modeling critical dependence structures between random inputs of computer models

: Uncertain information on input parameters of computer models is usually modeled by considering these parameters as random, and described by marginal distributions and a dependence structure of these variables. In numerous real-world applications, while information is mainly provided by marginal distributions, typically from samples, little is really known on the dependence structure itself. Faced with this problem of incomplete or missing information, risk studies that make use of these computer models are often conducted by considering independence of input variables, at the risk of including irrelevant situations. This approach is especially used when reliability functions are considered as black-box models. Such analyses remain weakened in absence of in-depth model exploration, at the possible price of a strong risk misestimation. Considering the frequent case where the reliability output is a quantile, this article provides a methodology to improve risk assessment, by exploring a set of pessimistic dependencies using a copula-based strategy. In dimension greater than two, a greedy algorithm is provided to build input regular vine copulas reaching a minimum quantile to which a reliability admissible limit value can be compared, by selecting pairwise components of sensitive influence on the result. The strategy is tested over toy models and a real industrial case-study. The results highlight that current approaches can provide non-conservative results.


Introduction
Many industrial companies, like energy producers or vehicle and aircraft manufacturers, have to ensure a high level of safety for their facilities or products. In each case, the structural reliability of certain so-called critical components plays an essential role in overall safety. For reasons related to the fact that these critical components are highly reliable, and that real robustness tests can be very expensive or even hardly feasible, structural reliability studies generally use simulation tools [18,43]. The physical phenomenon of interest being reproduced by a numerical model η (roughly speaking, a computer model), such studies are based on the calculation of a reliability indicator based on the comparison of y = η(x) and a safety margin, where x corresponds to a set of input parameters. In the framework of this article, such models are considered as black-boxes and can be explored only by simulation means.
While the problems of checking the validity of η and selecting inputs x ∈ χ ⊆ R d are addressed by an increasing methodological corpus [3,15], a perennial issue is the modeling of x. Di ering from the specication of η itself, this input vector is known with uncertainty, either because the number of experiments to estimate is limited, or because some inputs re ect intrinsically variable phenomena [59]. In most cases, these epistemic and aleatory uncertainties are jointly modeled by probability distributions [34]. Consecutively, the reliability indicator is often de ned as the probability that y be lower than a threshold (failure probability), or a limit quantile for y. This article focuses on this last indicator, which provides an upper or lower bound of the mean e ect of the output variable uncertainty.
Therefore the modeling of x stands on the assessment of a joint probability distribution with support χ, divided between marginal and dependencies features. Though information on each dimension of x can often be accessible experimentally or using physical or expert knowledge [5], the dependence structure between the component of x remains generally unknown. Typically, statistical data are only available per dimension, but not available for two or more dimensions simultaneously. For this reason, most of robustness studies are conducted by sampling within independent marginal distributions. Doing so, reliability engineers try to capture input situations that minimize the reliability indicator.
However, the assumption of independence between inputs has been severely criticized since the works by [30] and [62], who showed that output failure probabilities of industrial systems can signi cantly vary and be underestimated if the input dependencies are neglected. More generally, [60,61] showed that tail dependencies between inputs can have major expected e ects on the uncertainty analysis results. A noticeable contribution to this analysis was made by [64], who highlighted, using vine copula models, the strong impact of input dependence structures on reliability indicators computed by well-known methods of uncertainty quanti cation.
In this paper, we consider a general probabilistic and statistical approach at modelling dependence with a supervised target. In particular, we study the speci c problem of structural reliability as an application of such method. In decision-making problems, [40] proposed a general de nition of the worst penalizing distribution as the minimizer of an excepted cost among a set of possible distributions. More recently, [2] extended this approach to account for incomplete dependence information. These theoretical works, that propose selection rules over the in nite set of all possible joint distributions, remain hard to apply in practice. Recent applied works made use of copulas [50] to model dependencies between stochastic inputs [60,61], following other researchers confronted to similar problems in various elds: nance [13], structural safety [29], environmental sciences [57] or medicine [4]. These studies mainly consider bivariate copulas, which makes theses analysis e ective only when two random variables are correlated. Cases where a greater number of variables is involved were explored by [37] then [64], who used vine copulas to approach complex multidimensional correlation problems in structural reliability, and by [63] in a speci c context of regression by polynomial chaos expansion. A vine copula is a graphical representation of the pair-copula construction (PCC), proposed by [39], which de nes a multidimensional dependence structure using conditional bivariate copulas. Various class of vines exist (see [17] for a review), and among them the regular vines (R-vines) introduced by [6,7] are known for their appealing computational properties, while inference on PCC is usually demanding [20,33].
R-vine parametric copulas seem promising to improve the search for a worst-case dependence between stochastic inputs, while keeping the bene ts of a small number of parameters, as favoring inference and conducting simple sensitivity analyses a posteriori. To our knowledge, however, no practical methodology has been yet proposed to this end for which the notion of worst case is de ned by the minimization of an output quantile. This is the subject of this article. More precisely, the aim of this research is to determine a parametric copula over x, close to the worst case dependence structure, which is associated to a minimum value of the quantile of the distribution of y. Given a vine structure de ned by a parameter vector, the optimization problem involves to conduct empirical quantile estimations for each value of this vector in a nite set of interest (chosen as a grid). The proposed methodology stands on an encompassing greedy algorithm exploring copula structures, which integrates several sub-algorithms of increasing complexity and is based on some simplifying assumptions. These algorithms are made available in the Python library dep-impact [8].
The article is therefore organized as follows. Section 2 introduces the framework and studies the consistency of a statistical estimation of the minimum quantile, given an input copula family and a growing sequence of grids. The wider problem of selecting copulas in high-dimensional settings using a sequential algorithm of quantile minimization is considered in Section 3.
While the choice of R-vines is defended, similarly to [64] a sparsity hypothesis is made to diminish the computational burden, according to which only a limited number of pairwise dependencies is in uent on the result. A greedy algorithm is proposed to carry out the complete procedure of optimization and modeling. This heuristic is tested in Section 4 over toy examples, using simulation, and a real industrial case-study. The results highlight that worst-case scenarios produced by this algorithm are often bivariate copulas reaching the Fréchet-Hoe ding bounds [27,35] (describing perfect dependence between variables), as it could be expected in monotonic frameworks, but that other nontrivial copulas can be exhibited in alternative situations. Results and avenues for future research are extensively discussed in the last section of this article. We also refer to Appendix A and B for supplementary material on consistency proofs, on R-vine copulas and on R-vine iterative construction.

Minimization of the quantile of the output distribution
This section introduces a general framework for the calculation of the minimum quantile of the output distribution of a computational model, when the input distribution can be taken from a large family of distributions, each one corresponding to a particular choice of dependencies between the input variables.

. A general framework for the computation of the minimum quantile
To be general, let us consider a computer code which takes a vector x ∈ χ ⊆ R d as an input and produces a real quantity y in output. This code is represented by a deterministic function η : R d → R such that η(x) = y. The sets R and R d are endowed with their Borel sigma algebras and we assume that η is measurable. The general expression of the function η is unknown but for some vector x ∈ R d it is assumed that the quantity η(x) can always be computed. In particular, the derivatives of η, when they exist, are never assumed to be known.
Let P , . . . , P d be a xed family of d distributions, all supported on R. We introduce the set D(P , . . . , P d ) of all multivariate distributions P on R d such that the marginal distributions of P are all equal to the (P j ) j= ...d .
Henceforth, we use the shorter notation D for D(P , . . . , P d ).
For some P ∈ D, let G be the cumulative distribution function of the model output. In other terms dG is the push-forward measure of P by η. For α ∈ ( , ), let G − be the α-quantile of the output distribution: For the rest of this document, we denote as output quantile the α-quantile of the output distribution.
In many real situations, the function η corresponds to a known physical phenomenon. The input variables x of the model are subject to uncertainties and are quanti ed by the distribution P. The propagation of these uncertainties leads to the calculation of the output quantile. Due to the di culties to gather information, it is common to have this distribution incompletely de ned and only known through its marginal distributions. Therefore, the set D corresponds to all the possible distributions that are only known through their marginal distributions (P j ) j= ...d . Detecting very unfavorable correlations between inputs is similar to detecting correlations leading to low quantile values. We de ne as the worst quantile, the minimum value of the quantile by considering all the possible input distributions P ∈ D. This conservative approach consists in minimizing G − (α) over the family D such as Since the function η has no closed form in general, it is not possible to give a simple expression of G − (α) in function of the distribution P, and consequently the minimum G − * (α) does not have a simple expression Figure 1: Illustration of the link between the dependence parameter θ and the quantile function G − θ . The joint CDF F θ is obtained using (3) from a copula C θ and marginal CDF's (F j ) d j= . The push-forward of F θ through the model η leads to the CDF G θ and quantile function G − θ of the output distribution.
too. In this paper we propose to study a simpler problem than (2), by minimizing G − (α) over a subset of D. This subset is a family of distributions (P θ ) θ∈Θ associated to a parametric family of copula (C θ ) θ∈Θ , where Θ is a compact set of R p and p is the number of copula parameters.

. Copula-based approach
We introduce the real-values random vector X = (X , . . . , X d ) ∈ R d associated to the distribution P θ . Each component X j , for j = , . . . , d, is a real-value random variable with distribution P j . A copula describes the dependence structure between a group of random variables. Formally, a copula is a multidimensional continuous cumulative distribution function (CDF) linking the margins of X to its joint distribution. Sklar's Theorem [58] states that every joint distribution F θ associated to the measure P θ can be written as with some appropriate d-dimensional copula C θ with parameter θ ∈ Θ and the marginal CDF's F j (x j ) = P[X j ≤ x j ]. If all marginal distributions are continuous functions, then there exists a unique copula satisfying where u j = F j (x j ). For F θ absolutely continuous with strictly increasing marginal distributions, one can derive (3) to obtain the joint density of X: where c θ denotes the copula density function of C θ and f j (x j ) are the marginal densities of X. Numerous parametric copula families are available and are based on di erent dependence structures. Most of these families have bidimensional dependencies, but some can be extended to higher dimensions. However, these extensions have a lack of exibility and cannot describe all types of dependencies [50]. To overcome these di culties, tools like vine copulas [38] (described in Section 3) combine bivariate copulas, from di erent families, to create a multidimensional copula. Let G θ and G − θ be respectively the CDF and quantile function of the push-forward distribution of P θ by η (see Figure 1). For a given parametric family of copula (C θ ) θ∈Θ and a given α ∈ ( , ), the minimum output quantile for a given copula is de ned by and if it exists, we consider a minimum We call this quantity the minimum quantile parameter or worst dependence structure. Note that there is no reason for G − θ (α) to be a convex function of θ. The use of gradient descent algorithms is thus not straightforward in this context. Moreover, the gradient of θ → G − θ is unknown and only zeroorder optimization methods can be applied to solve (6). For this reason, in the following of this section, we analyze the basic approach which consists in estimating θ * C by approximating Θ with a nite regular grid Θ N of cardinality N. Therefore, for a given parametric copula (C θ ) θ∈Θ and a given α ∈ ( , ), we restrict the problem (6) to . Estimation with a grid search strategy In the restricted problem (7), the greater N, the closer θ * N to the minimum θ * of Θ; obviously the convergence rate should depend on the regularity of the function η and on the regularity of the quantile function θ → G − θ (α). Because η has no closed form, the quantile function G − θ (α) has no explicit expression. The minimizer θ * N can be estimated by coupling the simulation of independent and identically distributed (i.i.d) data (Y , . . . , Yn), de ned as realizations of the model output random variable Y := η(X) with distribution dG θ , with a minimization of the empirical quantile over Θ N . Note that this problem does not deal with tting vines on a set of observations: here we generate data according to various vine distributions and we estimate quantiles under these distributions.
For θ taking a value over the grid Θ N , the empirical CDF of Y is de ned for any y ∈ R by The corresponding empirical quantile function G − θ (α) is de ned as in (1) by replacing G with its empirical estimate. For a given α, the worst quantile on the xed grid Θ N is given by and can be estimated by replacing the quantile function with its empirical function: Finally the estimation of the minimum quantile parameter over the grid Θ N is denoted by The construction of the grid Θ N can be di cult because Θ can be unbounded (e.g. Θ = [ , ∞] for a Gumbel copula). To tackle this issue, we chose to construct Θ N among a normalized space using a concordance measure, which is bounded in [− , ] and does not rely on the marginal distributions. We chose the commonly used Kendall rank correlation coe cient (or Kendall's tau) [41] as a concordance measure to create this transitory space. This non-linear coe cient τ ∈ [− , ] is related to the copula function as follows: For many copula families, this relation is much more explicit (see for instance [28]). Therefore, the nite grid is created among [− , ] p and each element of this grid is converted to the copula parameter θ. Moreover, the use of concordance measures gives a normalized expression of the strength of dependencies for all pairs of variables, independently of the used copula families. The consistency of estimators (9) and (10) is studied in next section, under general regularity and geometric assumptions on η and the functional θ → P θ .

. Consistency of worst quantile-related estimators
In this section, we give consistency results of the estimators min θ∈Θ N G − θ (α) and θ N , for a growing sequence of grids on the domain Θ. For easier reading, we skip some de nitions needed for our assumptions. Section A in Appendix provides a more complete presentation, including the formal de nition of the modulus of increase of the quantile function.
Let α be a xed value in ( , ). To approximate Θ, we consider a sequence of nite discrete grids (Θ N ) N≥ on Θ where N is the cardinal of Θ N and such that sup θ∈Θ, θ ∈Θ N θ − θ → as N tends to in nity.
We rst introduce technical hypotheses required for the consistency result which are commented further in the text.
Assumption A. For all θ ∈ Θ, the distribution P θ admits a density f θ for the Lebesgue measure and the copula C θ admits a density c θ for the Lebesgue measure on [ , ] d such that is a continuous function.

Assumption B.
For all θ ∈ Θ, G θ is a continuous function.
Assumption C. For all θ ∈ Θ, G θ is strictly increasing and the modulus of increase of G θ at G − θ (α) is lower bounded by a positive function ϵ Θ .
Let (Nn) n≥ be a sequence of integers such that Nn n β for some β > . For every n ≥ we consider the grid Θ Nn and for every θ ∈ Θ Nn we compute the empirical quantile G − θ (α) from a sample of n i.i.d variables Y , . . . , Yn with Y i = η(X i ), where the X i s are i.i.d. random vectors with distribution P θ . We then introduce the extremum estimatorθ := θ Nn .
Moreover, if Assumption D is also satis ed, then for all h > we have It would be possible to provide rates of convergence for this extremum quantile and for θ * at the price of more technical proofs, by considering also the dimension metric of the domain Θ and the modulus of increase of the function θ → G θ (α) (see for instance the proofs of Theorems 1 and 2 in [12] for an illustration of such computations). It would be also possible to derive similar results for alternative extremum quantities. One rst example, useful in many applications, would be to estimate some "limit probability" by determining an extremum inf θ∈Θ G θ (y) of the CDF for a xed y.
This consistency result could also be extended for regular functional of G θ or G − θ , such that inf θ∈Θ y≥y G θ (y)dy or inf Output CDF and PDF Input Probabilistic Space for some xed values y and α . Extending our results for such quantities is possible essentially because the Dvoretsky-Kiefer-Wolfowitz (DKW) inequality [21], used in the proof, gives an uniform control on the estimation of the CDF and the quantile function. We now discuss the three rst assumptions and provide some geometric and probabilistic interpretations of them. Assumption A requires some regularity of the input distribution with respect to θ. This is indeed necessary to locate the minimum of the quantile. Assumption B and C ensure that the output quantile function G − θ has a regular behavior in a neighborhood of the computed quantile G − θ (α). Assumption B ensures that the output distribution dG θ has no Dirac masses whereas Assumption C ensures that there is no area of null mass inside the domain of dG θ . Figure 2 illustrates Assumption B with a possible con guration of the input distribution. For θ ∈ Θ an δ > , we consider a small neighborhood [G − θ (α)−δ, G − θ (α)+δ] of G − θ (α), and the pre-image of this neighborhood. The two right gures are the CDF G θ (top) and PDF g θ (bottom) of the output variable Y for a given θ. The gure at the left hand represents the contours of the pre-image in the input space. The red plain line is the level set η − (G − θ (α)) and the dot blue line is the perturbed level set η − (G − θ (α) ± δ)). The blue area in the right gure corresponds to [G − θ (α) − δ, G − θ (α) + δ] and the pre-image of this neighborhood is the blue area in the left gure. Assumption B requires that the mass of the blue domain is lower bounded by a positive function ε Θ (δ) that does not depend on θ.
It is possible to give su cient conditions on the input distribution F θ and on the geometry of η to obtain Assumptions B and C. Using the de nition of the modulus of continuity from Equation (A.1) in Appendix, it comes where H d− is the d − dimensional Hausdor measure (see for instance Chapter 2 in [26]). If the copula and the code are such that there exists a constant I such that for any θ ∈ Θ and any u in the support of dG θ Note that ∇η ∞ < ∞ since η is assumed to be Lipschitz. We have proved that Assumption C is satis ed in this context. Finally, by rewriting again the co-area formula for G θ (y), we nd that Assumption B is satis ed as soon as the set of stationary points ( ∇η(x) = ) of all level set η − {u} has null mass for the Hausdor measure.
In conclusion, we see that for smooth copulas, Assumptions C and B mainly depend on the regularity of the function η, by requiring on one side that η does not oscillate to much and on the other side that the set of stationary points does not have a positive mass on the level sets of η.

Quantile minimization and choice of penalized correlation structure
A preliminary study of the in uence of the dependence structure, speci c to quantile minimization, is conducted in Appendix C. It also illustrates that the worst quantile is not always reached for a trivial copula (i.e., reaching the Fréchet-Hoe ding bounds).
The current section rst provides a rationale for choosing the so-called R-vine structure as a preferential copula structure for modeling the variety of correlations between inputs. Then, the search for a minimum quantile is presented in two times. Subsection 3.2 proposes an exhaustive grid-search algorithm for estimating this quantile when the R-vine copula structure is xed with a given pair-copula families and indexed by the parameter vector θ. Subsection 3.3 extends this rigid framework by permitting the search of particular subcopula pairwise structures, such that the minimization be more signi cant. In each situation, examples are provided to demonstrate the feasibility of the approach.

. A rationale for R-vine copula structures
Representing dependence structures in high dimensional settings is a challenging problem. For the following de nition, we simplify the expressions by omitting the use of θ: f = f θ , F = F θ and c = c θ . By recursive conditioning, the joint density can be written as a product of conditioning distributions such as For clarity reason, we now simplify the expression with f | , = f | , (x |x , x ) and so on for other orders. From (4), the conditioning densities of (14) can be rewritten as products of conditioning copula and marginal densities. For example, in a case of three variables and using (14), one possible decomposition of the the joint density can be written as Using (4), the reformulation of f | , leads to where . By developing f | , in the same way, we nd that Thus, by replacing the expression of f | , in (15) and doing the same procedure for f | , the joint density can be written as This nal representation of the joint density based on pair-copulas has been developed in [39] and is called the pair-copula construction (PCC). The resulting copula represented by the product of conditional copulas in (18) o ers a very exible way to construct high-dimensional copulas. However, it is not unique; indeed, (14) has numerous decomposition forms and it increases with the dimension.
To describe all such possible constructions in an e cient way, [6,7] introduced the vine models. This graphical tool, based on a sequence of trees, gives a speci c way to decompose the multivariate probability distribution. Basically, a vine model is de ned by • a structure of trees which can be represented by a matrix [48], • a copula family for each pair of the structure, • a parameter for each pair-copula.
A R-vine is the general construction of a vine model, but particular cases exists such as the D-vines and Cvines, described in Appendix B. Vine models were deeply studied in terms of density estimation and model selection using maximum likelihood [1], sequential estimation [20,42], truncation [11] and Bayesian techniques [31]. Inspired by recent studies like [64], their popularity and well-known exibility led us to use R-vines in this article, despite the fact that in our context we are looking for a conservative form and not to select the most appropriate form with given data, in absence of correlated observations.

. Estimating a minimum quantile from a given R-vine . . Grid-search algorithm
Let Ω = {(i, j) : ≤ i, j ≤ d} be the set of all the possible pairs of X, in a d-dimensional problem. The number of pairs p is associated to the size of Ω such as p = |Ω| = d = d(d− )/ . We de ne V as the vine structure and we consider xed copula families for each pair. In this article, we only consider single parameter pair-copulas, such that the parameter θ is a p-dimensional parameter vector with a de nition space Θ := (i,j)∈Ω Θ i,j where Θ i,j is the parameter space for the pair-copula of the pair (i, j). However, the methodology can easily be extended to multi-parameter pair-copulas. Note that a pair-copula can be conditioned to other variables, depending on its position in the vine structure V. Thus, the input distribution dF θ (V) is de ned by the vine structure V, the copula families and the parameter θ. Also note that the copula parameter θ is associated to the R-vine structure V (i.e., θ = θ V ), see Section 3.2.2. For the sake of clarity, we simplify the notation to θ only.
The most direct approach to estimate the minimum quantile is the Exhaustive Grid-Search algorithm, described by Algorithm 1.
Algorithm 1: Exhaustive Grid-Search algorithm to minimize the output quantile.
For a given vine structure V, copula families, a grid Θ N and a sample size n, three steps are needed for each θ ∈ Θ N . The rst step simulates an input sample {X i } n i= according to the distribution dF θ (V) for a given sample size n. The second evaluates the sample through the model η. The third estimates the output quantile from the resulting sample {Y i = η(X i )} n i= . The minimum quantile is took among the results of each loop.

. . Influence of the vine structure
Using R-vines, the dependence parameter θ is associated to the vine structure V. Due to the hierarchy of the vine structure, some pair-copulas are conditioned to other variables and thus for their parameters. As an illustration, let us consider two vine structures with the two following copula densities, with the same simpli ed expressions as for (18): The di erence between these densities is the conditioning of some pairs, the dependence parameters of theses vines are Applying the same grid for these two vines may give di erent results due to the conditioning order from the vine structure. For example, if the pair X -X is very in uential on minimizing the output quantile, it would be more di cult to nd a minimum with V than V due to the conditioning of the pair with X and X in V . However, if the grid is thin enough, the minimum from these two vines should be equivalent.
To counter this di culty, one possible option consists in randomly permuting the indexes of the variables and repeating the algorithm several times to visit di erent vines structures.

. . Computational cost
For one given R-vine structure and one xed copula family at each pair, the overall cost of the method is equal to nN. However, as explained in § 2.4, the nite grid Θ N , should be thin enough to reasonably explore Θ. Therefore, N should increase with the number of dimensions d and more speci cally with the number of pairs p = d . A natural form for N would be to write it as N = γ p , where γ ∈ R + . Thus, the overall cost of the exhaustive grid-search would be equal to nγ ( d ) . The cost is in O(γ d ) which makes the method hardly scalable when the dimension d increases. Ideas about alternative optimizations methods are discussed in Section 5.
. Iterative search for a penalizing R-vine structure: a greedy heuristic based on pairwise copula

. . Going further in quantile minimization
With Algorithm 1, the previous subsection proposes an exhaustive grid-search strategy to determine a R-vine copula Cθ such that the associated output quantile G − θ (α) be the smallest (and also the most conservative in a structural reliability context). This approach remains however limited in practice since Cθ for xed paircopula families (e.g., Archimedean or max-stable copulas) and V which is a member of the set F d of all the possible d−dimensional R-vine structure. Intuitively, a more reliable approach to quantile minimization should be based on mixing this estimation method with a selection among all members of the nite set F d , as well for the copula families. It is indeed likely that searching within an associative class of copulas like Archimedean ones, allowing modeling dependence in arbitrarily high dimensions, be a too rigid choice for estimating the minimum G − θ (α).
A minimum quantile can probably be found using a R-vine structure de ned by conditional pairwise subcopulas (according to (18)) that are not part of the same rigid structure. However, a brute force exploration of F d would be conducted at an exponential cost increasing with d [49]. If we also consider the large computational cost of an exhaustive grid-search for a large number of dependent variables (as explained in § 3.2), this approach is not feasible in practice for high dimensions.
For this reason, it is proposed to extend Algorithm 1 by a greedy heuristic that dynamically selects the most in uential correlations between variables while limiting the search to pairwise correlations. Doing so, minimizing the output quantile can be conducted in a reasonable computational time. Therefore the selected d−dimensional vine structure would be lled with independent pair-copulas except for the pairs that are in uential on the minimization.
This working limitation, interpreted as a sparsity constraint, is based on the following assumption: it is hypothesized that only few pairs of variables have real in uences on the minimization. It is close in spirit to the main assumption of global sensitivity analysis applied to computer models, according to which only a limited number of random variables has a major impact on the output [36,56].

. . General principle
The method basically relies on an iterative algorithm exploring pairwise correlations between the uniform random variables U j = F − j (X j ) and progressively building a non-trivial R-vine structure, adding one pair of variable to the structure at each iteration. Starting at step k = from the simple independent copula the algorithm nally stops at a given step k = K while proposing a new copula C θ (K) associated to a R-vine structure V K mostly composed of independent pair-copulas.
At each iteration k, we denote by Ω k the selected pairs which are considered non-trivial (nonindependent) due to their in uence on the quantile minimization. Let Ω −k = Ω\Ω k be the candidate pairs, which were not the remaining pairs, which in uence on the minimization is still to be tested and are still considered independent. We also consider B as a set of candidate copula families. The pseudo-code of Algorithm 2 shows in detail how this iterative exploration and building is conducted. More algorithms in Appendix B.2 described how to construct a vine structure with a given list of indexed pairs of variable.

. . Example
Consider the four-dimensional (d = ) situation such as X = (X , . . . , X ) where, for to the sake of simplicity, all marginal distributions of X are assumed to be uniform on [ , ]. We consider a simple additive model described by For an additive model and uniform margins, the output quantile is monotonic with the dependence parameters (see Subsection 2 in Appendix C) which locates the minimum quantile at the edge of Θ. Thus, Step 1.b. of Algorithm 2 is simpli ed by considering only Fréchet-Hoe ding copulas in the exploration.
In this illustration we consider α = . and we select n = , large enough in order to have a great quantile estimation and the algorithm stops at K = . Figure 3 shows, for each iteration k, the p − k vine structures that have been created by the algorithm. The red nodes and edges are the candidate pairs (i, j) ∈ Ω −k and the blue nodes and edges are the selected pairs Ω k . At iteration k = , the selected pair is ( , ) with an estimated minimum quantile of − . . At iteration k = , the second selected pair is ( , ) with an estimated minimum quantile of − . . At iteration k = , the third selected pair is ( , ) with an estimated minimum quantile of − . . ( 3. Check the stopping condition; Extend the list of selected pairs: Ω k = Ω k ∪ (i, j) (k) and families: if k < K and computational budget not reached then 21 New iteration: k = k + ; We observe that X appears in all the selected pairs. This is not surprising since X is the most in uential variable with the largest coe cient in (21). The algorithm considers D-vines by default, but this is important for the rst iterations since most of the pairs are independent. When it is possible, the algorithm creates a vine such as the selected pairs and the candidate pair are in the rst trees. For example, the fourth vine at iteration k = with the candidate pair ( , ) shows a R-vine structure that respects the ranking of the listed pairs. However, the third vine at iteration k = for the candidate pair ( , ) along with the selected pairs {( , ), ( , ) could respect the ranking and set all the pairs in the rst tree altogether. Thus, using Algorithm 3 in Appendix, a valid vine structure is determined by placing the candidate pair ( , ) in the next tree.

. . Computational cost
The number of model evaluations is in uenced by several characteristics from the probabilistic model and from the algorithm. Let |B| be the number of family candidates. The total number of runs is The sum corresponds to the necessary iterations to determine the in uential pairs. The maximum possible cost is if all the pairs are equivalently in uential (i.e., K = p = d(d − )/ ), which would be extremely high. The term N k is the cost from the grid-search quantile minimization at step 2. of the algorithm. The greater N k is and the better the exploration of Θ k ∪ Θ i,j . Because the dimension of Θ k increases at each iteration k, it is normal that N k should also increases with k (e.g. N k = γ βk , where γ and β are constants). Also, the greater n is and the better the quantile estimations. The second term is the cost from the input dimension d which in uences the number of candidate pairs Ω −k at each iteration k. Table 1 illustrates the number of quantile estimations in function of the dimension d and the number of in uential pairs K. We observe that the computational cost increases signi cantly faster when K = than when K = .  Extensions can be implemented to reduce the computational cost such as removing from Ω, the pairs that are not su ciently improving the minimization.

Applications
The previously proposed methodology is applied to a toy example and a real industrial case-study. It is worth to mention that these experiments (and future ones) can be conducted again using the Python library dep-impact [8], in which are encoded all the procedures of estimation and optimization presented here.

. Numerical example
We consider and extend the portfolio example provided and illustrated in Subsection 2 of Appendix C . The numerical model η is now de ned by the weighted sum where the β = (β , . . . , β d ) is a vector of constant weights. The margins of the random vector X follow the same generalized Pareto distribution with scale σ and shape parameter ξ . Note that the bivariate example in Subsection 2 of Appendix C considers β = and the distribution parameters as σ = and ξ = . In the following examples, we aim at minimizing the median (α = . ) of the output distribution. We chose to x the marginal distribution's parameters at σ = and ξ = . , and we set the constant vector β to a base-10 increasing sequence such that β = ( /d , /d , . . . , ). This choice of weights aims to give more in uence to the latest components of X on Y. Thus, some pairs of variables should be more important in the minimization of the output quantile, as required by the sparsity constraint. We also took n large enough to estimate the output quantile with correct precision on this particular case-study (i.e. n = , ). For all these experiments the results from the di erent methods can be compared.
• Method 1: the grid-search approach with an optimized LHS sampling [46] inside Θ and a random vine structure, • Method 2: the iterative algorithm with an increasing grid-size of N k = (k + ) .
Method 1 is established with the same computational budget as Method 2.

. . Dimension 3
In a three dimensional problem, only three pairs of variables (p = ) are involved in the dependence structure. The sampling size of Θ in Method 1 is set to , which is great enough to explore a three dimensional space. The results are displayed on Figure 4: the estimated quantiles from Method 1 (blue dots) with a convex hull (blue dot line) and the quantile at independence (dark point) are provided. It also highlights the minimum estimated quantiles from Methods 1 and 2 which are respectively represented in blue and red points. We also show in green point, the minimum quantile by considering only the Fréchet-Hoe ding bounds. For each minimum, the 95 % bootstrap con dence intervals is displayed in dot lines.
This low dimensional problem con rms the non-monotonic form of the quantile with the dependence parameter, in particular for the variation of the quantile in function of τ , . As expected, the pair X -X is more in uential on the output quantile due to the large weights on X and X . The minimum values obtained by each method are still lower that the results given by an independent con guration. The minimum using Fréchet-Hoe ding bounds is also provided to show that the minimum is not at the boundary of Θ. Method 1 and 2 have very similar minimum results.

. . Dimension 10
To illustrate the advantages of the iterative procedure, we now consider d = .
In this example, we chose to only consider a Gaussian family for the set of pair-copula family candidates. The sampling size for the exploration of Θ in Method 1 is set to , . Experimental results are summarized over Figure 5, by displaying the minimum quantiles in function of the iteration k of Method 2. The quantile at independence is shown in dark line, the minimum estimated quantile from Method 1 is shown in blue line and the other lines are the minimum quantiles at each iteration of the algorithm, all with their 95% bootstrap con dence interval. We display at each iteration the minimum quantiles of each candidate pair in small dots.  The minimum result from Method 1 is even higher than the quantile at independence. This is due to the very large number of pairs (p = ) that makes the exploration of Θ extremely di cult. On the other hand, Method 2 (iterative algorithm) is de nitely better and signi cantly decreases the quantile value even at the rst iteration (for only one dependent pair). The results are slightly improved with the iterations. We observe at the last iteration that the results from the candidate pairs are slightly higher than the minimum from the previous iteration. It is due to the choice of N k which does not increases enough with the iterations to correctly explore Θ k , which also increases with the iterations.

. . Using multiple pair-copula family candidates
To show the importance of testing multiple copula families, we consider d = and three tests of Method 2 (iterative procedure). The gure 6 shows the minimum from the iterative results using three sets of family candidates: a set of Gaussian and Clayton in red (B = {G, C}), Gaussian only in green (B = {G}), and Clayton only in yellow (B = {C}). We also display below the iteration number, the selected family for B .
At iteration k = , the algorithm with the set B has selected the Gaussian copula as the selected pair and the result is as expected equivalent as for the set B . The next iteration, a Clayton copula has been selected for algorithm with the set B , which slightly improves the minimization compared to the others. The improvement start at iteration k = where the Algorithm with the set B minimizes more the output quantile than the other sets with only one copula family. At the last iteration, the algorithm with set B selected a mix between Gaussian an Clayton families. This diversity seems to lead to better results than using only one family for every pairs. Testing multiple families is an interesting feature of the algorithm and is something that cannot be feasible for the grid-search approach. However, the cost for B is twice larger than for the other methods.

. Industrial Application . . Context
We consider an industrial component belonging to a production unit. This component must maintain its integrity even in case of an accidental situation. For this reason, it is subject to a justi cation procedure by regulation authorities, in order to demonstrate its ability to withstand severe operating conditions. This undesirable event consists in the concomitance of three di erent factors: • the potential existence of small and undetectable manufacturing defects ; • the exposition of the structure to an ageing phenomenon harming the material which progressively diminishes its mechanical resistance throughout its lifespan ; • the occurrence of an accidental event generating severe constraints on the structure.
If combined, these three factors might lead to the initiation of a crack within the structure. Since no failure was observed until now, a structural reliability study should be conducted to check the safety of the structure. To do so a thermal-mechanical code η : R d → R + was used, which calculates the ratio between the resistance and the stress acting on the component during a simulated accident. The numerical model depends on parameters a ected by uncertainties quanti ed throughout numerous mechanical tests. Nevertheless, these experiments are mostly established individually and only few experiments involve simultaneously two parameters.

. . Probabilistic model
For this problem, we introduce d = random variables with prede ned marginal distributions (P j ) j= ...d . The dependence structure is however unknown. From the 15 pairs of variables, only the dependencies of two pairs are known: one is independent and the other follows a Gumbel copula with parameter 2.27. Therefore, we consider p = pairs of variables with unknown dependencies. Given expert feedbacks, we restricted the exploration space Θ by de ning bounds for each pair of variables (i, j) ∈ Ω such that where τ − i,j and τ + i,j are respectively the upper and lower Kendall's correlation coe cient bounds for the dependence of the pair (i, j) and Tc i,j is the transformation from Kendall's tau value to the copula parameter for the associated copula c i,j . This choice enables to explore only realistic dependence structures. For these experiments we only considered Gaussian copulas.

Remark 1.
Note that such bounds can arise from expert knowledge based on a phenomenologic law linking two variables: based on a physical property as the ideal or real gas laws, or the weakest link fracture theory linking material toughness to temperature, a variable is known to be either negatively or positively correlated to another one. Methodologies have been developed in various elds to re ne such assessments [14,67], but were not considered in this research work.

. . Results
We consider the quantile at α = . as a quantity of interest. A rst experiment is established with the incomplete probability structure: only the two pairs with a known dependence structure and all others at independence. Two other experiments are established: an exhaustive grid-search approach with a given vine structure and an iterative procedure with a maximum budget equivalent to the grid-search. A grid-size of is chosen with n = , . The results are displayed in Figure 7. and has the same description as Figure 5. The quantile for the incomplete probability structure is approximately at 1.8. The grid-search and the iterative approaches found dependence structures leading to output quantile values close to 1.2 and 1.1 respectively. The minimum quantile from the iterative procedure is slightly lower than the grid-search approach. The problem dimension is not big enough to create make a signi cant di erence between the methods. However, the resulting dependence structure from the iterative method is greatly simpli ed with only four pairs of variables, in addition to the already known pair. This result highlights the risk of having an incomplete dependence structure in a reliability problem. In this application, the critical limit (safety margin) of the considered industrial component is 1. With the incomplete distribution of X, the output quantile is very high compared to the critical limit and states a high reliability of the component. Unfortunately, if we consider worst-case dependence structures, the output quantile is signi cantly minimized and becomes closer to the critical limit. Thus, if the true dependence structure is close to the obtained worst case dependence structure, the risk of over estimating the output quantile can be important.

Conclusion and discussion
This article considers the situation where one wishes to detect and determine a penalizing correlation structure between the input random variables of a deterministic computer model, when no information is available, through data, on this correlation structure. The meaning given to the notion of penalization is as follows: a low-order quantile of the output distribution of the computer model must be as small as possible, under the constraint that the marginal distributions of the inputs are known. This work takes a rst step towards this goal, by proposing a greedy, heuristic algorithm that explores a set of possible dependencies, taking advantage of the pair-copula construction (PCC) of multivariate probability distributions. Results of experiments conducted on toy and a real models illustrate the good behavior of the procedure: in situations where the monotonicity of the output quantile with respect to the inputs is postulated, a minimum quantile value is obtained using Fréchet-Hoe ding bounds. Nonetheless, it is possible to exhibit situations where the algorithm detect more penalizing dependence structures.
The tools proposed in this paper may be particularly useful in some structural risk and reliability studies, where the inputs typically correspond to independently measured physical parameters, the computer model implements complex structural equations, and the resulting quantile is compared to a safety limit (or a historical quantile value). The correlation structure obtained, associated with the marginal distributions, then allows the generation of scenarios exploring the behavioural limits of the global model. We believe that the methodology proposed here can also be applied to other supervised types of decision-helping studies conducted in the more general context of uncertainty quanti cation. For instance, many problems of natural resource management and sustainable development (where the correlations between in uential parameters are little or poorly known) use such quantiles as key indicators and dimensional quantities of resource extraction [10,32].
However, it is essential that the correlation structure obtained remains phenomenologically relevant and not biased by the interactions between input variables within the model. In the case of structural reliability problems, the input parameters can be roughly separated between loads (forcing conditions, e.g. an injection of cooling water) and capacities (structural resistances, e.g. material toughness) [25], which are independent. This independence is known a priori and this constraint must be integrated into the exploration of correlation structures. The use of similar, trivial bounds on the values of correlation coe cients ( § 4.2.2) makes it possible to respect such a constraint.
The methodology presented in this article requires a number of hypotheses and approximations that pave the way for future research. Besides, some perspectives arise from additional technical results.
It would be interesting to improve the statistical estimation of the minimum quantile, given a dependence structure, by checking the hypotheses underlying the convergence results of Theorem 1. Checking and relaxing these hypotheses should be conducted in relation with expert knowledge on the computer model η and, possibly, a numerical exploration of its regularity properties. The grid search estimation strategy promoted in Section 2.3 arises from the lack of information about the convexity and the gradient of θ → G − θ (α). However, the method remains basic and stochastic recursive algorithms, such as the Robbins-Munro algorithm [54], can be proposed and tested as possibly more powerful (faster) alternatives. Moreover, the empirical quantile was considered for theoretical reasons but can be costly for large α. Another estimator could be consider to reduce the estimation cost. Methods like importance sampling (see for instance [55]) are commonly used to estimate small failure probabilities but its use can be more challenging in the case of quantile estimation.
A signi cant issue is the computational cost of the exploration of possible dependence structures. Reducing this cost while increasing the completeness of this exploration should be a main concern of future works. Guiding the exploration in the space of conditional bivariate copulas using enriching criteria and possible expert knowledge can facilitate the minimization. The Algorithm 2 can also be improved using nonparametric bootstrap. This would quantify the estimation quality of the selected minimum quantile of each iteration. Note however that a seducing feature of an iterative procedure is the a priori possibility of its adaptation to situations where the computational model η is time-consuming. In such cases, it is likely that Bayesian global optimization methods based on replacing the computer model by a surrogate model (e.g., a kriging-based meta-model) [52] should be explored, keeping in mind that nontrivial conservative correlations -losses of quantile monotonicity-can be due to edge e ects (e.g., discontinuities) characterizing the computational model itself. In this sense, mixing dynamic exploration designs with surrogate approaches, in the continuity of the active learning approach introduced in [22,23], seems a relevant research avenue.
We noticed in our experiments on real case-studies that expert knowledge remains di cult to incorporate otherwise that using association and concordance measures. Precise bounds on association measures between two variables are sometimes complex to obtain from the experts due to the lack of bivariate observations. However it is often possible to derive crude upper and lower bounds for association measures, for instance positive/negative correlation or a weak/strong association, and reducing greatly the exploration of the space of correlations. The guidance provided by such expert knowledge can certainly be improved, depending on the situation considered, by adapting dedicated methodologies evoked in Remark 1. Especially, some of these methodologies are based on representation tools of the features of sampled multivariate distributions, that provide intelligible diagnostics, and participate in expert training, a key methodological requirement for Bayesian analysis [5]. Doing so, these new bounds should be considered as probabilistic rather than deterministic, which would lead to adapt the algorithmic core proposed in the present article.
Finally, another approach to consider could be to address the optimization problem (2) within the more general framework of optimal transport theory, and to take advantage of the many ongoing works in this area. Indeed, the problem (2) can be seen as a multi-marginal optimal transport problem (see [53] for an overview). When d = , it corresponds respectively to the classical optimal transport problems of Monge and Kantorovich [66]. However, the multimarginal theory is not as well understood as for the bimarginal case, and developing e cient algorithms for solving this problem remains also a challenging issue [53].

A Proof of the consistency result
The consistency of the estimatorθ requires some regularity of the function θ → G θ . This regularity can be also expressed in term of modulus of increase of the function θ → G − θ (α), on which some useful de nitions and connections with the modulus of continuity are reminded.

A. Modulus of increase of a cumulative distribution function
Let us recall that a modulus of continuity is any real-extended valued function ω : [ , ∞) → [ , ∞) such that lim δ→ ω(x) = ω( ) = . The function f : R → R admits ω as modulus of continuity if for any (x, x ) ∈ R , Similarly, for some x ∈ R, the function f admits ω as a local modulus of continuity if for any x ∈ R , To control the deviation of the empirical quantile in the proof of Proposition 1 further, we consider the modulus of continuity of the quantile functions G − : [ , ] → R where G is a distribution function on R. The quantile function being an increasing function, the exact local modulus of continuity of the quantile function G − at α ∈ ( , ) can be de ned as In the proof of Proposition 1, we note that the continuity of a quantile function G − can be connected to the increase of the distribution function G (see also for instance Section A in [9]). Using the fact that the distribution function is increasing, we introduce the local modulus of increase of the distribution function ϵ G at y = G − (α) ∈ R as: ϵ G (δ, y) := min G(y + δ) − G(y), G(y) − G(y − δ) .

A. Proofs
The estimatorθ de ned in (12) is an extremum-estimator (see for instance Section 2.1 of [51]). The main ingredient to prove the consistency of this estimator is the uniform convergence in probability of the families of the empirical quantiles ( G − θ (α)) θ∈Θ Kn over the family of grids Θ Kn . Thus, We have shown that any ϵ ∈ ( , max(( − α), α)), for any δ > and for any θ ∈ Θ, We now prove the proposition. For any n ≥ and any ϵ > , we have Let ξ , . . . , ξn be n i.i.d. uniform random variables. The uniform empirical distribution function is de ned by

The inverse uniform empirical distribution function is the function
The empirical distribution function G θ can be rewritten as (see for instance [65]): and as well for the quantile function, From Inequality (24), we obtain By de nition of the local modulus of continuity ω G − θ of the quantile function G − θ at α, we have Therefore, by replacing (26) in (25) and using (23), we obtain Assumption C then yields The DKW inequality [21] gives an upper bound of the probability of an uniform empirical process {|Un(α) − α|}.
As well for an uniform empirical quantile process {|U − n (α) − α|} (see for example Section 1.4.1 of [16]), such as ∀λ > : Moreover, [45] proved that one can take C = . Therefore, Equation (27) can be bounded using the DKW and A second requirement to get the consistency of the extremum estimator is the regularity of θ → G − θ (α). This is shown in the next proposition.

Proposition 2. Under Assumptions A, B and C, the function
is continuous in θ over Θ.
Proof 2 (Proof of Proposition 2). According to Assumption A, for any θ ∈ Θ, the distribution P θ admits a density function f θ with respect to the Lebesgue measure on R d such that where f j is the marginal density function of X j , for j = , . . . , d and the Lebesgue measure on R. Moreover, for any x ∈ R d , the function θ → f θ (x) is continuous in θ over Θ.
The domain Θ × [ , ] p is a compact set and according to Assumption A, there exists a constantc such that ∀(θ, u) ∈ Θ × [ , ] d , c θ (u) ≤c. Consequently, we have For θ ∈ Θ and for any h > , we denote y h = G − θ+h (α). According to Assumption B we have α = G θ+h (y h ) and thus, Now, using Assumption C, we have that G θ is strictly increasing in the neighborhood of G − θ (α) and thus G − θ is continuous in the neighborhood of α. Note that We then apply a standard dominated convergence theorem using (28) to get that This, with (29) and with the continuity of θ → G θ , shows that We are now in position to prove Theorem 1.
We now assume that Assumption D is also satis ed. Let θ * be the unique minimizer of θ → G − θ (α). Let h > such that B(θ * , h) c := {θ ∈ Θ : θ − θ * ≥ ε} is not empty. According to Proposition 2 and using the fact that Θ is compact, we have sup Consequently, for any ∀h > small enough, there exists ϵ > such that Let h > and take ϵ such that (33) is satis ed for h. According to Proposition 1, G − θ (α) − G − θ (α) tends to zero in probability as n tends to in nity. This, with (13), shows that We conclude using (33).  The pair-copula in the rst tree characterize pairwise unconditional dependencies, while the pair-copula in higher order trees model the conditional dependency between two variables given a set of variables. The number of conditioning variables grows with the tree order. Note that a PCC where all trees have a path-like structure de ne the D-vine subclass while the star-like structures correspond to C-vine subclass. All other vine structures are called regular vines (R-vines) [6].
We illustrate the concept of a vine model with a d = dimensional example. For clarity reasons, we use the same simpli cations as in Section 3.2 which consider for instance f = f (x ), f = f (x ) and so on for higher order and conditioning. One possible PCC can be written for this 5-dimensional con guration:  Figure 8: R-vine structure for d = .
The vine structure associated to (34) is illustrated in Figure 8. This graphical model considerably simplify the understanding and we observe that this model is a R-vine because there is no speci c constraints on the trees.
A re-labeling of the variables can lead to a large number of di erent PCC. [49] calculated the number of possible vine structures with the dimension d and shows that it becomes extremely large for high dimension problems. We illustrate below, using the same d = dimensional example, two other PCC densities: where (35) and (36)  An e cient way to store the information of a vine structure is proposed in [48] and is called a R-vine array. The approach uses the speci cation of a lower triangular matrix where the entries belong to , . . . , d. Such matrix representation allows to directly derive the tree structure (or equivalently the associated PCC distribution). For more details, see [48].

B. Generating R-vine from an indexed list of pairs
The iterative procedure proposed in Section 3.3, described by Algorithm 2, minimizes the output quantile by iteratively determining the pairs of variables that in uences the most the quantile minimization. At each iteration of the algorithm (step 1.a), a new vine structure is created by considering the list of in uential pairs. The speci city of this vine creation is to consider the ranking of the list by placing the most in uential pairs in the rst trees of the R-vine. Thus, we describe in this section how to generate vine structure with the constraint of a given list of indexed pairs to ll in the structure.
(a) D-vine structure for d = .
(b) C-vine structure for d = .

B. . The algorithm
We consider the same notation as in Algorithm 2. Creating a vine structure from a given indexed list of pairs Ω k is not straightforward. The di culties come from respecting the ranking of Ω k and the respect of the R-vine conditions. Indeed, the pairs cannot be append in the structure easily. The vine structure must respect these conditions, which can be sometime very restrictive. The procedure we proposed is detailed by the pseudocode of Algorithm 3 and can be greatly simpli ed in these few key steps: 1. ll V with the list Ω k , 2. ll V with a permutation of Ω −k , 3. if V is not a R-vine, then permute Ω k and restart at step 1.
In step 1 and 2, the lling procedure, detailed in Algorithm 4, successively adds the pairs of a list in the trees of a vine structure. Adding a pair (i, j) in a tree T l associates (i, j) with the conditioned set and determine a possible conditioning set D from the previous tree such as a possible edge is i, j|D.
In step 2, because the ordering of Ω −k is not important in the lling of V, the permutation of Ω −k aims at nding a ranking such as V leads to a R-vine.
In step 3, when the previous step did not succeeded and the resulting V is not a R-vine structure, then the ranking of Ω k is not possible and must be changed. The permutation of some elements of Ω k must be done such as the ranking of the most in uential pairs remains as close as possible to the initial one.

B. . Example
For illustration, let's create a d = dimensional vine structure with the given list of pairs Ω k = (( , ), ( , ), ( , ), ( , ), ( , ), ( , )) using Algorithm 3. Using the original list Ω k , the Fill function may fail at line 7 of Algorithm 3, and more precisely, at line 15 of Algorithm 4. Indeed, the rst tree of V does not validate the R-vine conditions. The tree is illustrated in Figure 10 and as we can see, the nodes are not all connected into one single tree. Therefore, we permuted the list Ω k by exchanging the pairs ( , ) and ( , ), as shown in Figure 11. This permutation now leads to a vine structure that respects the new ranked list Ω k = (( , ), ( , ), ( , ), ( , ), ( , ), ( , )) .

C Supplementary material on copulas
This supplementary material is dedicated to a preliminary exploration of the in uence of copula structure on the behavior of the worst quantile, illustrated with toy examples. Especially, while it could be expected that G − θ (α) is a monotonic function with θ, and that the minimum can be reached for a trivial copula (i.e., reaching the Fréchet-Hoe ding bounds). Our experiments show that this behavior is not systematic.

C. About the copula choice
One of the most common approaches to model the dependence between random variables is to assume linear correlations feeding a Gaussian copula. In this case, the problem is reduced by determining the correlation matrix of X that minimizes G − θ (α). However, the positive semi-de nite constraint on the correlation matrix makes the exploration di cult and the minimization harder when the problem dimension increases. Moreover, such a Gaussian assumption is very restrictive and is inappropriate for simulating heavy tail dependencies [44]. Still in this elliptical con guration, the t-copulas [19] can be used to counterpart these problems. Nevertheless, tail dependencies are symmetric and with equal strengths for each pair of variables. Another alternative is to consider multivariate Archimedean copulas [47] which are great tools to describe asymmetric tail dependencies. However, only one parameter governs the strength of the dependence among all the pairs, which is very restrictive and not exible in high dimension. For a same correlation measure between two random variables, multiple copulas can be tted and lead to a di erent distribution of Y.
It is clear that the copula choice of X has a strong impact on the distribution of Y (see for instance [60]). Therefore, various copula types should be tested to determine the most conservative con guration. In the following, we may consider a exible approach setting by modeling the input multivariate distribution using regular vine copulas (R-vines). The necessary basics of R-vines are introduced in Section 3 of the article and detailed in its Appendix B.

C. About the monotony of the quantile
For many simple case studies case studies, the worst quantile is reached for perfect dependencies (Fréchet-Hoe ding bounds). More generally, when the function has a monotonic behavior with respect to many variables, it is likely that the minimum output quantile is reached at the boundary of Θ. This phenomenon is observed for various physical systems. To illustrate this phenomenon, we consider a simpli ed academic model that simulates the over ow of a river over a dike that protects industrial facilities. The river over ow S is described by such as, when S < , a ooding occurs. The involved parameters of (37) are physical characteristics of the river and the dike (e.g., ow rate, height of the dike) which are described by random variables with known marginal distributions. See [36] for more information. For a given risk α, we aim at quantifying the associated over ow's height describe by the α-quantile of S. We extend this model by supposing that the friction (Strickler-Manning) coe cient Ks and the maximal annual ow rate Q are dependent with an unknown dependence structure. To show the in uence of a possible correlation between Ks and Q on the quantile of S, we describe their dependence structure with multiple copula families. Figure 12 shows the variation of the estimated quantile of S (with a large sample size) in function of the Kendall coe cient τ between Ks and Q for di erent copula families. We observe di erent slopes of variation for the di erent copula families, with lower quantile values for the copulas with heavy tail dependencies (i.e., Clayton, Joe). At independence (τ = ) and for the counter-monotonic con guration (τ = − ),the quantile values of these families are obviously equivalent. This variation is slight and the quantile is still above zero, but this shows how the dependencies can in uence the results of a reliability problem. This illustration shows that the minimum is reached at the boundary of the exploration space, where the two variables are perfectly correlated.
We can take advantage of this observation to speed up the algorithms presented in the next sections by exploring only the boundaries of Θ. However, assuming that the minimum is reached on the boundary of Θ is a strong assumption that can be unsatis ed in some applications. See Fallacy 3 of [24] for a highlight of this pitfall.
To illustrate this statement, we now give a counter example in the bidimensional setting. We assume uniform marginal distributions for the input such that X ∼ U(− , ) and X ∼ U(− , ), and we consider the model function The same experience as for Figure 12 is established and the results are shown in Figure 13. The slopes of the quantile estimations with the Kendall coe cient, for each copula families, are quite di erent than the results of Figure 12. We observe that the quantile is not monotonic with the Kendall coe cient and its minimum is not reached at the boundary, but for τ ≈ . . Moreover, the Gaussian copula is the family that minimizes the most the quantile. It shows that copulas with tail dependencies are not always the most penalizing. A second example, inspired from Example 6 of [24], also shows that the worst case dependence structure in an additive problem is not necessary for perfectly correlated variables. We consider a simple portfolio optimization problem with two random variables X and X with generalized Pareto distributions such as F (x) = F (x) = x +x . We aim at maximizing the pro t of the portfolio, which is equivalent as minimizing the following additive model function η(X , X ) = −(X + X ).
We consider the median (α = . ) of the output as an e ciency measure. Figure 14 shows the output median in function of the Kendall coe cient τ between X and X . Just like the previous example, we observe a non-monotonic slope of the median in function of τ. The variation can be signi cant and the minimum is obtained at τ ≈ . for the heavy tail copula families (i.e., Clayton and Joe). The phenomenon can be explained by the marginal distributions of the random variables, which are close Pareto distributions. A large correlation seems to diminish the in uence of the tails, which gives a higher quantile value. This explains why the minimum is obtained for a dependence structure other that independence or the perfect dependence. Therefore, these examples show that the worst quantile can be reached for other con gurations than the perfect dependencies.