Vincent Moutoussamy
email: vincent.moutoussamy@edf.fr

Simon Nanty
email: simon.nanty@cea.fr

Benoît Pauwels
email: benoit.pauwels@ifpen.fr

EMULATORS FOR STOCHASTIC SIMULATION CODES

Numerical simulation codes are very common tools to study complex phenomena, but they are often time-consuming and considered as black boxes. For some statistical studies (e.g. asset management, sensitivity analysis) or optimization problems (e.g. tuning of a molecular model), a high number of runs of such codes is needed. Therefore it is more convenient to build a fast-running approximation -or metamodel -of this code based on a design of experiments. The topic of this paper is the definition of metamodels for stochastic codes. Contrary to deterministic codes, stochastic codes can give different results when they are called several times with the same input. In this paper, two approaches are proposed to build a metamodel of the probability density function of a stochastic code output. The first one is based on kernel regression and the second one consists in decomposing the output density on a basis of well-chosen probability density functions, with a metamodel linking the coefficients and the input parameters. For the second approach, two types of decomposition are proposed, but no metamodel has been designed for the coefficients yet. This is a topic of future research. These methods are applied to two analytical models and three industrial cases.

Résumé. Les codes de simulation numérique sont couramment utilisés pour étudier des phénomènes complexes. Ils sont cependant souvent coûteux en temps de calcul. Pour certaines études statistiques (e.g. gestion d'actifs, analyse de sensibilité) ou problèmes d'optimisation (e.g. réglage des paramètres d'un modèle moléculaire), un grand nombre d'appels à de tels codes est nécessaire. Il est alors plus approprié d'élaborer une approximation -ou métamodèle -peu coûteuse en temps de calcul à partir d'un plan d'expérience. Le sujet de cet article est la définition de métamodèles pour des codes stochastiques. Contrairement aux codes déterministes, les codes stochastiques peuvent renvoyer des résultats différents lorsqu'ils sont appelés plusieurs fois avec le même jeu de paramètres. Dans cet article deux approches sont proposées pour la construction d'un métamodèle de la densité de probabilité de la sortie d'un code stochastique. La première repose sur la régression à noyau et la seconde consiste à décomposer la densité de la sortie du code dans une base de densités de probabilité bien choisies, avec un métamodèle exprimant les coefficients en fonction des paramètres d'entrée. Pour la seconde approche, deux types de décomposition sont proposés, mais aucun métamodèle n'a encore été élaboré pour les coefficients: c'est un sujet de recherches futures. Ces méthodes sont appliquées à deux cas analytiques et trois cas industriels.

Introduction

Numerical simulation codes are very common tools to study complex phenomena. However these computer codes can be intricate and time-consuming: a single run may take from a few hours to a few days. Therefore, when carrying out statistical studies (e.g. sensitivity analysis, reliability methods) or optimization processes (e.g. tuning of a molecular model), it is more convenient to build an approximation -or metamodel -of the computer code. Since this metamodel is supposed to take less time to run than the actual computer code, it is substituted for the numerical code. A numerical code is said deterministic if it gives the same output each time it is called with a given set of input parameters. For this case, the design of metamodels has already been widely studied. The present work deals with the formulation of a metamodel for a numerical code which is not deterministic, but stochastic. Two types of numerical codes can be covered by this framework. First, the code under consideration may be intrinsically stochastic: if the code is called with the same input several times, the output may differ. Thus, the output conditionally to the input is itself a random variable. Second, the computer code may be deterministic with complex inputs. For example, these inputs may be functions or spatiotemporal fields. These complex inputs are not treated explicitly and are called uncontrollable parameters, while the other parameters are considered as controllable parameters. When the code is run several times with the same set of controllable parameters, the uncontrollable parameters may change, and thus different output values may be obtained: the output is a random variable conditionally to the controllable parameters.

Two types of approaches can be found in the literature for stochastic code metamodeling. In the first one, strong hypotheses are made on the shape of the conditional distribution of the output. In [START_REF] Rigby | Construction of reference centiles using mean and dispersion additive models[END_REF] the authors consider that the output is normally distributed and estimate its mean and variance, in [START_REF] Reich | Variable selection for high dimensional bayesian density estimation: application to human exposure simulation[END_REF] the authors assume that the distribution of the output of the code is a mixture of normal distributions and estimate its mean. In the second approach, the mean and variance of the output are fitted jointly. In [START_REF] Iooss | Global sensitivity analysis of computer models with functional inputs[END_REF][START_REF] Zabalza-Mezghani | Uncertainty management: From geological scenarios to production scheme optimization[END_REF] the authors use joint generalized linear models and joint generalized additive models respectively to fit these quantities. Joint metamodels based on Gaussian processes have been developed in [START_REF] Ankenman | Stochastic kriging for simulation metamodeling[END_REF][START_REF] Marrel | Global sensitivity analysis of stochastic computer models with joint metamodels[END_REF]. As we can observe, so far in the literature the metamodels proposed only provide an estimate for one or two of the output distribution moments, sometimes under strong hypotheses on the distribution of the output.

Our approach does not require such a priori hypotheses. Furthermore we study the whole distribution of the output, allowing to derive estimations for the mean, variance and other moments, and quantiles as well. The metamodeling is not carried out on the stochastic computer code under consideration itself, but rather on the code G that has same input and that returns the probability density function of the computer code output rather than just a single random value:

G : X ⊂ R d → F x → G(x, •),
where F is the space of probability density functions with support in an interval I of R:

F =    f : I → R + : I f = 1    .
Let X N = {(x i , f i ) : i = 1, . . . , N } ⊂ (X × F) N be a training set (N ∈ N \ {0}): we suppose that we have access to the probability density functions f i of the outputs associated with N given sets of input parameters x i , i.e.

f i = G(x i , •) (i = 1, . . . , N ).
In practice, for a given i in {1, . . . , N }, the output of a stochastic computer code is not a probability density function f i . We only have access to a finite sample of output values y i1 , . . . , y iM which can be written as a vector y i1 • • • y iM representing the function. In real applications (cf. Sections 3.1, 3.2 and 3.3) we chose to derive each probability density function f i thanks to kernel density estimation. In this non-parametric method developed by Rosenblatt [START_REF] Rosenblatt | Remarks on some nonparametric estimates of a density function[END_REF] and Parzen [START_REF] Parzen | On estimation of a probability density function and mode[END_REF], the estimator of the probability density function f i , i = 1, . . . , N at a point y in I is defined as follows:

fi (y) = 1 M M j=1 K H (y ij -y) ,
where H ∈ R is the bandwidth and K H is the kernel. Let x 0 ∈ X be a new set of input parameters, we want to estimate the probability density function f 0 of the corresponding output:

f 0 = G(x 0 , •) : I ⊂ R → R + t → f 0 (t) = G(x 0 , t).
For example, let us suppose that X = R × R * + , x 0 = (µ 0 , σ 0 ) and G((µ 0 , σ 0 ),

•) = e -1 2
•-µ 0 σ 0 2 2πσ 2 0
. Figure 1 shows the output of G for different values of µ 0 and σ 0 . Section 1 expounds two types of metamodeling for stochastic computer codes: the first one provides a metamodel with functional outputs based on kernel regression, the second one consists in choosing relevant probability density functions to express the probability density of the output as a convex combination of them, the coefficients being linked to the input parameters by another metamodel. In Section 2, the previous methods are applied to two analytical test cases. Finally, the same methods are carried out on three industrial applications in Section 3.

Metamodels for probability density functions

In this section we expound our two approaches of formulating a metamodel for G. The first relies on kernel regression (cf. Section 1.1). The second consists in writing f 0 as a convex combination of well-chosen probability density functions (cf. Section 1.2). In the following we denote • L 2 and •, • L 2 the norm and inner product of L 2 (I) respectively: for all u and v in L 2 (I),

u L 2 =   I u 2 (t) dt   1/2
and u, v L 2 = I u(t)v(t) dt.

Kernel regression-based metamodel

In this subsection, first we apply the classical kernel regression method to the problem under consideration, then a new kernel regression estimator involving the Hellinger distance is introduced and, finally, some perspectives of improvement are suggested.

Classical kernel regression

Let x 0 ∈ X and f 0 = G(x 0 , •). We suppose that the sample set X N = {(x i , f i ) : i = 1, . . . , N } ⊂ (X × F) N is available. We want to estimate f 0 by f0 ∈ F, knowing X N . However, forcing f0 to be in F can be difficult in practice. We propose a first estimator given by f0 = N i=1 α i f i where α i ∈ R (i = 1, . . . , N ). In order to have f0 in F, we impose α i ≥ 0 (i = 1, . . . , N ) ,

N i=1 α i = 1.
The non-parametric kernel regression verifies the previous constraints. This estimator, introduced in [START_REF] Elizbar | On estimating regression[END_REF][START_REF] Watson | Smooth regression analysis[END_REF], can be written as follows:

f0 = N i=1 K H (x i , x 0 ) N j=1 K H (x j , x 0 ) f i , (1)
where K H : R d × R d → R is a kernel function. In the following, the Gaussian kernel is used:

K H (x, y) = 1 2π det(H) e -(x-y) T H -1 (x-y) x, y ∈ R d ,
where

H = diag(h 1 , . . . , h d )
is the (diagonal) bandwidth matrix, with h 1 , . . . , h d ∈ R + . The higher h j is, the more points x i are taken into account in the regression in the direction j. The method relies on the intuition that if x 0 is close to x i then f 0 will be more influenced by the probability density function f i (see Figure 2). We want that the bandwidth matrix H * minimizes the global error of estimation. In order to get an approximation of H * , one uses the leave-one-out cross-validation, as proposed in Hardle and Marron [START_REF] Hardle | Optimal bandwidth selection in nonparametric regression function estimation[END_REF]. Then

H * ∈ arg min h1,...,h d ∈R+ N i=1 f-i -f i 2 L 2 , with f-i = N j=1 j =i α -i j f j α -i j = e -(xi-xj ) t H -1 (xi-xj ) N l=1 l =i e -(xi-x l ) t H -1 (xi-x l )
.

This optimization problem is solved with the L-BFGS-B algorithm developed by Byrd et al. [START_REF] Richard H Byrd | A limited memory algorithm for bound constrained optimization[END_REF]. Finally, the kernel regression metamodel is given by f0

= N i=1 α i,H * f i , (2)
with

α i,H * = e -(x0-xi) T (H * ) -1 (x0-xi) N j=1 e -(x0-xj ) T (H * ) -1 (x0-xj )
.

Kernel regression with the Hellinger distance

The kernel regression introduced in equation ( 1) can be considered as a minimization problem. The estimator of the function f 0 corresponding to the point x 0 is as follows:

f0 = arg min f ∈F N i=1 K H (x i , x 0 ) I (f i -f ) 2 . ( 3
)
The kernel regression is therefore a locally weighted constant regression, which is fitted by minimizing weighted least squares. The equivalence between estimators (1) and ( 3) is proven in Appendix A.1. We observe that the L 2 distances between the sample functions and the unknown function appear in the objective function of the minimization problem (3).

Another classical example of distance between probability density functions is the Hellinger distance. We derived a new kernel regression estimator with respect to this distance replacing the L 2 distance by the Hellinger distance in equation [START_REF] Ankenman | Stochastic kriging for simulation metamodeling[END_REF]. The kernel estimator thus takes the following form:

f0 = arg min f ∈F N i=1 K H (x i , x 0 ) I f i (t) -f (t) 2 dt. ( 4
)
The following analytical expression can be derived for this Hellinger distance-based estimator:

f0 = N i=1 K H (x i , x 0 ) f i 2 I N i=1 K H (x i , x 0 ) f i 2 . (5)
The calculation that leads to ( 5) is provided in Appendix A.2.

From now on, we will denote f0,L 2 and f0,He the estimators given by equations ( 2) and ( 5) respectively. Let us illustrate the Hellinger kernel regression on the four gaussian density functions given on Figure 2. On Figure 3 the same N = 4 learning curves (in dark blue, red, light blue and green) are represented, along with an unknown output f 0 (in dashed black) and its associated Hellinger kernel estimator f0,He (in orange). We can observe that, although the output is not very well estimated, most of the weight is given to the curve corresponding to the closest vector of parameters (in green).

Perspectives

One of the well-known major drawbacks of the kernel regression is due to the curse of dimensionality. The bandwidth estimation quality is deteriorating as the dimension d of the problem increases for a constant sample size. First, as the sample points are sparsely distributed in high dimension, the local averaging performs poorly. Second, the optimization problem to find the optimal bandwidth is more complex in higher dimension. In the presented work, kernel regression has been used with a fixed bandwidth for all sample points. The Adaptive or "Variable-bandwidth" kernel estimation techniques have been proposed in the literature to use varying bandwidth over the domain of the parameters X . See Muller and Stadtmuller [START_REF] Muller | Variable bandwidth kernel estimators of regression curves[END_REF] and Terrell and Scott [START_REF] George | Variable kernel density estimation[END_REF] for more details on these methods. Terrell and Scott [START_REF] George | Variable kernel density estimation[END_REF] review several adaptive methods and compare their asymptotic mean squared error, while Muller and Stadtmuller [START_REF] Muller | Variable bandwidth kernel estimators of regression curves[END_REF] propose a varying bandwidth selection method. In numerous cases, using non-constant bandwidths leads to better estimates as it enables to better capture the different features of the curves. Hence, the quality of the proposed methods could be improved by using adaptive kernel regression. However, the use of kernel regression is not advised in high dimension.

Functional decomposition-based metamodel

In this section, we aim at building a set of basis functions w 1 , . . . , w q in F to approximate f 0 by an estimator of the form

f0 = q k=1 ψ k (x 0 )w k , (6)
where ψ k are functions from X to R (k = 1, . . . , q), called coefficient functions, satisfying

     ψ k (x) ≥ 0 q k=1 ψ k (x) = 1 (x ∈ X ) .
Kernel regression and functional decomposition have different goals but have similar forms. The former provides an estimator which is a convex combination of all the N sampled probability density functions f 1 , . . . , f N . The latter gives a combination of fewer density functions built from the experimental data, thus lying in a smaller space, in order to reduce the dimension of the problem. We observe that the estimator provided by classical kernel regression (1) is of the form (6) with q = N (no dimension reduction), w k = f k for k = 1, . . . , q and

ψ k (x 0 ) = K H (x k , x 0 ) N j=1 K H (x j , x 0 ) (k = 1, . . . , q) .
We will build the basis functions and the coefficient functions so that they fit the available data X N = {(x i , f i ) : i = 1, . . . , N }. Thus the functions w 1 , . . . , w q will be designed so that, for all i in {1, . . . , N }, f i is approximated by

fi = q k=1 ψ ik w k where            ψ ik = ψ k (x i ) (k = 1, . . . , q) ψ ik ≥ 0 (k = 1, . . . , q) q k=1 ψ ik = 1. (7)
For all i in {1, . . . , N }, we require the approximation fi of f i to be a probability density function as well; therefore fi must be non-negative and have its integral equal to 1.

We study three ways of choosing the basis functions. In Section 1.2.1 we discuss the adaptation of Functional Principal Components Analysis (FPCA) to compute the basis functions: w 1 , . . . , w q are then orthonormal with integral equal to 1 and their negative values are interpreted as 0, hence the predictions on the experimental design are probability density functions. In the following two subsections we propose two new approaches to probability density function decomposition. The non-negativity constraint is taken into account directly rather than through a post-treatment of FPCA results (setting negative values to zero and renormalizing). In Section 1.2.2, we propose to adapt the Empirical Interpolation Method or Magic Points method [START_REF] Maday | A general, multipurpose interpolation procedure: the magic points[END_REF] to our problem: we obtain an algorithm to select the basis functions (not necessarily orthogonal) among the sample distributions f 1 , . . . , f N (while losing the interpolation property). It follows that the predictions on the experimental design are probability density functions. In Section 1.2.3 a third approach is proposed, consisting in computing directly the basis functions that minimize the L 2 -approximation error on the experimental design (without imposing orthogonality):

N i=1 f i - q k=1 ψ ik w k 2 L 2 . ( 8
)
We introduce an iterative procedure to ease the numerical computation. The basis functions obtained are probability density functions and so are the predictions on the experimental design. Finally, in Section 1.2.4 we discuss the formulation of the coefficient functions.

Constrained Principal Component Analysis (CPCA)

Among all the dimension reduction techniques, one of the best known is the Functional Principal Component Analysis (FPCA), developed by Ramsay and Silverman [START_REF] Ramsay | Functional Data Analysis[END_REF], which is based on Principal Component Analysis (PCA). Let us denote f the mean estimator of f 1 , . . . , f N :

f = 1 N N i=1 f i .
The goal of FPCA is to build an orthonormal basis w 1 , . . . , w q for the centered functions f c i = f i -f (i = 1, . . . , N ) so that the projected variance of every centered function onto the functions w 1 , . . . , w q is maximized. The following maximization problem is thus considered: Minimize w1,...,wq:I→R N i=1

f c i - q k=1 f c i , w k L 2 w k 2 L 2 subject to w k , w l L 2 = δ kl (k, l = 1, . . . , q),
where δ kl is 1 when k = l and 0 otherwise. Therefore the estimator of

f i (i = 1, . . . , N ) is fi = f + q k=1 ψ ik w k , with ψ ik = f c i , w k L 2 (k = 1, . . . , q).
To solve this minimization problem, Ramsay and Silverman [START_REF] Ramsay | Functional Data Analysis[END_REF] proposes to express the sample functions on a spline basis. Then PCA can be applied to the coefficients of the functions on the spline basis. They also propose to apply PCA directly to the discretized functions. The FPCA decompostion ensures that I fi = I f i for i = 1, . . . , N so that the integrals of the sample probability density functions approximations fi are equal to one. However, the FPCA decomposition does not ensure that the estimators are non-negative.

Delicado [START_REF] Delicado | Dimensionality reduction when data are density functions[END_REF] proposes to apply FPCA to the logarithm transform g i = log f i to ensure the positivity of the prediction fi = exp ĝi (i = 1, . . . , N ). However this approach does not ensure that the approximations are normalized.

Affleck-Graves et al. [START_REF] Affleck-Graves | A principal component index subject to constraints[END_REF] have proposed to put a non-negativity constraint on the first basis function w 1 and let the other basis functions free. Indeed they claim that in general it is not possible to put such a constraint on the basis functions w 2 , . . . , w q . Hence the first basis function is forced to be non-negative but not the other ones. Therefore, it cannot be guaranteed that the approximations are non-negative.

Another method has been proposed by Kneip and Utikal [START_REF] Kneip | Inference for density families using functional principal component analysis[END_REF] to ensure that the FPCA approximations are non-negative. The FPCA method is applied to the sample functions. Then the negative values of the probability density functions approximations are interpreted as 0 and they are normalized to one. In the numerical studies of Sections 2 and 3, we will refer to this method as CPCA, for Constrained Principal Component Analysis.

Modified Magic Points (MMP)

The Empirical Interpolation Method (EIM) or Magic points (MP) [START_REF] Maday | A general, multipurpose interpolation procedure: the magic points[END_REF] is a greedy algorithm that builds interpolators for a set of functions. Here the set of functions under consideration is G (X , •) = {G(x, •) : x ∈ X }. The algorithm iteratively picks a set of basis functions in G (X , •) and a set of interpolation points in I to build an interpolator with respect to these points. Let x be in X and f x = G (x, •). At step q -1, the current basis functions and interpolation points are denoted w 1 , . . . , w q-1 : I → R and t i1 , . . . , t iq-1 ∈ I respectively. The interpolator I q-1 [f x ] of f x is defined as a linear combination of the basis functions:

I q-1 [f x ] = q-1 k=1 ψ k (x)w k ,
where the values of the coefficient functions ψ 1 , . . . , ψ q-1 : X → R at x are uniquely defined by the following interpolation equations (see [START_REF] Maday | A general, multipurpose interpolation procedure: the magic points[END_REF] for existence and uniqueness):

I q-1 [f x ] (t i k ) = f x (t i k ) (k = 1, . . . , q -1) .
At step q, one picks the parameter x iq such that f xi q = G x iq , • is the element of G (X , •) whose L 2distance from its own current interpolator I q-1 [f xi q ] is the highest. The next interpolation point t iq is the one that maximizes the gap between f xi q and I q-1 [f xi q ]. Then the new basis function w q is defined as a particular linear combination of f xi q and I q-1 [f xi q ]. Let us be more specific by summarizing the algorithm [START_REF] Bebendorf | Comparison of some reduced representation approximations[END_REF] hereafter.

Algorithm 1.1 (MP).

• Set q ← 1 and I 0 ← 0.

• While ε > tol do (a) choose x iq in X :

x iq ← arg sup x∈X f x -I q-1 [f x ] L 2
and the associated interpolation point t iq :

t iq ← arg sup t∈I f xi q (t) -I q-1 [f xi q ](t) , (b)
define the next element of the basis:

w q ← f xi q (•) -I q-1 [f xi q ](•) f xi q (t iq ) -I q-1 [f xi q ](t iq ) ,
(c) compute the interpolation error:

ε ← err q L ∞ (X ) , where err q : x → f x -I q [f x ] L 2 , (d) set q ← q + 1.
This method has been successfully applied e.g. to a heat conduction problem, crack detection and damage assessment of flawed materials, inverse scattering analysis [START_REF] Nguyen | Reduced-basis approximations and a posteriori error bound for nonaffine and nonlinear partial differential equations: Application to inverse analysis[END_REF], parameter-dependent convection-diffusion problems around rigid bodies [START_REF] Tonn | A reduced-basis method for solving parameter-dependent convection-diffusion problems around rigid bodies[END_REF], in biomechanics to describe a stenosed channel and a bypass configuration [START_REF] Rozza | Reduced basis methods for stokes equations in domains with non-affine parameter dependence[END_REF].

In general, the interpolators provided by the MP algorithm are not probability density functions. Therefore we modified the MP method so that the approximations are non-negative with integral equal to one. In the derived new method the interpolation is lost, but the basis functions are picked in a similar greedy way. In the following we denote A q-1 [f x ] the interpolator at step q -1 of the function f x associated with a parameter x in X :

A q-1 [f x ] = q-1 k=1 ψ k (x)w k ,
where the values of the coefficient functions at x are now defined as solutions of the following convex quadratic program: Minimize ψ1(x),...,ψq-1(x)∈R

f x - q-1 k=1 ψ k (x)w k 2 L 2 subject to        ψ k (x) ≥ 0 (k = 1, . . . , q -1) q-1 k=1 ψ k (x) = 1. (9)
At step q, the parameter x iq is chosen in the same way as in the MP algorithm. The new basis function w q is then f xi q . The Modified Magic Points (MMP) algorithm is detailed below.

Algorithm 1.2 (MMP).

• Set q ← 1 and I 0 ← 0.

• While ε > tol do (a) choose x iq in X :

x iq ← arg sup x∈X f x -A q-1 [f x ] L 2 ,
(b) define the next element of the basis:

w q ← f xi q ,
(c) compute the estimation error:

ε ← err q L ∞ (X ) , where err q : x → f x -A q [f x ] L 2 ,
(d) set q ← q + 1.

Remark 1.3. In step (a) of the MMP algorithm, the L 2 norm can be replaced by the Hellinger distance. This yields a new algorithm which is tested along with the previous one in Section 2.

In practice, we do not apply the MMP algorithm to the whole set G (X , •) since the computation of the functions is time-consuming. It is rather applied to the available sample set {f 1 , . . . , f N }. Thus the MMP algorithm is a way to select the most relevant functions among the sample f 1 , . . . , f N . Then the step (a) of the algorithm is only a finite maximization. Furthermore we build a regular grid of M points t 1 , . . . , t M in the interval I:

t j = t 1 + (j -1)∆t (j = 1, . . . , M ) with ∆t ∈ R + , (10)
and we discretize the functions f i (i = 1, . . . , N ) and w k (k = 1, . . . , q) on this grid:

f ij = f i (t j ) w kj = w k (t j ) (j = 1, . . . , M ).
Thus we can rewrite the problem (9) in the following way for all i ∈ {1, . . . , N }:

Minimize ψi1,...,ψiq M j=1 f ij - q k=1 ψ ik w kj 2 ∆t subject to      ψ i1 , . . . , ψ iq ≥ 0 q k=1 ψ ik = 1. (11)
These N problems are convex quadratic programs with q unknowns and q + 1 linear constraints, they are easily solvable for basis sizes smaller than 10 4 .

Remark 1.4. Let us suppose q = 1, i.e. the basis contains only one function. The equality constraint in the minimization [START_REF] Hardle | Optimal bandwidth selection in nonparametric regression function estimation[END_REF] reduces to ψ 1i = 1: the approximation is the same for each function f i in the sample.

Alternate Quadratic Minimizations (AQM)

Our third functional decomposition approach consists in tackling directly the problem consisting in minimizing the L 2 norm of the approximation error [START_REF] Delicado | Dimensionality reduction when data are density functions[END_REF] with basis functions w 1 , . . . , w q in F and coefficients ψ ik satisfying the constraints [START_REF] Da | Gaussian process modeling with inequality constraints[END_REF]. It differs from PCA because no orthonormality condition is imposed, and the decomposition is carried out directly on the raw data, i.e. without centering. This functional minimization program can be written as follows.

Minimize

w k ,ψ ik k=1,...,q,i=1,...,N

1 2 N i=1 f i - q k=1 ψ ik w k 2 L 2 subject to            w k ∈ F (k = 1, . . . , q) ψ ik ≥ 0 (i = 1, . . . , N )(k = 1, . . . , q) q k=1 ψ ik = 1 (i = 1, . . . , N ). ( 12
)
Just as before, instead of working on actual functions, we consider discretizations of them. We discretize the interval I in the same way [START_REF] Goldfarb | A numerically stable dual method for solving strictly convex quadratic programs[END_REF] and we set

f i = f i (t 1 ) • • • f i (t M ) (i = 1, . . . , N ) = f i1 • • • f iM w k = w k (t 1 ) • • • w k (t M ) (k = 1, . . . , q) = w k1 • • • w kM .
The functional minimization program ( 12) is then replaced by the following vectorial minimization program (in this section • R M designates the euclidean norm on R M ).

Minimize

w k ,ψ ik k=1,...,q,i=1,...,N

1 2 N i=1 f i - q k=1 ψ ik w k 2 R M subject to                        w kj ≥ 0 (k = 1, . . . , q)(j = 1, . . . , M ) M j=1 w kj ∆t = 1 (k = 1, . . . , q) ψ ik ≥ 0 (i = 1, . . . , N )(k = 1, . . . , q) q k=1 ψ ik = 1 (i = 1, . . . , N ). ( 13
)
Let us introduce additional notations to reformulate the objective function and the variables in a more compact way.

Ψ = [ψ ik ] i=1,...,N ;k=1,...,q =    ψ 1 . . . ψ N    W = [w kj ] k=1,...,q; j=1,...,M =    w 1 . . . w q   
We recall that the Frobenius norm of a matrix A = [a ij ] i=1,...,N ;j=1,...,M ∈ R N ×M is defined as

A F = N i=1 M j=1 a 2 ij .
This being said, the vectorial minimization program ( 13) can be written in matrix form.

Minimize

Ψ∈R N ×q + ,W∈R q×M + O (Ψ, W) = 1 2 F -ΨW 2 F
subject to

ψ i 1 = 1 (i = 1, . . . , N ) w k 1 = 1/∆t (k = 1, . . . , q). (14)
The objective function of the program ( 13) is second-order polynomial, hence twice continuously differentiable (but not necessarily convex). Its first-order derivatives are expressed as

∂ ψ ik O (Ψ, W) = -w k (f i -ψ i W) T (i = 1, . . . , N )(k = 1, . . . , q), ∂ w kj O (Ψ, W) = -Ψ T •,k (F •,j -ΨW •,j ) (k = 1, . . . , q) (j = 1, . . . , M ) ,
where Ψ •,k is the k th column of Ψ, F •,j and W •,j are the j th columns of F and W respectively, and here are the second-order derivatives which will be useful further down:

∂ 2 ψ ik ψ ik ′ O (Ψ, W) = w k w T k ′ (i = 1, . . . , N )(k, k ′ = 1, . . . , q), ∂ 2 w kj w kj ′ O (Ψ, W) = δ jj ′ Ψ •,k 2 R N (k = 1, . . . , q) (j, j ′ = 1, . . . , M ) .
The program [START_REF] Kneip | Inference for density families using functional principal component analysis[END_REF] has q(M + N ) ≈ 10 4 design variables, q(M + N ) + q + N ≈ 10 4 constraints and may not be convex. Hence its resolution through the use of a numerical solver may be very time-consuming (it turns out we could not manage to get a satisfactory feasible solution). In order to circumvent this issue we chose to implement [START_REF] Kneip | Inference for density families using functional principal component analysis[END_REF] as successive convex quadratic minimization programs. The value of the program can be written as follows.

inf

(Ψ,W)∈(R N ×M + ) 2 ψ i 1=1 (i=1,...,N ) w k 1=1/∆t (k=1,...,q) O (Ψ, W) = inf w1∈R M + w11=1/∆t • • • inf wq∈R M + wq1=1/∆t inf ψ 1 ∈R q + ψ 1 1=1 • • • inf ψ N ∈R q + ψ N 1=1 O (Ψ, W) .
The idea is to minimize the criterion O (Ψ, W) working on one line-vector at a time: ψ 1 , then ψ 2 , and so on until ψ N , and then w 1 , w 2 and so forth until w q . This process being repeated as many times as necessaryand affordable -to get some kind of convergence.

Let us make one more remark before giving out the algorithm. Let i ∈ {1, . . . , N }. The two following programs are equivalent -in the sense that their feasible and optimal solutions are the same (we just reduced the criterion to minimize).

Minimize

ψ i O (Ψ, W) subject to ψ i ∈ R q + ψ i 1 = 1 ⇐⇒ Minimize ψ i 1 2 f i -ψ i W 2 R M subject to ψ i ∈ R q + ψ i 1 = 1.
The algorithm we implemented to numerically solve ( 13) is the following.

Algorithm 1.5 (AQM).

• Initialize Ψ and W with uniform values.

ψ ik = 1/q (i = 1, . . . , N )(k = 1, . . . , q) w kj = 1/(M ∆t) (k = 1, . . . , q)(j = 1, . . . , M ).
• While some stopping criterion is not reached (i.e. a maximum number iter max ≈ 20 of iterations), (a) do, for i = 1, . . . , N , Minimize

ψ i 1 2 f i -ψ i W 2 R M subject to ψ i ∈ R q + ψ i 1 = 1. (b) do, for k = 1, . . . , q, Minimize w k O (Ψ, W) subject to      w k ∈ R M + w k 1 = 1/∆t.
Each of the minimization problems addressed in this algorithm has a convex quadratic criterion depending on at most M ≈ 2000 variables constrained by the same amount of lower bounds plus one linear equality. These optimization problems were solved with a dual method of the active set type, detailed in Goldfarb and Idnani [START_REF] Goldfarb | A numerically stable dual method for solving strictly convex quadratic programs[END_REF], specifically designed to address strictly convex quadratic programs. For a given quadratic program P , the active set of a point x is the set of (linear) inequality constraints satisfied with equality (active constraints). The algorithm starts from a point x such that the set A of its active constraints (possibly empty) is linearly independent, and the following steps are repeated until all constraints are satisfied. First, a violated constraint c is picked. Second, the feasibility of the subproblem P (A ∪ {c}) of P constrained only by the constraint c and those of the current active set A is tested. If P (A ∪ {c}) is infeasible then P is infeasible as well. Otherwise x is updated with a point whose active set is a linearly independent subset of A ∪ {c} including c. The loop terminates with an optimal solution.

Metamodelling of functional decomposition coefficients

In Sections 1.2.1, 1.2.2 and 1.2.3, three methods have been studied to approximate a sample of probability density functions on basis such that:

fi = q k=1 ψ ik w k .
The functions are characterized by their coefficients. The problem dimension is therefore reduced to the basis size q. A metamodel must now be designed to link the input parameters and the coefficients: a function ψ : X → R q , where ψ = ( ψ1 , . . . , ψq ). The metamodel is built thanks to the known points:

ψ ik = ψ k (x i ) (i = 1, . . . , N ,k = 1, . . . , q).
For a new x ∈ X , we look for the estimation fx of the output f x of G:

fx = q k=1 ψk (x)w k .
For MMP and AQM methods, the searched coefficients must respect the following constraints: for all x ∈ X ,

     q k=1 ψ k (x) = 1 ψ k (x) ≥ 0 (k = 1, . . . , q).
Therefore, the outputs of the metamodel to be built must be positive and have their sum equal to one. So far, the main approach investigated has been to build separate metamodels on each coefficient without constraining the output of the metamodel. When the metamodel is used to make predictions on new points in the input space, the predicted coefficients which are lower than zero are put equal to zero, and the coefficients are then renormalized such that their sum is equal to one. Any type of metamodel can be used in this approach. Polynomial regression, generalized additive models [START_REF] Trevor | Generalized Additive Models[END_REF] and Gaussian process metamodels have been tested, but the obtained results on different test cases are not convincing. This approach on the presented test cases does not lead to good results. The metamodel error is bigger than the error due to the decomposition on a functional basis.

Future research should be dedicated on the construction of more efficient metamodels taking into account these constraints. Metamodels which ensure the positivity of the output exist in the literature (see for example [START_REF] Da | Gaussian process modeling with inequality constraints[END_REF] for Gaussian process with inequality constraints). The main difficulty is to ensure that the sum of the coefficients is equal to one. A solution can be to consider the coefficients as realizations of a random variable on a simplex space. Aitchison [START_REF] Aitchison | The statistical analysis of compositional data[END_REF] proposes to apply a bijective transformation from this simplex space to R q . A metamodel without constraint could then be built between the input parameters and the transformed coefficients.

Numerical tests on toy functions

In this section, we apply the methods presented on two toy examples respectively for d = 1 in Section 2.1 and d = 5 in Section 2.2. We compare the relative error of the estimation of different quantities of interest: L 1 , L 2 and the Hellinger distances between two densities, mean, variance and the (1%, 25%, 75%, 99%)-quantiles. We recall that the Hellinger distance between f and g is the L 2 distance between √ f and √ g. The relative error is defined below for the three norms and the other scalar quantities of interest:

100 I f (t) -f (t) 2 dt I f (t) 2 dt 100 I f (t) -f (t) 2 dt I f (t) 2 dt 100 I f (t) -f (t) dt I |f (t)| dt 100 |u -û| |u| ,
where u is the mean, variance and the studied quantiles. This study is separated in two independent parts for the two toy examples. First, we compare the error obtained from the two estimators fL 2 and fHe based on kernel regression. We show the relative error in function of different sizes N 1 < • • • < N k of the design of experiments such that a design with N i points is included in the designs with N j points (i < j). For each N i , relative error for different quantities of interest are averaged on the same 1000 test points chosen uniformly in X . These computations are repeated 25 times with different designs of experiments.

In the second part, we compare the relative error obtained from the construction of basis obtained in Sections 1.2.1, 1.2.2 and 1.2.3 to reconstruct the N density probability functions of the design of experiments. These two methods decompose the known probability density functions on a functional basis and approximate them. However, no metamodel has been designed between the coefficients of the probability density functions on the basis and the parameters of the computer code, so that no estimation of an unknown probability density function can be done. Therefore, both methods cannot be compared to the kernel regression.

Each probability density function is discretized on 512 points.

A one-dimensional example (Toy example 1)

The first toy example is defined as follows:

G(x, ξ 1 , ξ 2 , U ) = (sin(x(ξ 1 + ξ 2 )) + U )1 {sin(x(ξ1+ξ2))+U ≥-1}
where Let x 0 ∈ [0, 1], we want to estimate the probability density function f 0 of the random variable G(x 0 , ξ 1 , ξ 2 , U ).

x ∈ X = [0, 1] ξ 1 ∼ N (1, 1) ξ 2 ∼ N (2, 1) U ∼ U ([0, 1]).

Kernel regression

In this section, one compares the estimators based on Hellinger distance and L 2 norm given respectively by equations ( 2) and [START_REF] Richard H Byrd | A limited memory algorithm for bound constrained optimization[END_REF]. Figure 4 represents the output for N = 100 points simulated uniformly and independently in X . Figure 5 represents for different values of x 0 the mean (in plain line) and the standard deviation (in dashed line) of the estimation obtained from the two estimators of the kernel regression in red for L 2 norms and in blue for the Hellinger distance. For these four parameters, the two estimators give approximately the same estimations. Figure 6 shows the boxplot of the relative error for L 2 norms and the Hellinger distance. In this first example the estimator based on the Hellinger distance seems slightly better than the one based on the L 2 norm. The errors are very close for the norms and the mean but are more different for other quantities of interest. The variance of the error decreases very quickly for norms. Moreover, the errors are low for most quantities of interest except for 1% and 25% quantiles.

Functional decomposition methods

In this section, the four functional decompositions presented in 1.2, MMP with L 2 and Hellinger distances, AQM and CPCA are compared. The decompositions are applied to learning samples whose sizes are N = 50 and 100 and for decomposition basis sizes ranging from 1 to 20. The relative errors between the learning sample functions and their approximations are computed for the 9 quantities of interest. Figure 7 represents, in logarithmic scale, these relative errors in function of the basis size for the four decompositions. First, the relative errors are low for all quantities of interest and especially for the norms, the mean and the higher quantiles. For the four decompositions, the relative error for norm decreases very evenly. The decrease is quite regular too for the mean, variance, 75% and 99% quantiles. The quality of approximation of the lower quantiles behaves more irregularly, as it sometimes increases with the basis size. CPCA method outperforms the three others for the three distances. This method seems better for all other quantities of interest except the 25% quantile, but the gap between the decompositions is smaller. The relative error for the learning sample of 100 densities is for most of the quantities higher than the relative error for the sample of 50 densities. For the same number of components q, it is indeed more difficult to approximate a higher number N of functions, so that the relative error increases with the size of the design of experiments for a constant number of components. Figure 8 represents the L 2 relative error of MMP (in red) with N = 200 and the relative error of twenty independent basis (in black) containing probability functions chosen randomly in the learning sample f 1 , . . . , f N . The MMP outperforms the random strategies. This result validates the way to choose a new curve among the sample of probability density functions, as this choice is better than randomly picking functions in the sample.

A five-dimensional example (Toy example 2)

Let x = (x 1 , x 2 , x 3 , x 4 , x 5 ) ∈ X = [0, 1] 5 .
A second toy example is defined as follows: Quantile 99% Figure 9 represents the output for N = 100 points simulated uniformly and independently in X .

G(x, N, U 1 , U 2 , B) = (x 1 + 2x 2 + U 1 ) sin(3x 3 -4x 4 + N ) + U 2 + 10x 5 B + 5 i=1 ix i where N ∼ N (0, 1) U 1 ∼ U ([0, 1]) U 2 ∼ U ([1, 2]) B ∼ Bern(1/2).
L2 • • • • • • • • • • • • • • • • • • • • • • • • • • 50
• • • • • • • • • • • • • • • • 50

Kernel regression

In this section, the kernel regression is tested on this second toy example with isotropic (Figure 11) and anisotropic bandwidth (Figure 12). As in the first example, Figure 10 represents for different values of x 0 the true probability density function (blue dashed line) and the estimation obtained from the two estimators of the kernel regression (red plain line and orange dots line). In this example, the two estimators give different results. In red is represented the relative error for the estimator based on Hellinger distance, in blue for the estimator based on L 2 norm. The Hellinger estimator gives a better approximation of the quantities presented in Figure 11, except for the mean. Contrary to the case of the first toy example, the errors are quite high. The variance of the error on the norms, presented in Figure 11, is low but does not seem to decrease steadily. The use of an anisotropic bandwidth does not seem to improve the quality of the estimation, while it is much longer to compute. The errors are approximately the same with isotropic or anisotropic bandwidth. The greater precision brought by the anisotropic bandwidth is compensated by the difficulty to estimate it.

Functional decomposition methods

In this section, the four functional decompositions MMP with L 2 and Hellinger distances, AQM and CPCA are compared. The decompositions are applied to learning samples whose sizes are N = 50 and 100 and for decomposition basis sizes ranging from 1 to 20. The relative errors between the learning sample functions and

• • • • • • • • • • • • • • • • • • • 5

L1

Basis size

Relative error

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 5

L2

Basis size

Relative error

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 5

Hellinger

Basis size

Relative error

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 5

Mean

Basis size

Relative error

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 5

Variance

Basis size

Relative error

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 5

Quantile 1%

Basis size

Relative error

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 5

Quantile 99%

Basis size

Relative error

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 5 10

Quantile 25%

Basis size

Relative error

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 5 10

Quantile 75%

Basis size

Relative error their approximations are computed for the 9 quantities of interest. Figure 13 represents, in logarithmic scale, these relative errors in function of the basis size for the three decompositions. First, the relative errors are low for all quantities of interest, as in the previous example. The relative errors decrease is less even than in the one-dimensional case. The behavior of the errors on 1% and 99% quantiles is particularly unsteady. CPCA method gives better results than the other methods for L 1 , L 2 norms, 25% and 75% quantiles. However, for all studied quantities, the AQM and CPCA method seem to give quite equivalent results. MMP method errors are slightly higher than AQM errors for most of the quantities.

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •

Industrial applications

In this section we apply the methods we developed in Section 1 to three industrial numerical codes of CEA, EDF, and IFPEN. We compute the same relatives errors as for the toy examples of Section 2.

CEA application: CASTEM test case

Code description

In the framework of nuclear plant risk assessment studies, the evaluation of component reliability during accidental conditions is a major issue required for the safety case. Thermal-hydraulic (T-H) and thermalmechanic (T-M) system codes model the behaviour of the considered component subjected to highly hypothetic accidental conditions. In the study that we consider here, the T-M code CASTEM takes as input 13 uncertain parameters, related to the initial plant conditions or to the safety system characteristics. Three of them are functional T-H parameters which depend on time: fluid temperature, flow rate and pressure. The other ten parameters are T-M scalar variables. For each set of parameters, CASTEM calculates the absolute mechanical strength of the component and the thermo-mechanical actual applied load. From these two elements, a safety margin (SM) is deduced. The objective is to assess how these uncertain parameters can affect the code forecasts and more specifically the predicted safety margin. However, CASTEM code is too time expensive to be directly used to conduct uncertainty propagation studies or global sensitivity analysis based on sampling methods. To avoid the problem of huge calculation time, it can be useful to replace CASTEM code by a metamodel. One way to fit a metamodel on CASTEM could be to discretize the functional inputs and to consider the values of the discretization as scalar inputs of CASTEM code. Nevertheless, this solution is often intractable due to the high number of points in the discretization. To cope with this problem, in [START_REF] Iooss | Global sensitivity analysis of computer models with functional inputs[END_REF][START_REF] Marrel | Global sensitivity analysis of stochastic computer models with joint metamodels[END_REF] a method was proposed to treat implicitly these "uncontrollable" parameters functional parameters, while the other ten scalar parameters are considered as "controllable". CASTEM output is then a random variable conditionally to "controllable" parameters.

A latin hypercube sampling method [START_REF] Michael D Mckay | Comparison of three methods for selecting values of input variables in the analysis of output from a computer code[END_REF] is used to build a learning sample of 500 points in dimension 10. For each set of controllable parameters, CASTEM has been run 400 times with different uncontrollable functional parameters, randomly chosen in an available database. The probability density function f i (i = 1, .., 500) of the safety margin is computed by kernel estimation with the 400 outputs of CASTEM for each set of parameters.

A few examples of the obtained probability density functions are represented on Figure 14. In the following,

• • • • • • • • • • • • • • • • • • • 5

L1

Basis size

Relative error

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 5 10 15
5e-02 5e-01 5e+00 5e+01

L2

Basis size

Relative error

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 5

Hellinger

Basis size

Relative error

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 5

Mean

Basis size

Relative error

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 5

Variance

Basis size

Relative error

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 5

Quantile 1%

Basis size

Relative error

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 5

Quantile 99%

Basis size

Relative error

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 5

Quantile 25%

Basis size

Relative error

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 5

Quantile 75%

Basis size

Relative error

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • Figure 13
. Toy example 2: comparison of the relative error for different quantities in function of the size of the basis q with a design of experiments of size N = 50 (circles) and 100 (filled circles). Blue: MMP decomposition given by L 2 norm. Red: MMP decomposition given by the Hellinger distance. Black: AQM decomposition. Green: CPCA decomposition.

the two kernel regression metamodels, then MMP, AQM and CPCA decomposition methods are applied on CASTEM test case.

Kernel regression

In this section, the kernel regression method is applied with isotropic and anisotropic bandwidths. Figure 15 represents for four different values of x 0 the true probability density function (blue dashed line) and the estimation obtained from the two estimators of the kernel regression (L 2 estimator in red plain line and Hellinger estimator in orange dotted line), with isotropic bandwidth. One can see that the estimation by the L 2 estimator of the four probability density functions is very far from the real function. In particular, at the bottom left of the figure, the mean of the predicted probability function in red is very far from the mean of the real one in dashed blue. On these four probability density functions, the Hellinger estimator gives much better results than the L 2 one.

To verify this first graphical analysis, the Leave-One-Out method has been used to assess the efficiency of the kernel regression-based metamodel. The bandwidth of the kernel regression is estimated thanks to all probability density functions except the i th one. The function f i is estimated with the corresponding metamodel and the relative errors on the quantities of interest are computed between f i and its estimation. This process is repeated for each probability density function in the dataset and all the computed relative errors are averaged. The mean relative errors of the estimators based on Hellinger distance and L 2 norm are given in 1. CASTEM test case: the mean relative errors on the quantities of interest computed by the Leave-One-Out method for the kernel regression estimators based on the Hellinger distance and L 2 norm with an isotropic bandwidth.

Table 2 for the different quantities of interest respectively for an isotropic and anisotropic bandwidths. The two estimators perform poorly on this dataset. The errors for L 1 and L 2 norms are particularly high. The Leave-One-Out validation confirms that the kernel regression method is not adapted to CASTEM test case. This can be explained by the important influence of the controllable T-M parameters compared to the one of the uncontrollable T-H parameters.

Functional decomposition methods

In this section, the four functional decompositions MMP with L 2 and Hellinger distances, AQM and CPCA are compared. Figure 16 represents, in logarithmic scale, the relative errors on the different norms, modes and quantiles versus the basis size (from 1 to 20). First, the relative errors are low for all quantities of interest, except for the variance. For the variance, the errors are over 10% with a decomposition basis with 20 functions. The errors for small bases are high but decrease very quickly. The decrease is especially quick and even for errors on L 1 , L 2 and Hellinger distances. For modes and quantiles, the errors of the four methods do not decrease steadily for all quantities of interest, and especially for small basis sizes. CPCA clearly outperforms other methods for the L 1 , L 2 and Hellinger distances. For other quantities of interest, it gives good results. For 2. CASTEM test case: the mean relative errors on the quantities of interest computed by the Leave-One-Out method for the kernel regression estimators based on the Hellinger distance and L 2 norm with an anisotropic bandwidth. the variance, 1% and 99% quantiles, AQM method performs poorly compared to the three others. The errors of AQM and CPCA methods seem more stable than for MMP method for most of the quantities of interest, especially for modes and quantiles. Overall, the results are quite good for the four methods.

• • • • • • • • • • • • • • • • • • • 5

L1

Basis size

Relative error

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 5

L2

Basis size

Relative error

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 5

Hellinger

Basis size

Relative error

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 5

Mean

Basis size

Relative error

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 5

Variance

Basis size

Relative error

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 5

Quantile 1%

Basis size

Relative error

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 5

Quantile 99%

Basis size

Relative error

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 5 10 15
5e-02 5e-01 5e+00 5e+01

Quantile 25%

Basis size

Relative error

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 5

Quantile 75%

Basis size

Relative error

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • Figure 16
. CASTEM test case: comparison of the relative error for the 9 quantities in function of the size of the basis q. Blue: MMP decomposition given by L 2 norm. Red: MMP decomposition given by the Hellinger distance. Black: AQM decomposition. Green: CPCA decomposition.

EDF application: VME test case

Code description

To optimize the whole life cost of its nuclear fleet, EDF has developed an asset management methodology [START_REF] Fessart | Optimisation d'un programme d'investissements : méthode IPOP. Conférence lambdamu[END_REF]. A part of this methodology deals with exceptional maintenance tasks strategies. To help the decision maker to choose the best strategy (how many times do we need to carry out exceptional tasks, when,...?), EDF has developed a dedicated tool called VME (described in Figure 17) based on Monte-Carlo simulation to compute many technical economic indicators among which the density function of the Net Present Value (called VAN for "Valeur Actuelle Nette") is the most relevant.

This tool leads to an important simulation time and requires an important amount of input data that are surrounded with uncertainties:

• Reliability data: generally there is not enough (or sometimes not any) feedback data to precisely evaluate reliability model parameters; • Economic data: economic indicators and duration of maintenance tasks remain on several hypothesis that can be modified; • Other data: uncertainty on operating times of power plants, maintenance tasks dates, etc. In our work, we study a particular scenario simulation with VME concerning the life cycle management of four large transformers, adapted from the EPRI study made for CENG [START_REF] Lonchampt | An optimization approach for life cycle management applied to large power transformers[END_REF]. The uncertainty on the inputs are 3. VME: uncertainty inputs used in the model. Each distribution depends of three parameters: the shape, the scale and the location. represented by four Weibull distributions which are summarized in Table 3. Define for i = 1, . . . , 4,

D i = {x = (x 1 , x 2 , x 3 ) ∈ R 3 : E[Weibull(x 1 , x 2 , x 3 )] ∈ [E[W i ] ± 0.2E[W i ]]},
and denote sh i , sc i , loc i respectively the shape, the scale and the location of W i . For i = 1, . . . , 4 let x i,1 , x i,2 , x i,3 , x i,4 be parameters simulated uniformly on D i . Note that W 2 and W 4 follow the same law but sh 2 , sc 2 , loc 2 and sh 4 , sc 4 , loc 4 are considered as different parameters. That means there are three parameters for each Weibull distribution. Moreover, there are four other parameters which represent some costs and two other parameters which represent durations of services. All these 18 parameters are summarized in Table 4.

In this case, we are particulary interested by the VAN indicator. It is defined by

V AN = P(G(x, W) ≤ 0),
where W = (W 1 , . . . , W 4 ) and x the eighteen fixed parameters. Figure 18 represents 75 outputs from the VME model. The kernel regression is applied to designs of experiments of increasing sizes such that smaller designs are included in greater ones. A test sample of 500 points is used to assess the efficiency of the regression. In this section, the kernel regression is tested on this real case study with isotropic and anisotropic bandwidth (Figure 20). In red is represented the relative error for the estimator based on the Hellinger distance, in blue for the estimator based on the L 2 norm. The estimation seems equivalent for the two estimators. The results are equivalent for anisotropic and isotropic bandwidth. For the anisotropic bandwidth, the estimation is less stable and the estimation of H becomes much more time-consuming.

Kernel regression

The relative error in norm is smaller for VME (with isotropic and anisotropic bandwidth) than for the toy example 2, whereas the dimension is higher. That shows that the quality of estimation depends not only on the dimension but also on the regularity of the model, if the bandwith is sufficiently well estimated. Bad results are found for the variance, 99% and 99.9% quantiles. The relative error for the mean is also very high. The real mean is quite close to zero, so that in the relative error computation, the numerator is divided by a quantity close to zero.

Figure 21 represents the estimation of the P(V AN < 0) for isotropic and anisotropic bandwidths. In this case, choosing an anisotropic bandwidth does not bring information.

Functional decomposition methods

In this section, the four functional decompositions MMP with L 2 and Hellinger distances, AQM and CPCA are compared. Figure 22 represents, in logarithmic scale, the relative errors on the different norms, modes and quantiles versus the basis size (from 1 to 20). First, the relative errors are low for all quantities of interest, except for the mean and variance. For the variance, the errors are over 20% with a decomposition basis with 20 functions. The very high errors on the mean can be explained by the fact that the mean to be estimated is close to zero for most of densities in the VME case. The errors on mean, variance have moreover a very chaotic behavior. The errors on distances decrease very quickly. For modes and quantiles, the errors of the four methods do not decrease steadily for all quantities of interest. CPCA outperforms other methods for the L 1 , L 2 , Hellinger distances, the variance and the 99% and 25% quantiles. The performances of AQM and MMP are very close for modes and quantiles. CPCA seems more stable than MMP and AQM on most of the quantities of interest.

In Figure 23, the errors on the probability for the VAN to be negative P (V AN < 0) for the different decomposition methods are represented in function of the basis size. The errors are quite low and decrease quickly. The CPCA method outperforms the others for both learning sets sizes. The other three methods give similar results.

• • • • • • • • • • • • • • • • • • • 5

L1

Basis size

Relative error

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 5 10 15
5e-03 5e-02 5e-01 5e+00

L2

Basis size

Relative error

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 5

Hellinger

Basis size

Relative error

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 5 10

Mean

Basis size

Relative error

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 5

Variance

Basis size

Relative error

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 5

Quantile 1%

Basis size

Relative error

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 5

Quantile 99%

Basis size

Relative error

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •

Quantile 25%

Basis size

Relative error

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 5

Quantile 75%

Basis size

Relative error

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • Figure 22
. VME: comparison of the relative error for the 9 quantities in function of the size of the basis q. Blue: MMP decomposition given by L 2 norm with a learning set size of N = 50 (circles) and 100 (filled circles). Red: MMP decomposition given by the Hellinger distance. Black: AQM decomposition. Green: CPCA decomposition. Relative error

• • • • • • • • • • • • • • • • • • • 5
• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • Figure 23
. VME: comparison of the relative error for P (V AN < 0) in function of the size of the basis q. Blue: MMP decomposition given by L 2 norm with a learning set size of N = 50 (circles) and 100 (filled circles). Red: MMP decomposition given by the Hellinger distance. Black: AQM decomposition. Green: CPCA decomposition. At IFP Energies Nouvelles, the researchers of the Department of Thermodynamics and Molecular Modelling elaborate models for the structure of molecules such as hydrocarbons and alcohols. In the following numerical study we consider a system of pentane molecules at a reduced temperature of 0.75. Pentane is an alkane with five carbon atoms (C 5 H 12 , cf. Figure 24.) Here we use the Anisotropic United Atoms model (AUA) [START_REF] Ungerer | Optimization of the anisotropic united atoms intermolecular potential for n-alkanes[END_REF] to represent the molecular structure of pentane. The carbon (C) and hydrogen (H) atoms in each of the two terminal methyls (CH 3 ) and each of the three methylene bridges (CH 2 ) are treated as a single interaction center. To each of these two types of interaction centers correspond three parameters of the model: an energy parameter (ε CH3 and ε CH2 ), a size parameter (σ CH3 and σ CH2 ), and a displacement parameter (δ CH3 and δ CH2 ). These six parameters will be the components of the set of input parameters of the numerical codes under consideration:

x = (ε CH3 , ε CH2 , σ CH3 , σ CH2 , δ CH3 , δ CH2 ) ∈ R 6 .
Here the matter of interest is the prediction of two macroscopic properties of a pentane system in chemical equilibrium: the volumetric mass density ρ eq liq of the liquid phase and the vapor pressure P eq gas (pressure of the gas phase). The algorithm carried out to numerically predict ρ eq liq and P eq gas -for a given set of parameters x -relies on the fundamental postulate of statistical mechanics: an isolated system in equilibrium is found with equal probability in each of its accessible microstates. Thus, ρ eq liq (respectively P eq gas ) can be approximated by the average of the values of the volumetric mass density of the liquid phase ρ liq (respectively the vapor pressure P gas ) of a large number of accessible microstates of the system. Therefore the algorithm is essentially a loop, each step consisting of the following two instructions: (1) randomly generate an accessible microstate of the system (characterized by the positions of the molecules and the volume of the system), (2) compute the value of ρ liq (respectively P gas ) given by the model for this particular microstate. Let us denote G ρ liq (x) (respectively G Pgas (x)) the probability density function of the random value ρ liq (respectively P gas ) computed at each loop step. The two instructions above are iterated a hundred million times. After the execution of the loop, the arithmetical mean of these millions of values is computed -it is a Monte Carlo approximation of the integral of G ρ liq (x) (respectively G Pgas (x)) -and is returned by the algorithm as the prediction of ρ liq (respectively P gas .) The whole algorithm is carried out in about twenty-four hours on a supercomputer (for each set x of input parameters.)

In order to circumvent this time-consuming process, we want to build a metamodel of the numerical codes mapping the input parameters x to the probability density functions G Pgas (x) and G ρ liq (x) of the random outputs. We had at our disposal a sample of size N = 50. Figure 25 and Figure 26 show the corresponding 50 probability functions for volumetric mass density and vapor pressure respectively.

Kernel regression

We applied the kernel regression method on the numerical codes providing the volumetric mass density of the liquid phase and the vapor pressure. The results provided by the L 2 estimator with isotropic bandwidth are provided in Table 5 for the volumetric mass density and in been computed by the Leave-One-Out methodology as in Section 3.1.2. We can see that most of the results are disappointing: all relative errors are higher than 30%, excepted -surprisingly -the relative Hellinger error of the approximation. The latter amounts to nearly 7% for the volumetric mass density, and only nearly 0.05% for the vapor pressure. The errors for quantile 1% and quantile 25% are particularly high. The Hellinger estimator give almost the same results: the difference between the error for the L 2 estimator and the error for the Hellinger estimator is lower than 10 -5 in every case. That is why they are not charted here.

In Figure 27 and Figure 28 we plot the probability density functions of the volumetric mass density and the vapor pressure with their kernel regression approximations, for arbitrarily chosen sets of input parameters. We can see how poorly the densities are approximated. 6. GIBBS test case: KR relative errors of the L 2 estimator for the vapor pressure.

Functional decomposition methods

In this section, the four functional decompositions MMP with L 2 and Hellinger distances, AQM and CPCA are compared. Figure 29 represents, in logarithmic scale, the relative errors on the different norms, modes and quantiles versus the basis size (from 1 to 20) for the volumetric mass density of the liquid phase. The error curves are globally decreasing excepted the one corresponding to the variance relative error for AQM. Moreover, the relative error on variance is high even for basis with 20 functions and especially for AQM. The behavior of the AQM approximation is more erratic for the 99% quantile. Both MMP methods perform worse on the 25% and 75% quantiles than other methods. Other relative errors are low. CPCA gives better results for L 1 , L 2 , Hellinger distances and the 25% quantile. For 25% and 75% quantiles, errors on both MMP bases are higher than for others methods.

The results for the vapor pressure are presented in Figure 30. The errors on all quantities of interest are more steady than for the volumetric mass. In this case, CPCA and AQM have very different behavior compared to both MMP methods. For all quantities of interest except L 1 , CPCA and AQM errors have a plateau for higher basis sizes. For higher basis sizes, MMP method performs better than AQM and CPCA for all studied quantities, except the 25% quantile. Thus, contrary to what happens in the other presented numerical examples, CPCA method does not perform better than AQM and MMP for distances. For MMP decompositions, the errors on the mean and quantiles display a jump for basis sizes of 12 or 13.

• • • • • • • • • • • • • • • • • • • 5

L1

Basis size

Relative error

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 5

L2

Basis size Relative error

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 5

Hellinger

Basis size

Relative error

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 5

Mean

Basis size

Relative error Relative error

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 5
• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 5

Quantile 1%

Basis size

Relative error

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 5

Quantile 99%

Basis size

Relative error

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 5 10 15
5e-01 5e+00 5e+01 5e+02

Quantile 25%

Basis size

Relative error

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •

Quantile 75%

Basis size

Relative error

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • Figure 29
. GIBBS test case: comparison of the relative errors on different quantities for the volumetric mass density in function of the size of the basis q. Blue: MMP decomposition given by L 2 norm. Red: MMP decomposition given by the Hellinger distance. Black: AQM decomposition. Green: CPCA decomposition.

To this mean, we first adapted the well-known kernel regression developed by [START_REF] Elizbar | On estimating regression[END_REF] and [START_REF] Watson | Smooth regression analysis[END_REF]. Second, we proposed to approximate the sample probability density functions on a functional basis in order to reduce the problem dimension. In this way, the probability density functions are characterized by their coefficients on the basis, so that the problem becomes finite-dimensional. Two methods have been proposed to build a functional basis. The first method is adapted from Magic Points Method [START_REF] Maday | A general, multipurpose interpolation procedure: the magic points[END_REF] and builds iteratively the basis by adding the sample function that maximizes the approximation error on the previous basis. The second one, called Alternate Quadratic Minimization (AQM), aims at minimizing the L 2 approximation error under the defined constraints. Both methods require the resolution of constrained optimization problems. Then, in future work, a metamodel could be adjusted on the coefficients, to link them with the uncertain inputs, and provide a global metamodel for probability density function outputs.

Our methods have been tested on two analytical test cases and three industrial applications proposed by CEA, EDF and IFPEN. The kernel regression-based metamodel performs well on the analytical test cases, but shows its limits on the three numerical codes. Indeed, kernel regression is efficient in low dimension, but the difficulty to estimate the parameters increases with the dimension and the size of the sample. Moreover, this method is not adapted for numerical codes for which the influence of the controllable parameters is much more important than one of the uncontrollable parameters. Both Modified Magic Points (MMP) and AQM methods give good results on all the test cases. However, in most cases, the method proposed by Kneip and Utikal [START_REF] Kneip | Inference for density families using functional principal component analysis[END_REF] performs better than the proposed ones.

L1

Basis size

Relative error

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 5 10 15
5e-01 5e+00 5e+01 5e+02

L2

Basis size

Relative error

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 5

Hellinger

Basis size

Relative error

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •

Mean

Basis size

Relative error

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 5

Variance

Basis size

Relative error

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •

Quantile 1%

Basis size

Relative error

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •

Quantile 99%

Basis size

Relative error

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •

Quantile 25%

Basis size

Relative error

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •

Quantile 75%

Basis size

Relative error

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • Figure 30
. GIBBS test case: comparison of the relative errors on different quantities for the vapor pressure in function of the size of the basis q. Blue: MMP decomposition given by L 2 norm. Red: MMP decomposition given by the Hellinger distance. Black: AQM decomposition. Green: CPCA decomposition.

Metamodels linking the uncertain inputs of the code and the functional basis coefficients remains to be defined. The main difficulty to build this metamodel lies in ensuring that the sum of the functional basis coefficients is equal to one. One of the limits of our methods is that all quantities of interest are not as well approximated. In particular, quantiles are often poorly approximated. Indeed, all the proposed methods are based on norm minimization, and the norms used are not well suited for the study of extreme values. It could be interesting to adapt the proposed methods to norms giving more weights on the tails of the densities. The use of the Wasserstein metric between probability measures will be also studied.

A. Appendix

A.1. L 2 distance-based kernel estimator

We prove here the equivalence between the classical kernel estimator and the L 2 distance-based kernel estimator introduced in Sections 1.1.1 and 1.1.2 respectively. Indeed the classical kernel estimator (1) can be retrieved by writing the optimality condition of the minimization problem (3) defining the L 2 distance-based kernel estimator:

N i=1 -2K H (x i , x 0 ) I f i -f0 = 0 ⇐⇒ I N i=1 K H (x i , x 0 ) f i -f0 = 0 ⇐⇒ N i=1 K H (x i , x 0 ) f i -f0 = 0 ⇐⇒ f0 = N i=1 K H (x i , x 0 ) N j=1
K H (x j , x 0 ) f i .

A.2. Hellinger distance-based kernel estimator

We prove here the equivalence between expressions (4) and ( 5) of the Hellinger distance-based kernel estimator given in Section 1.1.2. The constraint that the integral of f0 is 1 in the optimization problem (4) is handled using the associated Lagrangian function. The problem becomes f0 = arg min

f,λ L(f, λ), with L(f, λ) = N i=1 K H (x i , x 0 ) I f i -f 2 -λ   I f -1   .
The first order optimality conditions are

               0 = ∂L ∂f = N i=1 K H (x i , x 0 ) I   1 - √ f i f0   -λ I 1, 0 = ∂L ∂λ = 1 - I f0 . (15)

10 Figure 1 .

 101 Figure 1. Probability density function given by the numerical code G for different values of the input parameters.

2 L 2 Figure 2 .

 222 Figure 2. Left: the spots correspond to design points in X . Right: the associated probability density functions given by the code G defined in the introduction. The probability density function associated to the black square (dashed line) in left Figure seems to be closer to the probability density function represented in green than to the others.

Figure 3 .

 3 Figure3. The orange curve is the Hellinger kernel prediction for the (unknown) output in black, built from the dark blue, red, light blue and green learning probability density functions (N = 4).

Figure 4 .

 4 Figure 4. Toy example 1: representations of outputs for N = 100 different parameters.

Figure 5 .

 5 Figure 5. Toy example 1: curves in red (resp. blue) plain line represents the estimation of a density for a fixed x 0 from the first (resp. second) method and dashed red (resp. blue) line represent the standard deviation of the estimation obtained from 25 independent design of experiments. The blue and the red lines are superposed.

Figure 6 .

 6 Figure 6. Toy example 1: boxplot of the errors for different sizes of N . Blue: estimator given by L 2 norm. Red: estimator given by the Hellinger distance. Results have been averaged with 25 independent experiments. Blue: estimator given by L 2 norm. Red: estimator given by the Hellinger distance.

Figure 7 .

 7 Figure 7. Toy example 1: comparison of the relative error for different quantities in function of the size of the basis q with a design of experiments of size N = 50 (circles) and 100 (filled circles). Blue: MMP decomposition given by L 2 norm. Red: MMP decomposition given by the Hellinger distance. Black: AQM decomposition. Green: CPCA decomposition.

Figure 8 .

 8 Figure 8. Toy example 1: comparison of the L 2 relative error obtained with MMP (in red) and a random choice in the design of experiments (in black).

Figure 9 .Figure 10 .Figure 11 .

 91011 Figure 9. Toy example 2: representations of outputs for N = 100 different parameters.

Figure 12 .

 12 Figure 12. Toy example 2: boxplot of the errors for different sizes of N . Left: estimator given by L 2 norm. Right: estimator given by the Hellinger distance. Results have been averaged with 25 independent experiments. These estimations were made with h anisotropic. Blue: estimator given by L 2 norm. Red: estimator given by the Hellinger distance.

Figure 15 .

 15 Figure 15. CASTEM test case: curves in dashed gray lines are probability density functions, red plain line represents the estimation of a density for a fixed x 0 and dashed blue line is the true probability density function.

Figure 17 . 1

 171 Figure 17. VME: EDF valorization of maintenance strategies tool description.

Figure 19

 19 Figure 19 represents for different values of x 0 the true probability density function (blue dashed line) and the estimation obtained from the two estimators of the kernel regression (red plain line and orange dots line).

Figure 18 .

 18 Figure 18. VME: representations of outputs for N = 75 different parameters.

Figure 19 .

 19 Figure 19. VME: curves in dashed gray lines are probability density functions, red plain line represents the estimation of a density for a fixed x 0 and dashed blue line is the true probability density function for N = 100.

Figure 24 .

 24 Figure 24. The structural formula of pentane. (Wikipedia)

Table 1 and

 1 Figure 14. CASTEM test case: representations of outputs for N = 75 different parameters.

	1.4
	1.2
	1.0
	0.8
	0.6
	0.4
	0.2
	0.0
	0	100	200	300	400	500
			Hellinger distance estimator L 2 norm estimator
	L 1 norm	114.2%		114.9%
	L 2 norm	167.7%		171.9%
	Hellinger distance	17.0%		17.0%
	Mean		22.2%		22.4%
	Variance	86.1%		86.9%
	1% quantile	72.0%		74.3%
	99% quantile	39.7%		40.9%
	25% quantile	30.9%		31.5%
	75% quantile	23.7%		23.9%
	Table

Table 4 .

 4 VME: distribution of the eighteen parameters.

	Uncertainty inputs distribution	support
	(sh 1 , sc 1 , loc 1 ) (sh 2 , sc 2 , loc 2 ) (sh 3 , sc 3 , loc 3 ) (sh 4 , sc 4 , loc 4 ) Cost1 Cost2 Cost3 Cost4 Duration1 Duration2	U nif orm U nif orm U nif orm U nif orm U nif orm U nif orm U nif orm U nif orm U nif orm U nif orm	D 1 D 2 D 3 D 4 [796 , 1194] [49.6 , 74.4] [5.6 , 8.4 ] [4 , 8 ] [0.8 , 1.2 ] [0.8 , 1.2]

  Table 6 for the vapor pressure. The errors have

		0.25				1 2	14 15	27 28	40 41
						3	16	29	42
		0.20				4 5 6	17 18 19	30 31 32	43 44 45
						7	20	33	46
	t(densite)	0.10 0.15				13 8 9 10 12 11	26 21 22 23 25 24	39 34 35 36 38 37	47 48 49 50
		0.05
		0.00
		0	200	400	600	800	1000
				t
		Figure 25. GIBBS test case: distributions of the volumetric mass density for different sets of
		input parameters.
		4e-05				1 2	14 15	27 28	40 41
						3	16	29	42
						4	17	30	43
		3e-05				5 6 7	18 19 20	31 32 33	44 45 46
						8	21	34	47
	t(densite)	2e-05				9 10 12 11	22 23 25 24	35 36 38 37	48 49 50
						13	26	39
		1e-05
		0e+00
		0e+00	1e+06	2e+06		3e+06	4e+06
				t
		Figure 26. GIBBS test case: distributions of the vapor pressure for different sets of input
		parameters.

Table 5 .

 5 GIBBS test case: kernel regression relative errors of the L 2 estimator for the volumetric mass density.

	L1	L2	Hellinger
	151.1%	375.9%	7.089%
	Mean	Variance	Quantile 1%
	50.77%	98.07%	513.8%
	Quantile 99% Quantile 25% Quantile 75%
	30.44%	198.4%	41.52%
	L1	L2	Hellinger
	123.7%	254.1%	0.04878%
	Mean	Variance	Quantile 1%
	69.86%	80.72%	667.63%
	Quantile 99% Quantile 25% Quantile 75%
	74.82%	156.1%	73.5%
	Table

We would like to thank our supervisors Nicolas Bousquet (EDF), Frédéric Delbos (IFPEN), Bertrand Iooss (EDF), Jérôme Lonchampt (EDF), Rafael Lugo (IFPEN), Amandine Marrel (CEA) and Nadia Pérot (CEA). We would also like to thank Anthony Nouy ( École Centrale de Nantes) and Sébastien Da Veiga (IFPEN) for their guidance. Furthermore, we would like to thank Nicolas Champagnat (Institut Élie Cartan de Lorraine) and Tony Lelièvre ( École des Ponts ParisTech) who organized the CEMRACS 2013 in collaboration with Anthony Nouy. Finally, we would like to thank CEA, EDF and IFPEN for funding this research project.

Hellinger

Sample size

Relative error

Quantile 99%

Sample size

Relative error Quantile 25%

Sample size

Relative error

Conclusion

The aim of this work was to build a metamodel for code outputs which are probability density functions. A first idea to design such a model was to build a convex combination of the sample probability density functions.

The first condition in [START_REF] Lonchampt | An optimization approach for life cycle management applied to large power transformers[END_REF] gives

Then we replace f0 in the second condition of ( 15) by the expression given in ( 16):

The combined equations ( 16) and ( 17) entail the formula given by ( 5) for the kernel estimator associated to the Hellinger distance: