Clement Gauchy

Jerome Stenger

Roman Sueur

Bertrand Iooss
email: bertrand.iooss@edf.fr

An Information Geometry Approach to Robustness Analysis for the Uncertainty Quantification of Computer Codes

Keywords: Computer experiments, Density perturbation, Fisher metric, Importance sampling, Quantile, Sensitivity analysis

come

Introduction

Over the last few decades, two major trends in industrial and research practices have led to a rise in the importance of uncertainty quantification (UQ) methods (de Rocquigny et al., 2008;[START_REF] Smith | Uncertainty quantification[END_REF]Ghanem et al., 2017). The first is the replacement of full-scale physical experiments, considered costly and difficult to implement, by numerical models. This choice raises the issue of a potential mismatch between computer codes and the physical reality they aim to simulate. The second trend involves accounting for risk in an increasing number of industrial activities, and in particular, doing it quantitatively.

In both situations, UQ can be performed by considering the imperfectly known inputs of a computer code (given as a function G(•)) as a vector of random variables X = (X 1 , . . . , X d ). The most widespread approach consists of running G(•) with different combinations of inputs in accordance with their range of plausible values, in order to study the associated uncertainty found in the output Y = G(X 1 , . . . , X d ), or, in other settings, to estimate a specific quantity of interest (QoI). The latter means some statistical quantity derived from Y , e.g,. the performance, via the mean of Y , or some risk criterion in terms of a high-level quantile.

For instance, the nuclear industry faces major issues as facilities age and regulatory authorities' requirements toughen [START_REF] Bucalossi | Comparison between best-estimate-plus-uncertainty methods and conservative tools for nuclear power plant licensing[END_REF][START_REF] Mousseau | Uncertainty quantification in a regulatory environment[END_REF]. One precise case is where operators have to study loss of coolant accidents (LOCA), which result in a break in the primary loop of pressurized water nuclear reactors. This scenario can be simulated using system thermal-hydraulic computer codes which involve dozens of physical parameters including such things as condensation and heat transfer coefficients [START_REF] Mazgaj | Comparison of CATHARE results with the experimental results of cold leg intermediate break LOCA obtained during ROSA-2/LSTF test 7[END_REF][START_REF] Sanchez-Saez | Uncertainty analysis of large break loss of coolant accident in a pressurized water reactor using non-parametric methods[END_REF]. However, the values of many such parameters are known with only limited precision [START_REF] Larget | How to bring conservatism to a BEPU analysis[END_REF] as they are typically calculated by way of other quantities measured via small-scale physical experiments. Certain other variables can only be measured during periodic inspections, e.g., the characteristics of pumps in hydraulic systems.

Various methods from the UQ domain can be useful when dealing with such uncertainties in system safety analyses. First, certain methods aim to improve the exploration of the input domain X by using design of experiments, e.g., space filling designs [START_REF] Fang | Design and modeling for computer experiments[END_REF]. These make it possible to cover an input domain as evenly as possible for a fixed number of code runs, and to limit unexplored areas. For the estimation of certain specific QoI such as the probability of threshold exceedance or an α-order quantile of the output, Monte Carlo-type methods are often preferred. In particular, accelerated Monte Carlo methods (e.g., importance sampling or subset simulation) target the most informative areas of X in the sampling algorithm in order to estimate the QoI while controlling its estimation error [START_REF] Morio | Estimation of rare event probabilities in complex aerospace and other systems[END_REF]. As a preliminary or concomitant stage, global sensitivity analysis is also essential in order to eliminate non-influential parameters and rank influential ones according to their impact on the QoI [START_REF] Iooss | A review on global sensitivity analysis methods[END_REF][START_REF] Iooss | Advanced methodology for uncertainty propagation in computer experiments with large number of inputs[END_REF].

All of these approaches are useful for dealing with the existence of uncertainty in applied problems. However, industrial (e.g., nuclear facility) operators face the difficulty of justifying their risk assessment methods over and above providing simulation results. Risk assessment methods must demonstrate that they provide reliable safety margins. Doing so means a probable overestimation of the likely risk after accounting for all known sources of uncertainty. This principle of conservatism, which can be easily implemented when dealing with very simple monotonic physical models, can be hard in practice to adapt to computer codes simulating complex and non-monotonic physical phenomena. Indeed, it is rarely simple to apply this principle when implementing UQ methods based on a set of computer experiments that provide a whole range of values for the output quantity Y .

To address this issue, a new UQ branch of robustness analysis has emerged in recent years in the field of sensitivity analysis. It consists of evaluating the impact of the choice of the inputs' distributions and, more precisely, involves an analysis of how values of the QoI vary with respect to these choices.

An initial strategy is to consider a whole set of input distributions and analyze the related output distributions. For global sensitivity analysis, [START_REF] Hart | Robustness of the Sobol' indices to marginal distribution uncertainty[END_REF] use "optimal perturbations" of probability density functions to analyze the robustness of variance-based sensitivity indices (called Sobol indices [START_REF] Sobol | Sensitivity estimates for non linear mathematical models[END_REF]). [START_REF] Meynaoui | New statistical methodology for second level global sensitivity analysis[END_REF] and [START_REF] Chabridon | Reliability-based sensitivity estimators of rare event probability in the presence of distribution parameter uncertainty[END_REF] propose approaches that deal with "second-level" uncertainty, i.e., uncertainty in the parameters of the input distributions. An alternative approach known as "optimal uncertainty quantification" avoids specifying the input probability distributions, transforming the problem into the definition of constraints on moments [START_REF] Owhadi | Optimal Uncertainty Quantification[END_REF][START_REF] Stenger | Optimal Uncertainty Quantification of a risk measurement from a thermal-hydraulic code using Canonical Moments[END_REF]. The latter is out of the scope of the present work; here we consider that the initial probability distributions input by the user remain important.

In real-world uncertainty quantification studies in engineering, input distributions are typically truncated, as we are dealing with physical parameters that in reality belong to finite-dimensional compact sets. It is therefore natural to assume that there is no uncertainty in the support of input random variables. In this paper, we also assume that the inputs are mutually independent. Keeping in mind that our goal is to directly deal with input distributions (and not to consider second-level uncertainty), one particularly interesting solution has been proposed in the context of reliability-oriented sensitivity analysis by [START_REF] Lemaître | Analyse de sensibilité en fiabilité des structures -Sensitivity analysis in structural reliability[END_REF] (see also [START_REF] Lemaître | Density modification based reliability sensitivity analysis[END_REF]; [START_REF] Sueur | Perturbed-Law based sensitivity indices for sensitivity analysis in structural reliability[END_REF]): perturbed-law based indices (PLI). A density perturbation consists of replacing the density f i of one input X i by a perturbed one f iδ , where δ ∈ R represents a shift in the value of a specific moment (e.g., the mean or variance). Among all densities which have an equivalent shift in the chosen moment, f iδ is defined as that which minimizes the Kullback-Leibler divergence from f i . This method has been applied to compute a QoI that corresponds to a failure probability [START_REF] Iooss | Uncertainty and sensitivity analysis of functional risk curves based on Gaussian processes[END_REF][START_REF] Perrin | Efficient estimation of reliability-oriented sensitivity indices[END_REF], a quantile (Sueur et al., 2017;[START_REF] Larget | How to bring conservatism to a BEPU analysis[END_REF] and a superquantile [START_REF] Iooss | BEPU robustness analysis via perturbedlaw based sensitivity indices[END_REF][START_REF] Larget | Increasing conservatism in BEPU IB LOCA safety studies using complementary and industrially cost effective statistical tools[END_REF].

However, this method is not fully satisfactory. Indeed, the interest of perturbing moments is arguable: perturbing the mean or variance of input distributions is essentially an arbitrary decision. Moreover, it has been observed that applying the same perturbation to two different parametric distributions can yield significant differences in the value of the Kullback-Leibler divergence (between the initial and perturbed density).

Another possibility described in [START_REF] Perrin | Efficient estimation of reliability-oriented sensitivity indices[END_REF] is to use an iso-probabilistic operator to transform all input random variables into standard Gaussian ones. This makes perturbations comparable when applied in this new space, but it remains difficult to trans-late this interpretation back into the initial physical space, that which is of interest to users. We note in passing another type of robustness analysis that has been proposed in quantitative finance by [START_REF] Cont | Robustness and sensitivity analysis of risk measurement procedures[END_REF]. These authors investigate whether the estimated QoI is sensitive to small perturbations of the empirical distribution function. To this end, they define the robustness of a QoI in terms of its continuity with respect to the Prokhorov distance on the set of integrable random variables.

The goal of the present paper is to propose a novel approach to perturbing probability distributions. It relies on density perturbations based on the Fisher distance [START_REF] Costa | Fisher information distance: a geometrical reading[END_REF] as a measure of dissimilarity between an initial density f i and a perturbed one f iδ .

This distance defines a geometry on spaces of probability measures known as information geometry [START_REF] Nielsen | Cramér-Rao Lower Bound and Information Geometry[END_REF]. The statistical interpretation of the Fisher distance provides an equivalence between perturbations of non-homogeneous quantities, and consequently a coherent framework for robustness analysis. Before presenting this approach, we first review existing density perturbation methods in Section 2. Section 3 is dedicated to the description of our method and a discussion of numerical methods. Section 4 illustrates our density perturbation method for the practical robustness index (PLI). We apply this to a toy example function and an industrial case study in Section 5. This is followed by conclusions and research options going forward. The online supplementary materials contain four further sections describing: 1) some numerical results on Fisher spheres computations; 2) the reverse importance sampling method; 3) the theoretical properties of the PLI-quantile estimator; and 4) an application of the PLI to another toy model.

Previous approaches to density perturbation for uncertainty quantification robustness analysis

The method of [START_REF] Lemaître | Density modification based reliability sensitivity analysis[END_REF], later known as PLI by [START_REF] Sueur | Perturbed-Law based sensitivity indices for sensitivity analysis in structural reliability[END_REF], is based on the idea of perturbing input densities. It aims to provide a practical counterpart to the general idea of analyzing the output QoI of a model in a UQ framework when one or several parameters of the input probabilistic model (considered as the reference) is/are modified. This strategy can be seen as a way to take into account an "error term" that could be added to an imperfectly-known input distribution.

Kullback-Leibler divergence minimization

To construct a perturbed distribution f iδ from a distribution f i , the approach of [START_REF] Lemaître | Density modification based reliability sensitivity analysis[END_REF] is nonparametric. This method essentially aims to analyze perturbations in the most common characteristics of input distributions such as the mean or variance. To illustrate this for the mean, let us assume the random variable

X i ∼ f i has mean E[X i ] = µ.
By definition, the perturbed density will have a mean of µ + δ. However, this is clearly insufficient to fully determine the perturbed distribution and in particular explicitly access the value of f iδ over the whole domain of X i . Amongst all densities with a mean of µ + δ, f iδ is defined as the solution of the minimization problem:

f iδ = arg min π∈P, s.t Eπ[X i ]=E f i [X i ]+δ KL(π||f i ) , (1)
where P is the set of all probability measures which are absolutely continuous with respect to f i . This approach basically consists of perturbing the chosen parameter while chang-ing the initial model as little as possible. Under this definition, "changing" the model corresponds to increasing the entropy, where the Kullback-Leibler divergence between two densities f and π is given by:

KL(π||f ) = log π(x) f (x) f (x) dx . (2)
This method can be applied to higher-order moments (for instance moments of order 2-to define a variance perturbation) and, more generally, to constraints that can be expressed as a function of the perturbed density, such as quantiles [START_REF] Lemaître | Analyse de sensibilité en fiabilité des structures -Sensitivity analysis in structural reliability[END_REF]. Notice that in the case of an initial Gaussian distribution, the perturbed distribution remains Gaussian with simply a change in the mean or variance.

In the general case, this method has several drawbacks. First, the likelihood ratio between f iδ and f i (see section 4.1) may not have an analytic form, which leads to numerical difficulties. Second, this method requires that moments be defined for the initial density.

Third, the main difficulty concerns the interpretation of results. Indeed, each input containing uncertainty in the UQ model is perturbed using a range of values of δ. However, the choice of the moment to be perturbed remains arbitrary. Such problems are particularly worrisome when uncertainty in the input parameter is epistemic [START_REF] Kiureghian | Aleatory or epistemic? Does it matter?[END_REF] (meaning it comes from a lack of knowledge). Perturbing only the mean or variance of such an input distribution may be overly limiting when trying to satisfy regulatory authorities focused on being overly conservative.

We also recall that all input random variables are assumed to be mutually independent.

Nonetheless, the effect of perturbations can be considered only for each variable individually, and in absolute terms (the same δ shift may have a quite different impact on different input densities). This method thus proves challenging when one wants to compare the relative impact of perturbations on different inputs.

Standard space transformations

To interpret a δ shift on an input distribution and in particular to enable inputs to be compared with respect the impact on the QoI of the same perturbation, an equivalence criterion between inputs is required. An idea developed by [START_REF] Perrin | Efficient estimation of reliability-oriented sensitivity indices[END_REF] consists of applying perturbations in the so-called standard space (instead of the initial physical space) in which all input distributions are identical, thus making all perturbations equivalent. Then the perturbed densities are obtained by applying a reverse transform.

In the case of independent inputs, the required distribution transform is a simple inverse probability one. Given a real-valued random vector X = (X 1 , . . . , X d ) T with cumulative distribution function (cdf) F , the transform is the random vector S = Φ -1 (F (X)), where Φ is the cdf of the standard Gaussian distribution N (0, I d ). Consequently, S follows a standard Gaussian distribution whatever the initial distribution F . In the Kullback-Leibler divergence minimization framework (see Section 2.1), a perturbation of the mean in the standard space simply consists of a mean shift without the standard deviation changing.

For instance, for a mean shift of δ ∈ R d , the perturbed distribution is defined by the inverse transformation F δ = F -1 (Φ(S + δ)) in order to get back from the standard space to the physical space. This therefore leads to an analytic expression for the perturbed density f iδ in the one-dimensional case via the change of variable formula (Stirzaker, 2003, p.318):

f iδ (x) = e -δ 2 +2δΦ -1 (F i (x)) 2 f i (x) , (3)
where F i is the cdf of the variable X i for i ∈ {1, . . . , d}. This simple formula makes the perturbed density and likelihood ratio easy to compute.

However, similar perturbations in the standard space can lead to quite different ones in the physical space, depending on the initial distributions. As an example, Figure 1 depicts the Kullback-Leibler divergence (approximated with Simpson's rule (Abramowitz and Stegun, 1974)) of two specific distributions (the Triangular1 T (-1, 0, 1) and the Uniform U[-1, 1]) with their respective distributions in the standard space as the value of δ changes. We see that the Kullback-Leibler divergence can behave very differently for different distributions in the physical space, even though the exact same perturbation has been applied in the standard space. Furthermore, there is no general rule for estimating the mean of a perturbed physical input in its physical space for a given mean perturbation of the input in the standard space. Such difficulties are amplified when considering perturbations involving moments other than the mean. For instance, there is no general equivalence in the physical space between perturbations applied to the mean and standard deviation of the same input probability distribution in the standard space. Overall, it would seem generally difficult, if not impossible, to have a simple way to convert results provided by this method into a relationship between input and output physical quantities, making such results difficult to interpret.

0.0 0.5 1.0 1.5 2.0 δ 0 10 20 30 KL divergence T (-1, 0, 1) U ([-1, 1])
Figure 1: The Kullback-Leibler divergence between the initial distribution and the perturbed one for perturbations of the mean, δ ∈ [0, 2]. Further details on perturbed distributions can be found in [START_REF] Lemaître | Analyse de sensibilité en fiabilité des structures -Sensitivity analysis in structural reliability[END_REF] and [START_REF] Lemaître | Density modification based reliability sensitivity analysis[END_REF].

A perturbation method based on information geometry

The Kullback-Leibler divergence can be interpreted as the power of a hypothesis test with null hypothesis: "X i follows the distribution f i " and alternative hypothesis: "X i follows the distribution f iδ " [START_REF] Eguchi | Interpreting Kullback-Leibler Divergence with the Neyman-Pearson Lemma[END_REF]. For this reason, it would appear to be an appropriate tool to measure how far a perturbed density is from the initial one and thus provides a formal counterpart to the vague idea of "uncertainty in the distribution".

However, in a practical robustness analysis context, its use often implies embedding inputs in a standard Gaussian space to begin with, in order to allow a comparison of perturbations of different input distributions, as well as to enable easy computation of perturbed densities. This transformation of the inputs leads to manipulating perturbations in a non-physical input space, making it impossible to interpret them in terms of how they modify physical inputs' densities.

The Fisher distance

To try to form an intuitive understanding of the consequences of such perturbations on the distribution of the output variable Y , it is necessary to base our perturbation method on a metric which allows us to compare perturbations coming from different inputs distributions in the UQ model. For instance, if one input is Gaussian and another is log uniform, their respective perturbations (associated with the same level of perturbation) should be able to be interpreted in the same way. In particular, the perturbation method should result in identical perturbed densities when applied to two different parametric representations of the same input distribution. The Fisher distance [START_REF] Rao | Information and the accuracy attainable in the estimation of statistical parameters[END_REF] is consistent with these wishes. It is based on the local scalar product induced by the Fisher information matrix in a given parametric space, and defines a Riemannian geometry on the corresponding set of probability measures as on any Riemannian manifold with its associated metric.

Consider the family of parametric densities S = {f θ , θ ∈ Θ ⊂ R r }. We recall that every input variable represents a physical parameter with known domain of validity, and therefore, for all θ in Θ, the support of f θ is assumed to be a compact set of R. The metric associated with the coordinate function θ, called the Fisher (or Fisher-Rao) metric, is defined as:

I(θ 0 ) = E ∇ θ log f θ (X)| θ=θ 0 (∇ θ log f θ (X)| θ=θ 0 ) T ,
where I(θ 0 ) is the Fisher information matrix evaluated at θ 0 for this statistical model.

The Fisher information, well-known in the fields of optimal design, Bayesian statistics, and machine learning, is a way of measuring the amount of information that an observable random variable X carries about an unknown parameter θ of the distribution of X. The Fisher information matrix defines the following local inner product in S for u ∈ R r and

v ∈ R r : u, v θ = u T I(θ)v . (4)
Given two distributions f θ 0 and f θ 1 in a manifold S, a path from f θ 0 to f θ 1 is a piecewise smooth map q : [0, 1] → Θ satisfying q(0) = θ 0 and q(1) = θ 1 . Its length [START_REF] Sternberg | Lectures on Differential Geometry[END_REF] l(q) satisfies the following equation:

l(q) = 1 0 q(t), q(t) q(t) dt , (5)
where q is the derivative of q. Similarly, the energy [START_REF] Sternberg | Lectures on Differential Geometry[END_REF]) E(q) of a path is defined by the equation:

E(q) = 1 0 1 2 q(t), q(t) q(t) dt . (6)
The Fisher distance between f θ 0 and f θ 1 is defined as the minimal length over the set of

paths P(f θ 0 , f θ 1 ) from f θ 0 to f θ 1 : d F (f θ 1 , f θ 2 ) = inf q∈P(f θ 1 ,f θ 2 ) l(q) . (7)
The path γ minimizing this length-or equivalently minimizing the energy-is called a geodesic [START_REF] Costa | Fisher information distance: a geometrical reading[END_REF]. The specific choice of the Fisher information matrix for a Riemannian metric matrix leads to a very interesting statistical interpretation, as shown in Amari (1985, p.27); it is directly related to the Cramer-Rao lower bound [START_REF] Rao | Information and the accuracy attainable in the estimation of statistical parameters[END_REF] which states that, for any unbiased estimator θ of θ, the covariance matrix Var( θ) is bounded by

I(θ) -1 .
This means that the Fisher information is the maximum amount of information about the value of a parameter one can extract from a given sample. More formally, under some regularity conditions (given by (Newey and McFadden, 1994, Theorem 3.3)), if x 1 , . . . , x n are n independent observations distributed according to a density f θ , the maximum likelihood estimator θ n of θ converges weakly to a Gaussian distribution with mean θ and covariance

I(θ) -1 n .
The density of θ n , denoted p( θ n , θ) is written:

p( θ n , θ) = 1 ( (2π) n ) r det(I(θ) -1 ) exp - n( θ n -θ) T I(θ)( θ n -θ) 2 . ( 8
)
When n → +∞, this probability density is proportional to ( 1985), which is the local inner product defined in equation ( 4). This result can be interpreted as follows: the Fisher distance between two distributions f θ and f θ represents the separability of the two distributions by a finite sample of independent observations sampled from the distribution f θ [START_REF] Amari | Differential-Geometrical Methods in Statistics[END_REF].

θ n -θ) T I(θ)( θ n -θ) (Amari,
Let us illustrate the Fisher distance in a simple example. Consider the statistical manifold of univariate normal distributions S = {N (µ, σ 2 ), (µ, σ) ∈ R × R * + }. The Fisher information matrix has the analytical form [START_REF] Costa | Fisher information distance: a geometrical reading[END_REF]:

I(µ, σ) =     1/σ 2 0 0 2/σ 2     . ( 9
)
If we apply the change of coordinates φ(µ, σ) → ( µ √ 2 , σ), the related geometry is the hyperbolic geometry in the Poincaré half-plane [START_REF] Stillwell | Numbers and Geometry[END_REF] in which the geodesic and distance between two normal distributions are known analytically [START_REF] Costa | Fisher information distance: a geometrical reading[END_REF].

Geometrically, geodesics are the vertical lines and half-circle centered on the line σ = 0.

Further details on interpretations of information geometry can be found in [START_REF] Costa | Fisher information distance: a geometrical reading[END_REF]. and C and D. We observe that the Fisher distance between A and B is greater that that between C and D. This illustrates how information geometry provides a well-grounded framework to measure statistical dissimilarities in probability measure spaces. In summary, the Fisher distance provides a satisfactory and coherent grounding for our notion of density perturbation. Before continuing, let us define a perturbation of a density

µ σ A • B • C • D • σ 1 σ 2 µ 1 µ 2 µ 1 µ 2 C D A B
f to be of level δ if the Fisher distance between f and the perturbed density f δ is equal to δ. The set of all perturbations of f at level δ is then the Fisher sphere of radius δ centered at f , whenever this perturbation is applied to one or another of the parameters. This implies that, in this framework, we do not consider one specific perturbed distribution, but rather, a non-finite set of probability densities. The following section is dedicated to the development of a numerical method to compute Fisher spheres of radius δ centered at f .

Computing Fisher spheres

As detailed in Section 3.1, geodesics are defined as the solution of a minimization problem.

Specifically, a geodesic is a path with minimal length or energy (denoted E). Given a smooth map q : [0, 1] → S, we have

E(q) = 1 0 1 2 q(t), q(t) q(t) dt . (10)
In the following we denote L(t, q, q) = 1 2 q(t), q(t) q(t) the Lagrangian of the system. The energy of a path can be rewritten as

E(q) = 1 0 L(t, q, q)dt . (11)
A necessary condition for the path q to minimize the energy E is to satisfy the Euler-Lagrange equation (see [START_REF] Gelfand | Calculus of Variations[END_REF] for details):

∂L ∂q = d dt ∂L ∂ q . ( 12
)
We denote p = ∂L ∂ q and obtain by differentiation of the quadratic form L(t, q, q) = 1 2 qT I(q) q that p = I(q) q and q = I -1 (q)p. Then, inspired by Lagrangian mechanics theory (Arnold, 1997, p.65), the Hamiltonian H(p, q) defined by

H(p, q) = p T q -L(t, q, q) = p T I -1 (q)p -1 2 qT I(q) q = 1 2 p T I -1 (q)p (13)
is constant whenever q is a geodesic. Eq. ( 13) is derived from the Euler-Lagrange equation, and implies that (p, q) follows a system of ordinary differential equations (ODE) known as

Hamilton's equations:

       q = ∂H ∂p = I -1 (q)p , ṗ = - ∂H ∂q = ∂L(t, q, I -1 (q)p) ∂q . ( 14
)
The goal is to find any geodesics q that satisfy q(0) = θ 0 and d F (f, q(1)) = δ. This corresponds to computing the Fisher sphere centered at f θ 0 with radius δ. The only degree of freedom remaining to solve the ODE system ( 14) entirely is the initial velocity p(0). Notice that the Hamiltonian is equal to the kinetic energy, since p = I(q) q. As the Hamiltonian is constant on a geodesic, we have for all t:

1 2 q(t), q(t) q(t) = k , ( 15
)
where k is non-negative. The length of q is therefore equal to

1 0 q(t), q(t) q(t) dt = √ 2k , (16)
whereby δ = √ 2k. Therefore, Eq. ( 13) can be rewritten:

δ = √ 2k ⇐⇒ p T I -1 (q)p = δ 2 . ( 17
)
Taking equation ( 17) at the initial state t = 0, we can determine all of the initial velocities for which d F (q(0), q(1)) = δ. These are needed to solve the system of ODEs ( 14) and compute the geodesics.

Generally, computing the geodesic between two distributions is challenging. Methods relying on shooting algorithms have been developed for this problem (Le Brigant, 2019).

However, in our framework, we don't want to compute the distance between two given distributions, which implies to find the particular geodesic joining these, but we aim at screening a Fisher sphere as well as possible for a fixed distance (as detailed in section 4.2).

In the Section 1 of the online supplementary material, we focus on numerical methods for computing geodesics by solving the system of ODEs ( 14). We then illustrate these methods by computing Fisher spheres on the Gaussian manifold

S = {N (µ, σ 2 ), (µ, σ) ∈ R × R * + }.
4 Applying this to perturbed-law based indices

The UQ robustness analysis detailed in Section 1 and Section 2 aims to quantify the impact of a lack of knowledge about an input distribution on the UQ of model outputs.

In Section 3, a coherent formal definition of density perturbation was proposed. We now illustrate how to use this solution to define a practical robustness analysis method. First however, analyzing the effect of perturbing an input density requires the definition of some criterion or "index" to summarize the effects on a QoI.

Definition of perturbed-law based indices

A PLI aims to measure the impact of the modification of an input density on a QoI, such as a quantile or a threshold exceedance probability of the model output [START_REF] Lemaître | Density modification based reliability sensitivity analysis[END_REF][START_REF] Sueur | Perturbed-Law based sensitivity indices for sensitivity analysis in structural reliability[END_REF]. In the following, we focus on a quantile of order α, often used as a risk measure in real-world applications [START_REF] Mousseau | Uncertainty quantification in a regulatory environment[END_REF][START_REF] Delage | Robustness analysis of epistemic uncertainties propagation studies in LOCA assessment thermal-hydraulic model[END_REF][START_REF] Larget | How to bring conservatism to a BEPU analysis[END_REF].

If we consider X to be a compact set of R d , then given the random vector X = (X 1 , . . . , X d ) ∈ X of our d independent imperfectly known input variables, G(•) our numerical model, and Y = G(X) ∈ R the model output, the quantile of order α of Y is given by:

q α = inf{t ∈ R, F Y (t) ≥ α} , (18)
where F Y is the cdf of the random variable Y . In order to compute the i-th PLI (see below),

we first change the density f i of X i into a density f iδ , where δ ∈ R + represents the level of perturbation. The perturbed quantile is then given by:

q α iδ = inf{t ∈ R, F Y,iδ (t) ≥ α} , (19)
where F Y,iδ is the cdf corresponding to the input variable X i sampled from f iδ . The i-th PLI S i is then simply defined as the relative change in the output quantile generated by the perturbation, i.e.,

S i (f iδ ) = q α iδ -q α q α . ( 20
)
This definition slightly differs from that proposed in previous studies [START_REF] Lemaître | Density modification based reliability sensitivity analysis[END_REF]Sueur et al., 2017), but is preferred here as it has been found to be more intuitive in a number of engineering studies that used PLI. Indeed, it simply consists of the relative variation of the quantile when the i-th input is submitted to a density perturbation and, as such, allows for a clearer interpretation.

In many applications including nuclear safety exercises, computer models are costly in terms of CPU time and memory. This means that only a limited number N of code runs are available for estimating all of the PLIs. Thus, we have a sample

Y N = {y (n) } 1≤n≤N of N outputs of the model from a sample X N = {X (n) = (x (n) 1 , ..., x (n) d )} 1≤n≤N of N independent
realizations of X. The estimation of the quantile is then based on the empirical quantile

estimator q α N = inf{t ∈ R, F N Y (t) ≤ α}, where F N Y (t) = 1 N N n=1 1 (y (n) ≤t)
is the empirical estimator of the cdf of Y . In order to estimate the perturbed quantile q α N,iδ from the same sample X N , we use reverse importance sampling [START_REF] Hesterberg | Estimates and confidence intervals for importance sampling sensitivity analysis[END_REF] (see Section 2 of the online supplementary materials) to compute F N Y,iδ [START_REF] Delage | Robustness analysis of epistemic uncertainties propagation studies in LOCA assessment thermal-hydraulic model[END_REF]:

F N Y,iδ (t) = N n=1 L (n) i 1 (y (n) ≤t) N n=1 L (n) i , (21) with L
(n) i the likelihood ratio f iδ (x (n) i ) f i (x (n) i )
. The estimator of the PLI is then

S N,iδ = q α N,iδ -q α N q α N . (22)
The theoretical properties of S N,iδ are studied in Section 3 of the online supplementary materials.

As presented in Section 3, the Fisher sphere of radius δ centered on the initial input distribution f i and written ∂B F (f i , δ) = {g, d F (f i , g) = δ}, provides a good base for perturbing distributions. We do not consider one specific perturbation at level δ, but instead a whole set of them: ∂B F (f i , δ). Over this set, we compute the maximum S + iδ and the minimum S - iδ of the PLI:

S + iδ = max g∈∂B F (f i ,δ) S i (g) , (23)
S - iδ = min g∈∂B F (f i ,δ) S i (g) , (24)
where S i (g) is the PLI with g the perturbed density for the variable X i .

Among all perturbed distributions at level δ, we look for those which most deviate the quantile from its original value. These two quantities S + iδ and S - iδ , called OF-PLI (for "Optimal Fisher-based PLI"), are measures of the robustness of the numerical code while taking into account uncertainty in the input distribution.

Practical implementation

The main problem when estimating the OF-PLI arises from the available sample size being finite, and possibly not as large as it needs to be. Indeed, due to inherent computational budget constraints in practice, there might not be enough sample points to correctly compute the perturbed quantile (and its confidence interval) at certain levels of perturbation.

Thus, the key issue is to determine the extent to which an input distribution should be perturbed. To address this, we propose adapting the empirical criterion from Iooss et al.

(2020) to help establish a maximal perturbed level δ max . For a proper estimation of OF-PLI indices, the number of points N Y in the output sample Y N exceeding the δ-perturbed quantile needs to be sufficiently high. As the sample size is finite, the δ-perturbed quantile can often not be computed accurately enough for high values of δ. In practice, a value of N Y = 10 has been decided on (after numerous numerical tests) as the smallest allowed for computing a PLI-quantile. As soon as a distribution on the Fisher sphere exceeds the criterion in [START_REF] Iooss | BEPU robustness analysis via perturbedlaw based sensitivity indices[END_REF], the corresponding radius is taken as δ max .

The estimation of the quantities of interest S + iδ and S - iδ is described in Algorithm 1.

We emphasize that this approach only applies to expensive computer models. Indeed, the bootstrap variance of the estimated quantile with reverse importance sampling tends to be very large, as illustrated in [START_REF] Iooss | BEPU robustness analysis via perturbedlaw based sensitivity indices[END_REF]. This is due to the likelihood ratio, which explodes locally. Thus, when dealing with a cheap code, one can directly resample over the perturbed distribution in order to estimate the output quantile. In this situation, there is no limiting level of perturbation δ max .

Algorithm 1 Estimation of S + iδ and S - iδ 1. Initialisation: δ (level of perturbation), i ∈ {1, . . . , d} (input number), a sample of K points on the Fisher sphere of radius δ centered at f i (using the numerical method of Section 2 of the online supplementary material).

2. For k = 1, . . . , K:

(a) estimate the α-quantile q α,(k) iδ

for the k-th perturbed density {f In future work, OF-PLI confidence intervals (computed via bootstrap) will provide valuable additional information like for instance confidence intervals, but further work on these is required and they are not included in the industrial application (section 5.2). The code for the OF-PLI is available at https://github.com/JeromeStenger/PLI-Technometrics.

The code for computing the older version of the PLI, called E-PLI (entropy-based PLI) in the following, is available in the sensitivity package of the software R.

Perturbed-law based indices in engineering studies

The PLI as defined here allows us to assess to what extent an output quantile can be affected by an error of magnitude δ in the characterization of an input distribution. In the following section, we compare, using a toy example, our new method (OF-PLI) to the earlier one (E-PLI). Following that, we illustrate the use of OF-PLI in a nuclear safety analysis of a pressurized water nuclear reactor. Moreover, as PLIs are based on changes in input distributions, they differ from global sensitivity measures [START_REF] Iooss | A review on global sensitivity analysis methods[END_REF], which evaluate the effect of input variability for fixed probabilistic models. To study potential coherence and/or divergence between the two approaches, we compare Sobol indices and OF-PLI on an analytic model in Section 4 of the online supplementary material.

A toy example: the Ishigami function

The Ishigami function [START_REF] Ishigami | An importance quantification technique in uncer-tainty analysis for computer models[END_REF] is often used as an example for uncertainty and sensitivity analysis methods, in particular because it exhibits strong nonlinearity and non-monotonicity. In this section, we apply the methodology introduced in Section 4.2 to estimate the OF-PLI and compare our results to E-PLI. The Ishigami function, which takes three input random variables (X 1 , X 2 , X 3 ), each distributed normally as N (0, 1), is defined with the following analytic formula:

G(x 1 , x 2 , x 3 ) = sin(x 1 ) + 7 sin(x 2 ) 2 + 0.1x 4 3 sin(x 1 ) . ( 25
)
We intend to evaluate the impact of a perturbed input distribution on the value of the 95%-quantile. In this simple example where the function is cheap to evaluate, we do not use the reverse importance sampling estimator of the quantile as proposed in Section 4.2.

Since the computational burden of running new code simulations is negligible, we instead draw new samples of size N = 2000 directly from the perturbed input distributions in order to compute the output perturbed quantile. Hence, this toy example does not follow the reverse important sampling estimation procedure detailed in Section 4.2. The OF-PLI are computed for perturbation levels δ from the interval [0, 0.9]. Notice that the choice δ max = 0.9 is arbitrary here since the reverse importance sampling estimator has been put aside. Indeed, here there is no actual limit for the maximal perturbation level as the OF-PLI are computed by resampling from the perturbed distribution. We chose a value of K = 100 trajectories over each Fisher sphere for computing the minimum and maximum of the OF-PLI. We also compute the 95%-confidence intervals calculated from 50 values of S + N,iδ and S - N,iδ .

The OF-PLI results are depicted in Figure 3. It appears that the third input has the most impact in shifting the quantile to the right. On the other hand, the second input has more impact in shifting the quantile to the left. These results are consistent with the well-known behavior of the Ishigami function in terms both of non-linearity of the model and primary influence of the third input.

Because the maximum and minimum of the OF-PLI are taken over the Fisher sphere, we depict in Figure 4 the distribution of the OF-PLI over the Fisher sphere with radius δ = 0.9 for the third input. One can see that in this case, the maximum and minimum are trajectories for δ varying in [0, 0.9], and their 95%-confidence intervals. For the Ishigami model, we directly resample from the perturbed distribution to compute the perturbed quantile instead of using the reverse important sampling estimator.

found respectively for high and low variance, with seemingly no change with respect to the mean.

These results can be compared to the E-PLI ones (see Section 2). As in this case, inputs follow normal distributions, applying an inverse cdf transformation is not necessary.

Therefore, perturbing the mean (resp. the variance) of the input variable is equivalent to drawing straight horizontal (resp. vertical) trajectories in the parameter space (see Figure 1 of the online supplementary material). Results are depicted in Figure 5; the mean of the Gaussian is perturbed in [-1, 1] and its variance in [0,4]. This corresponds to the range of variation of these parameters for the Fisher sphere radius varying in [0, 0.9]. We then compare the third input for the two methods. We saw in Figure 4 that the maximal OF-PLI was obtained with high variance and no mean perturbation, which is coherent with the results in Figure 5. However, we miss the true impact a perturbed density can induce in situations where the maximal and minimal OF-PLI are not located in these two axes, such as for instance is the case for the first variable. Hence, the E-PLI, restricted to two directions in the Fisher sphere, has limited interpretation.

A nuclear safety application

This industrial application concerns the study of the peak cladding temperature (PCT) of fuel rods in the case of a loss of coolant accident caused by an intermediate-sized break in the primary loop (IB-LOCA) in a nuclear pressurized water reactor. According to operating rules, this temperature must remain below a threshold to prevent any deterioration of the state of the reactor. The thermal-hydraulic transient caused by this accidental scenario is simulated with the CATHARE2 code [START_REF] Geffraye | CATHARE2 V2.5 2: A single version for various applications[END_REF], providing a temperature profile throughout time for the surface of the nuclear core assemblies [START_REF] Mazgaj | Comparison of CATHARE results with the experimental results of cold leg intermediate break LOCA obtained during ROSA-2/LSTF test 7[END_REF].

-1 0 1 µ -0.10 -0.05 0.00 0.05 E-PLI X 1 X 2 X 3 1 2 3 4 σ 2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 E-PLI X 1 X 2 X 3
The thermal-hydraulic model involves boundary and initial conditions, and many physical parameters (heat transfer coefficient, friction coefficient, etc.) whose exact values are unknown. The probability distributions of these inputs can be obtained from data and expert knowledge, or by solving inverse problems on experimental data [START_REF] Baccou | Development of good practice guidance for quantification of thermal-hydraulic code model input uncertainty[END_REF].

Input uncertainties are propagated inside this model and the goal of UQ is to estimate a high-order quantile of the PCT (model output). This α-quantile is interpreted as a pessimistic estimate of the PCT. Like any scientific approach, this methodology is based on hypotheses; regulatory authorities require an evaluation of the effect of these hypotheses on claimed results. Nuclear power operators are required to conduct studies in such a way to ensure that actual risks are overestimated. Under this "conservative principle" they are required to choose the most pessimistic assumption each time a modeling decision has to be made. In deterministic studies, this simply consists of taking the most penalizing value for each input variable. In this way, the resulting computation is seen as simulating the worst case scenario for the risk being examined. It is, however, not straightforward to implement such a principle when the numerical code is complex, with both non-monotonic effects of inputs and interactions between inputs. It is even harder to extend this rationale to a UQ framework that aims to represent all potential scenarios with related occurrence likelihoods. Nevertheless, recent works [START_REF] Larget | How to bring conservatism to a BEPU analysis[END_REF] have shown that the E-PLI can be useful as evidence in discussions on the choice of input distributions.

In the present application, we study a reduced-scale mock-up of a pressurized water reactor with the 7 imperfectly known inputs given in Table 1 [START_REF] Delage | Robustness analysis of epistemic uncertainties propagation studies in LOCA assessment thermal-hydraulic model[END_REF]. To compute the OF-PLI, an input-output sample of size N = 1000 is available, coming from a space filling design of experiments [START_REF] Fang | Design and modeling for computer experiments[END_REF] (whose points in [0, 1] d have been transformed to follow the input probability distributions). More precisely, a Latin hypercube sample minimizing the L 2 -centered discrepancy criterion [START_REF] Jin | An efficient algorithm for constructing optimal design of computer experiments[END_REF] has been used. The OF-PLI (with respect to a quantile of order α = 0.95) will then be estimated without any additional code run (see Section 4.1).

Figure 6 presents the maximum and minimum values of our two estimators S + N,iδ and S - N,iδ . We compute Fisher spheres with radius δ sampled uniformly in [0.1, 0.5], all respectively centered on the initial input distributions. The maximal perturbation level δ max = 0.5 is derived from the stopping criterion introduced in Section 4.2. On every sphere, K = 100 perturbed densities are sampled. The OF-PLIs are finally estimated on a dataset of size 1000. The stopping criterion of Section 4.2 gives a maximal admissible OF-PLI of 4%; this value is determined from the maximal admissible quantile for which there are N Y = 10 sample points above it. In fact, we see that S + N,7δ is close to this maximal admissible value.

Studies conducted earlier for the same application [START_REF] Delage | Robustness analysis of epistemic uncertainties propagation studies in LOCA assessment thermal-hydraulic model[END_REF] lead to similar results concerning the most influential inputs on the PCT's quantile: strong influence of variables 3 and 4, weak influence of variables 1, 2, and 5. With respect to these previous studies based on the standard space transformation, our information geometry perturbation method evaluates the influence of variable 7 as being smaller. In fact, as it is the only Gaussian distribution, the reverse transformation from the standard space to the physical one operates differently for this input than for the others. Overall, according to the values of S + N,3δ and S + N,7δ , the variables 3 and 7 appear to be the most influential inputs on the value of the PCT's quantile. This was not observed using the standard space transformation, and is probably due to the fact that the standard space approach allows for perturbing only one of the probability distribution parameters (for example the expected value) at a time.

-3 -2 -1 0 1 2 3 OF-PLI (%) S - N,1δ S - N,2δ S - N,3δ S - N,4δ S - N,5δ S - N,6δ S - N,7δ S + N,1δ S + N,2δ S + N,3δ S + N,4δ S + N,5δ S + N,6δ S + N,7δ
In contrast, our estimator corresponds to the maximal quantile deviation over an entire set of equivalent perturbations. This suggests two major advantages of our newly developed method: it (i) prevents interpretation bias induced by the standard space transformation, and (ii) allows for an exhaustive exploration of density perturbations for a given δ.

Conclusion

Based on the Fisher distance, we have defined an original methodology for perturbing input probability distributions in the specific case of mutually independent input random variables. The Fisher information is an intrinsic characteristic of probability measures and, in particular, does not depend on the specific choice of parametric representation. This fundamental property makes it the proper mathematical tool to compare perturbations on different imperfectly known physical inputs of a computer model, and also on different parameters of the same input distribution. It is even possible to remove all reference to a parametric sub-domain of the set of probability measures on X , given the nonparametric extension of the Fisher distance proposed by [START_REF] Holbrook | The nonparametric Fisher geometry and the chi-square process density prior[END_REF]. However, the latter is limited by practical issues as it supposes a finite-dimensional representation of the densities, for example by means of projection onto an orthonormal basis of the probability space. This implies truncating the infinite sum of the projections of a given probability on all elements of the basis. This approximation will be poor for probabilities which are very different from those of the chosen basis. This shows that in practice it is not easy to eliminate the reference to a particular parametric model, even in a nonparametric framework.

Nevertheless, based on PLIs, our method provides useful information on the most influential uncertainties related to the distributions of input variables, so-called "epistemic uncertainties". This is in particular crucial not only in making decisions concerning further research programs aiming at gaining better knowledge about these variables, but also to bring strong backing arguments to operators' safety demonstrations. Indeed, we argue that this methodology is adequate for uncertainty studies with unreliable input distribution identification, or when an improved level of robustness is demanded for the choice of input distributions. In our target application (nuclear licensing), the aim is not only to exhibit safety margin values for simulated accidents, but also to prove that the methodology as a whole does not induce any risk of underestimating these values. Hence we are not only looking for a worst case assessment method, but also for a more global understanding of how a potential error in an input's distribution affects the output. In this sense, a practical option to increase conservatism in UQ studies is to replace one or several input distributions by penalized deterministic values, or by penalized versions of the distributions themselves.

This nevertheless implies justifying the choice of the variables for which this penalization is done (see, e.g., [START_REF] Larget | Increasing conservatism in BEPU IB LOCA safety studies using complementary and industrially cost effective statistical tools[END_REF]).

Further investigations remain to be completed, due in part to this method increasing the numerical complexity and computational time (by several hours) required compared to the previous method of [START_REF] Lemaître | Density modification based reliability sensitivity analysis[END_REF]. Several Monte Carlo loops are needed to compute the maximal and minimal PLI over Fisher spheres and we have ongoing work on the improvement of the estimation of the maximum and the minimum of the PLI on a Fisher sphere. In particular, using a better-adapted programming language, a large computational gain (of hours) could be expected. There is a known numerical issue with the reverse importance sampling strategy whereby the likelihood ratio tends to explode, as do the confidence intervals. Moreover, the method consists in sampling trajectories over the Fisher sphere, but one could benefit from a more advanced strategy by directly optimizing the PLI over the sphere via, for instance, gradient descent along this manifold.

The crucial problem of probabilistic dependency between inputs needs also to be studied in order to extend our framework to the non-independent input case; works in robustness analysis dealing with dependent inputs can be found for instance in [START_REF] Pesenti | Cascade Sensitivity Measures[END_REF].

Moreover, using a distance based on the Fisher information might be less intuitive for non-statisticians than simply shifting moments. In order to promote our methodology in engineering applications, it is essential to explicit its link with statistical test theory, thus helping practitioners interpreting the results it provides. Finally, future works will improve the estimation accuracy of the OF-PLI, especially for high values of δ where numerical difficulties occurred, as mentioned in Section 1 of the online supplementary material. 7. The relative error for the Adams-Moulton method is negligible, while the maximum relative error for the explicit Euler scheme is around 0.3%. Hence, in the Gaussian case the Adams-Moulton scheme is preferable. However, instabilities have been observed in practice. Symplectic methods [START_REF] Amari | Methods of Information Geometry[END_REF][START_REF] Leimkuhler | Simulating Hamiltonian Dynamics[END_REF], and in particular a symplectic Euler algorithm, could help to alleviate this problem by forcing the Hamiltonian to be constant. This will be the subject of future work. Note also that truncation can lead to other numerical errors when the radius δ is too large. Indeed, the normalization constant of some truncated distributions can sometimes become smaller than a computer's machine precision. 8 The reverse importance sampling method Importance sampling (IS) is a Monte Carlo technique to compute statistical quantities such as E[φ(X)] (X = (X 1 , ..., X d )) with a sample XN = { X(1) , . . . , X(N) } following a different distribution from that of X. If we call π the density of X and π the sampling one, we have

E π [φ(X)] = φ(x)π(x)dx = φ(x)π(x) π(x) π(x) dx = E π[φ(X )L(X)],
where L = π/π is the likelihood ratio allowing us to convert an integral under the distribution of X into an integral under the sampling one. This leads to estimate of

M = E[φ(X)] via: MIS = 1 N N n=1 φ( X(n) )L( X(n) ), (27)
with an estimation error given by Var

[ MIS ] = N -1 Var π[φL].
The use of an alternative density (called an instrumental density) is most often a vari-ance reduction method [START_REF] Kahn | Methods of Reducing Sample Size in Monte Carlo Computations[END_REF], as the variance of M can be much lower than that of the usual empirical mean for a well chosen π. Indeed, an optimal variance reduction can be obtained for a sampling distribution ∝ πφ, which however implies that we know the value of E[φ] we are aiming to estimate.

The principle is similar if we aim to estimate quantiles corresponding to a perturbed density, though from a sample generated from the initial density of X [START_REF] Hesterberg | Estimates and confidence intervals for importance sampling sensitivity analysis[END_REF].

Here the IS is called "reverse IS" as we do not change the sampling density to achieve a better estimation of the initial quantity, but instead change the quantity to be estimated with the same initial sample. Our purpose is to estimate a perturbed quantile of Y , e.g., a quantile of the distribution of Y iδ = G(X 1 , . . . , X i-1 , Xi , X i+1 , . . . , X d ), where the i-th input Xi is perturbed and has a density of f iδ instead of the initial one.

Hence the initial density of X can be seen as the instrumental density of a typical IS-Monte Carlo scheme, and the perturbed density is the one under which we must estimate the cdf F Y,iδ . If we denote π(x 1 , . . . , x d ) = Π d i=1 f i (x i ) the initial density and π(x 1 , . . . , x d ) =

f 1 (x 1 ) × ... × f i-1 (x i-1 ) × f iδ (x i ) × f i+1 (x i+1 ) × . . . × f d (x d
) the perturbed one as regards the i-th input, and recalling that F Y,iδ (y) = E π[1 (G(X)≤y) ], we have the following estimator for the perturbed cdf:

F N Y,iδ (y) = 1 N N n=1 L (n) i 1 (y (n) ≤y) . (28)
In this estimator, illustrated in Figure 9, we recognize that of formula (27) with φ(•) = 1 (G(•)≤y) and L(x

(n) 1 , . . . , x (n) d ) = π(x (n) 1 , . . . , x (n) d )/π(x (n) 1 , . . . , x (n) d ) = f i (x (n) i )/f iδ (x (n) i ) = L (n) i .
The estimator q α N,iδ of the α-order quantile of Y iδ is then obtained by plugging F N Y,iδ into the definition of the theoretical quantile q α iδ = inf{t ∈ R, F Y,iδ (t) ≥ α}, which gives:

q α N,iδ = inf{t ∈ R, F N Y,iδ (t) ≥ α}.
Note that the estimator (28) of F Y,iδ is normalized by 1/N , which slightly differs from the one presented here, which is normalized by 1/ N L (n)

i , as well as the one plotted in Figure 9. This guarantees F N Y,iδ (y) → 1 when y → +∞ for finite values of N without changing the asymptotic properties of the quantile estimator [START_REF] Glynn | Importance Sampling for Monte Carlo Estimation of Quantiles[END_REF].

9 Theoretical properties of the PLI-quantile estimator

In this section, we investigate some theoretical aspects of the PLI estimator S N,iδ . As it is based on quantile estimators, we first focus on the asymptotic properties of the estimator

q α N , q α N,iδ . Theorem 1. Suppose that F Y is differentiable at q α = F -1 Y (α) with F Y (q α ) > 0 and that F Y,iδ is differentiable at q α iδ = F -1 Y,iδ (
α) with F Y,iδ (q α iδ ) > 0. We define:

σ 2 i = α(1 -α) F Y (q α ) 2 , ( 29
) σ2 iδ = E f iδ (X i ) f i (X i ) 2 (1 (G(X)≤q α iδ ) -α) 2 F Y,iδ (q α iδ ) 2 , θi = E f iδ (X i ) f i (X i ) 1 (G(X)≤q α ) 1 (G(X)≤q α iδ ) -αE[1 (G(X)≤q α iδ ) ] F Y (q α )F Y,iδ (q α iδ )
.

Let us write the covariance matrix Σ as:

Σ =     σ 2 θi θi σ2 iδ     ,
and suppose that it is invertible and E

f iδ (X i ) f i (X i ) 3 < +∞. Then √ N         q α N q α N,iδ     -     q α q α iδ         L -→ N (0, Σ) .
Proof. We study the consistency and asymptotic normality of specific M and Z-estimators in order to establish the proof. We suppose the theory around these estimators is known so that details can be kept to the bare minimum. Further details can be found in Chapters 5.2 and 5.3 of Van der Vaart (2000). Given a sample (X (n) ) n∈(1,...,N ) where X is a d-dimensional random vector, we define

η = α 1 -α , m θ (x) = -(G(x) -θ)1 (G(x)≤θ) + η(G(x) -θ)1 (G(x)>θ) , M N (θ 1 , θ 2 ) = 1 N N n=1 m θ 1 (X (n) ) + f iδ (X (n) i ) f i (X (n) i ) m θ 2 (X (n) ) , θN = arg max M N (θ 1 , θ 2 ) . ( 30
)
θN is defined such that its two components correspond respectively to the estimators qα N and qα N,iδ of the quantile and the perturbed quantile. The map θ → ∇ θ M N (θ) with θ = (θ 1 , θ 2 ) T has two non-decreasing components (it is a sum of non-decreasing maps). Now, by definition of θN and concavity of M n (θ), it holds that ∇ θ M N ( θN ) = 0. Furthermore, we have that

∇ θ M N (θ) P -→ [(1 + η)F Y (θ 1 ) -η, ((1 + η)F Y,iδ (θ 2 ) -LN η)] T with LN = 1 N N n=1 f iδ (X (n) i ) f i (X (n) i )
, and this limit is a strictly non-decreasing function. Therefore, the assumptions of Lemma 5.10 in (Van der Vaart, 2000, p.47) are satisfied, proving the consistency of the estimator θN P -→ (q α , q α iδ ) T .

The asymptotic normality is studied via the map mθ (

x) → m θ 1 (x)+ f iδ (x) f i (x) m θ 2 (x) which is Lipschitz for the variable θ with Lipschitz constant h(x) = max(1, η) 1 + f iδ (x i ) f i (x i ) . The function h belongs in L 2 if E f iδ (X i ) f i (X i ) 2 < +∞.
The map mθ is also differentiable at )] with gradient:

θ 0 = arg max θ∈Θ E[ mθ (X
∇ θ 0 mθ 0 (x) = ((1 + η)1 (G(x)≤θ 1 ) -η, f iδ (x i ) f i (x i ) ((1 + η)1 (G(x)≤θ 2 ) -η)) T . (31)
Moreover, the map θ → E[ mθ (X)] has the following Hessian:

V θ 0 =     (1 + η)F Y (q α ) 0 0 (1 + η)F Y,iδ (q α iδ )     , (32)
which is non-negative definite symmetric whenever F Y (q α ) > 0 and F Y,iδ (q α iδ ) > 0. Hence, Theorem 5.23 in (Van der Vaart, 2000, p.53) applies. This proves the asymptotic normality of the estimator (q α , qα iδ ) T .

The PLI S iδ is a straightforward transformation of (q α , q α iδ ) T . To obtain the almost sure convergence of S N,iδ to S iδ , it suffices to apply the continuous mapping theorem to the function s(x, y) = yx x .

Theorem 2. Under the assumptions of theorem 1, we have Notice that the asymptotic variance relies on the initial α-quantile and the perturbed quantile, which are precisely what we want to estimate. Hence, Theorem 2 cannot be used for building asymptotic confidence intervals. However, its convergence properties are important for our method's credibility and acceptance. In practice, the estimation error can be measured using bootstrapping [START_REF] Efron | Bootstrap Methods: Another Look at the Jackknife[END_REF].

10 Applying PLI to an analytic flood risk model

The model of interest concerns a flooded river simulation, which is especially useful in assessing the risk of submergence of a dike protecting industrial sites near a river. To this total Sobol indices of the inputs of the flood model (Eq. ( 34)), the asymptotically efficient pick-freeze estimator [START_REF] Prieur | Variance-Based Sensitivity Analysis: Theory and Estimation Algorithms[END_REF] is used with an elementary Monte Carlo matrix of size 10 6 . This gives a total cost of N = 6 × 10 6 model runs and a standard deviation of the indices' estimation error smaller than 10 -3 . As shown in Table 3, in the central part of the distribution (conventional Sobol indices), we observe that the variable Q is clearly more influential than the variable K s , whereas Z v and Z m appear to have almost no influence on the output. From the target Sobol indices, we observe that in the extreme part of the distribution (close to the 95%-quantile), Q and K s have the same total effect (due to a strong interaction effect between them, which comes from the fact that threshold exceedence is principally sensitive to the joint values of these two inputs). indicate that the maximal values of the OF-PLI are obtained for densities allocating the largest probabilities to large input values (blue curves), whereas the minimal values are obtained for densities with the least weighted tails (green curves). This behavior was predictable here since the river height is an increasing function of river flow (see Eq. ( 34)), though this type of analysis can in practice provide substantial information in real-world engineering studies. Finally, we also observe that Fisher spheres sometimes have nontrivial shapes, as illustrated in Fig. 10a. In this plot, we see that Fisher spheres are non-closed curves for radius 0.3 ≤ δ, the parameter's domain being bounded by 0 for β and γ. Such behaviour is specific to each parametric family of densities and, for instance, is not observed for non-truncated normal distributions.

The results of the OF-PLI, displayed in Figure 11, confirm those of the target Sobol indices (see Table 3): the variables 3 and 4, corresponding to Z v and Z m , are much less influential on the output quantile of level α = 0.95 than the variables 1 and 2, corresponding to Q and K s . It turns out that in this case study, perturbations of Q and K s appear to have comparable effects on the 95%-quantile of H, although they contribute quite differently to the output variance. On the other hand, compared to target Sobol indices, OF-PLI provide

Figure 2

 2 Figure 2 shows the position of four Gaussian distributions in the ( µ √ 2 , σ) half-plane. It is clear that the distributions C and D are more difficult to distinguish than A and B, though

Figure 2 :

 2 Figure 2: Four Gaussian distributions represented in the parameter space (left) and their

  on the sample X N , (b) check whether the number of points in the output sample below or above the perturbed quantile q α,(k) iδ satisfies the stopping criterion N Y ≥ 10, (c) compute the PLI estimator S (k) N,iδ . 3. the estimators S + N,iδ and S - N,iδ of the quantities of interest S + iδ and S - iδ are taken as the maximal and minimal values of the PLI sampled on the Fisher sphere { S

Figure 3 :

 3 Figure 3: Minimum and maximum of the OF-PLI over the Fisher sphere with K = 100

Figure 4 :

 4 Figure 4: Value of the OF-PLI ŜN,iδ (red line) for the third input of the Ishigami model

Figure 5 :

 5 Figure 5: Computation of the E-PLI. Left: perturbation of the mean of the Gaussian

Figure 6 :

 6 Figure 6: Bootstrap mean of the maximum and minimum of the OF-PLI S iδ for the

Figure 7 :

 7 Figure 7: Geodesics in Gaussian information geometry computed with Euler's explicit

Figure 8 :

 8 Figure 8: Relative variation of the Hamiltonian ∆ along a geodesic for two different nu-

Figure 9 :

 9 Figure 9: Empirical (solid lines) and theoretical (dotted lines) cumulative distribution

√

  N ( S N,iδ -S iδ ) L -→ N (0, d T s Σd s ), with d s =

First

  OF-PLI (w.r.t. a quantile of order α = 0.95) for the flood model inputs by increasing the Fisher spheres radii δ ∈ [0, 1.4] by steps of 0.1. The spheres are respectively centered on the distributions shown in Table2. On each of these spheres, we compute the OF-PLI for K = 100 different perturbed distributions using a sample of N = 2000 points distributed according to the initial distribution. The maximal radius δ max = 1.4 comes from the stopping criteria explained in Section 4.2 ("Practical implementation") of the paper. More precisely, the criterion is hit for the first input Q at perturbation level δ > 1.4, meaning there are less than N Y = 10 sample points above the maximal perturbed quantile.

Figure 10

 10 Figure 10 depicts how the Fisher sphere centered at the variable Q deforms, and how

Figure 11 :

 11 Figure 11: Maximum and minimum estimated value of the OF-PLI S + N,iδ and S - N,iδ for the

Table 1 :

 1 Input variables of the CATHARE2 code with their associated probability distri-

	butions.
	Variable	Input
			Probability distribution
	number	name
	1	STMFSCO	Uniform U([-44.9, 63.5])
	2	STBAEBU	Truncated Log Normal LN (0, 0.76) on [0.1, 10]
	3	STOIBC1	Truncated Log Normal LN (0, 0.76) on [0.1, 10]
	4	STOIBC3	Truncated Log Normal LN (0, 0.76) on [0.1, 10]
	5	STOIDC	Truncated Log Normal LN (0, 0.76) on [0.1, 10]
	6	STOICO	Truncated Log Normal LN (-0.1, 0.45) on [0.23, 3.45]
	7	CLFBR	Truncated Normal N (6.4, 4.27) on [0, 12.8]

Table 3 :

 3 Sobol index estimates of the flood model inputs.

	Inputs	Q	K s	Z v	Z m

The triangular distribution T (-1, 0, 1) is parametrized by its minimum a, mode b, and maximum c.

Online Supplementary Materials for the Paper "An Information Geometry Approach to Robustness Analysis for the Uncertainty Quantification of Computer Codes"

7 Computing Fisher spheres: Numerical results

Hamilton's equations (see Section 3.2 of the main article),

can be solved using numerical approximation methods. Figure 7 illustrates our numerical solutions in the Gaussian case, i.e., when

We compare the solution given by the explicit Euler algorithm to that given by the Adams-Moulton algorithm. We recall that in the Gaussian case we have at our disposal an exact analytical expression of the Fisher sphere, as detailed in [START_REF] Costa | Fisher information distance: a geometrical reading[END_REF]. The Fisher sphere is centered at N (0, 1) with radius δ = 1. Notice that there is no observable difference between the two methods in Figure 7. Hence, the numerical error is estimated from the variation of the Hamiltonian value that should be theoretically conserved along the geodesics. As a consequence, it is possible to quantify the performance of the numerical approximation by

Here, ∆ represents the relative variation of the Hamiltonian along the path q computed with our numerical methods. Venant equations. This model computes H, the maximal annual water level of the river, from four parameters Q, K s , Z m and Z v , which are considered imperfectly known:

The inputs are modeled as random variables with the associated truncated distributions given in Table 2 [START_REF] Iooss | A review on global sensitivity analysis methods[END_REF]. In global sensitivity analyses, Sobol indices are the most popular sensitivity measures because they are easy to interpret; each Sobol index represents a share of the output variance, and all indices sum to 1 under the assumption of independent inputs [START_REF] Sobol | Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates[END_REF][START_REF] Saltelli | On the relative importance of input factors in mathematical models: Safety assessment for nuclear waste disposal[END_REF][START_REF] Prieur | Variance-Based Sensitivity Analysis: Theory and Estimation Algorithms[END_REF]. We will compare these to the results of our robustness analysis framework in order to illustrate differences. Note however that these conventional Sobol indices focus on the central part of the distribution (variance of the output). We then also compute the target Sobol indices [START_REF] Marrel | Statistical developments for target and conditional sensitivity analysis: application on safety studies for nuclear reactor[END_REF], i.e., Sobol indices applied to the indicator function of exceeding a given threshold (chosen here as the 95%-quantile of the output). To compute the first order and