Keywords: Gaussian process, Metamodel, Quantile, Screening, Sobol' indices

In the framework of the estimation of safety margins in nuclear accident analysis, a quantitative assessment of the uncertainties tainting the results of computer simulations is essential.

Accurate uncertainty propagation (estimation of high probabilities or quantiles) and quantitative sensitivity analysis may call for several thousand of code simulations. Complex computer codes, as the ones used in thermal-hydraulic accident scenario simulations, are often too cputime expensive to be directly used to perform these studies. A solution consists in replacing the computer model by a cpu inexpensive mathematical function, called a metamodel, built from a reduced number of code simulations. However, in case of high dimensional experiments (with typically several tens of inputs), the metamodel building process remains difficult. To face this limitation, we propose a methodology which combines several advanced statistical tools: initial space-filling design, screening to identify the non-influential inputs, Gaussian process (Gp) metamodel building with the group of influential inputs as explanatory variables. The residual effect of the group of non-influential inputs is captured by another Gp metamodel. Then, the resulting joint Gp metamodel is used to accurately estimate Sobol' sensitivity indices and high quantiles (here 95%-quantile). The efficiency of the methodology to deal with a large number of inputs and reduce the calculation budget is illustrated on a thermal-hydraulic calculation case simulating with the CATHARE2 code a Loss Of Coolant Accident scenario in a Pressurized Water Reactor. A predictive Gp metamodel is built with only a few hundred of code simulations and allows the calculation of the Sobol' sensitivity indices. This Gp also provides a more accurate estimation of the 95%-quantile and associated confidence interval than the empirical approach, at equal calculation budget. Moreover, on this test case, the joint Gp approach outperforms the simple Gp.

I. INTRODUCTION

Best-estimate computer codes are increasingly used to estimate safety margins in nuclear accident management analysis instead of conservative procedures [START_REF] Pourgol-Mohamad | Integrated methodology for thermal-hydraulic code uncertainty analysis with application[END_REF][START_REF] Bucalossi | Comparison between bestestimate-plus-uncertainty methods and conservative tools for nuclear power plant licensing[END_REF]. In this context, it is essential to evaluate the accuracy of the numerical model results, whose uncertainties come mainly from the lack of knowledge of the underlying physics and the model input parameters. The socalled Best Estimate Plus Uncertainty (BEPU) methods were then developed and introduced in safety analyses, especially for thermal-hydraulic issues, and even more precisely for the large-break loss-of-coolant accident [START_REF] Prosek | The state-of-the-art theory and applications of best-estimate plus uncertainty methods[END_REF][START_REF] Wilson | Historical insights in the development of Best estimate Plus Uncertainty safety analysis[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]. Its main principles rely on a probabilistic modeling of the model input uncertainties, on Monte Carlo sampling for running the thermal-hydraulic computer code on sets of inputs, and on the application of statistical tools (based for example on order statistics as the Wilks' formula) to infer high quantiles of the output variables of interest [START_REF] Nutt | Evaluation of nuclear safety from the outputs of computer codes in the presence of uncertainties[END_REF][START_REF] Wallis | Uncertainties and probabilities in nuclear reactor regulation[END_REF][START_REF] De Crécy | Uncertainty and sensitivity analysis of the LOFT L2-5 test: Results of the BEMUSE programme[END_REF][START_REF] Petruzzi | Approaches, relevant topics, and internal method for uncertainty evaluation in predictions of thermal-hydraulic system codes[END_REF][START_REF] Martin | Perspectives on the application of order-statistics in best-estimate plus uncertainty nuclear safety analysis[END_REF].

Not restricted to the nuclear engineering domain, the BEPU approach is more largely known as the uncertainty analysis framework [START_REF]Uncertainty in industrial practice[END_REF]. Quantitative assessment of the uncertainties tainting the results of computer simulations is indeed a major topic of interest in both industrial and scientific communities. One of the key issues in such studies is to get information about the output when the numerical simulations are expensive to run. For example, one often faces up with cpu time consuming numerical models and, in such cases, uncertainty propagation, sensitivity analysis, optimization processing and system robustness analysis become difficult tasks using such models. In order to circumvent this problem, a widely accepted method consists in replacing cputime expensive computer models by cpu inexpensive mathematical functions (called metamodels) based, e.g., on polynomials, neural networks, or Gaussian processes [START_REF] Fang | Design and modeling for computer experiments[END_REF]. This metamodel is built from a set of computer code simulations and must be as representative as possible of the code in the variation domain of its uncertain parameters, with good prediction capabilities. The use of metamodels has been extensively applied in engineering issues as it provides a multi-objective tool [START_REF] Forrester | Engineering design via surrogate modelling: a practical guide[END_REF]: once estimated, the metamodel can be used to perform global sensitivity analysis, as well as uncertainty propagation, optimization, or calibration studies. In BEPU-kind analyses, several works [START_REF] Cannamela | Controlled stratification for quantile estimation[END_REF][START_REF] Zio | Neural network and Order Statistics for Quantifying nuclear power plant Safety Margins[END_REF][START_REF] Lorenzo | Assessment of an isolation condenser of an integral reactor in view of uncertainties in engineering parameters[END_REF] have introduced the use of metamodels and shown how this technique can help estimate quantiles or probability of failure in thermal-hydraulic calculations.

However, the metamodeling technique is known to be relevant when simulated phenomena are related to a small number of random input variables (see [START_REF]Springer Handbook on Uncertainty Quantification[END_REF] for example). In case of high dimensional numerical experiments (with typically several tens of inputs), and depending on the complexity of the underlying numerical model, the metamodel building process remains difficult, or even impracticable. For example, the Gaussian process (Gp) model [START_REF] Santner | The design and analysis of computer experiments[END_REF] which has shown strong capabilities to solve practical problems, has some caveats when dealing with high dimensional problems. The main difficulty relies on the estimation of Gp hyperparameters. Manipulating predefined or well-adapted Gp kernels (as in [START_REF] Muehlenstaedt | Data-driven Kriging models based on FANOVA-decomposition[END_REF][START_REF] Durrande | ANOVA kernels and RKHS of zero mean functions for model-based sensitivity analysis[END_REF]) is a current research way, while several authors propose to couple the estimation procedure with variable selection techniques [START_REF] Welch | Screening, predicting, and computer experiments[END_REF][START_REF] Marrel | An efficient methodology for modeling complex computer codes with Gaussian processes[END_REF][START_REF] Woods | Design of experiments for screening[END_REF].

In this paper, following the latter technique, we propose a rigorous and robust method for building a Gp metamodel with a high-dimensional vector of inputs. Moreover, concerning the practical use of this metamodel, our final goal is twofold: to be able to interpret the relationships between the model inputs and outputs with a quantitative sensitivity analysis, and to have a reliable high-level quantile estimation method which does not require additional runs of the computer code.

In what follows, the system under study is generically denoted

Y = g (X 1 , . . . , X d ) (1)
where g(•) is the numerical model (also called the computer code), whose output Y and input parameters X 1 , . . . , X d belong to some measurable spaces Y and X 1 , . . . , X d respectively. X = (X 1 , . . . , X d ) is the input vector and we suppose that X = d k=1 X k ⊂ R d and Y ⊂ R. For a given value of the vector of inputs x = (x 1 , . . . , x d ) ∈ R d , a simulation run of the code yields an observed value y = g(x).

To meet our aforementioned objectives, we propose a sequential methodology which combines several relevant statistical techniques. Our approach consists in four steps: The obtained sample of inputs/outputs constitutes the learning sample. The goal is here to explore, the most efficiently and with a reasonable simulation budget, the domain of uncertain parameters X and get as much information as possible about the behavior of the simulator output Y . For this, we use a space-filling design (SFD) of experiments, providing a full coverage of the high-dimensional input space [START_REF] Fang | Design and modeling for computer experiments[END_REF][START_REF] Woods | Design of experiments for screening[END_REF].

2.

Step 2: Preliminary screening. From the learning sample, a screening technique is performed in order to identify the Primary Influential Inputs (PII) on the model output variability and rank them by decreasing influence. To achieve it, we use dependence measures with associated statistical tests. These measures have several advantages: they can be directly estimated on a SFD, the sensitivity indices that they provide are interpretable quantitatively and their efficiency for screening purpose has been recently illustrated by [START_REF] Veiga | Global sensitivity analysis with dependence measures[END_REF][START_REF] Lozzo | New improvements in the use of dependence measures for sensitivity analysis and screening[END_REF][START_REF] Iooss | An efficient methodology for the analysis and metamodeling of computer experiments with large number of inputs[END_REF].

3.

Step 3: Building and validation of joint metamodel. From the learning sample, a metamodel is built to fit the simulator output Y . For this, we propose to use a joint Gp metamodel [START_REF] Marrel | Global sensitivity analysis of stochastic computer models with joint metamodels[END_REF], by considering only the PII as the explanatory inputs while the other inputs (screened as non significantly influential) are considered as a global stochastic (i.e. unknown) input. Moreover, we use a sequential process to build the joint Gp where the ranked PII are successively included as explanatory inputs in the metamodel (ranking from Step 2). At each iteration, a first Gp model [START_REF] Marrel | An efficient methodology for modeling complex computer codes with Gaussian processes[END_REF], only depending on the current set of explanatory inputs, is built to approximate the mean component. The residual effect of the other inputs is captured using a second Gp model, also function of the explanatory inputs, which approximates the variance component. The accuracy and prediction capabilities of the joint metamodel are controlled on a test sample or by cross-validation.

4.

Step 4: Use of the metamodel for sensitivity analysis and uncertainty propagation. A quantitative sensitivity analysis (Step 4A) and an uncertainty propagation (Step 4B) are performed using the joint metamodel instead of the computer model, leading to a large gain of computation time [START_REF] Zabalza | Prediction and density estimation of a horizontal well productivity index using generalized linear models[END_REF][START_REF] Marrel | Global sensitivity analysis of stochastic computer models with joint metamodels[END_REF]. In this work, we are particularly interested in estimating variance-based sensitivity indices (namely Sobol' indices) and the 95%-quantile of model output.

For ease of reading and understanding of the methodology, Figure 1 gives a general workflow of the articulation of the main steps. For each step, the dedicated sections, the main notations that will be properly introduced later and the key equations are also referenced, in order to provide a guideline for the reader. The paper is organized as follows. In the following section, the nuclear test case which constitutes the guideline application of the paper is described. It consists in a thermal-hydraulic calculation case which simulates an accidental scenario in a nuclear Pressurized Water Reactor. Then, each step of the above statistical methodology is detailed in a dedicated section and illustrated as the same time on the use-case. A last section concludes the work.

II. THERMAL-HYDRAULIC TEST CASE

Our use-case consists in thermal-hydraulic computer experiments, typically used in support of regulatory work and nuclear power plant design and operation. Indeed, some safety analysis considers the so-called "Loss Of Coolant Accident" which takes into account a double-ended guillotine break with a specific size piping rupture. The test case under study is a simplified one, as regards both physical phenomena and dimensions of the system, with respect to a realistic modeling of the reactor. The numerical model is based on code CATHARE2 (V2.5 3mod3.1) which simulates the time evolution of physical quantities during a thermal-hydraulic transient. It models a test carried out on the Japanese mock-up "Large Scale Test Facility" (LSTF) in the framework of the OECD/ROSA-2 project, and which is representative of an Intermediate Break Loss Of Coolant Accident (IBLOCA) [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]. This mock-up represents a reduced scale Westinghouse PWR (1/1 ratio in height and 1/48 in volume), with two loops instead of the four loops on the actual reactor and an electric powered heating core (10 MWe), see Figure 2. It operates at the same pressure and temperature values as the reference PWR. The simulated accidental transient is an IBLOCA with a break on the cold leg and no safety injection on the broken leg. The test under study reproduces a PWR 17% (of cold leg cross-sectional area) cold leg IBLOCA transient with total failure of the auxiliary feedwater, single failure of diesel generators and three systems only available in the intact loop (high pressure injection, accumulator and low pressure injection).

CATHARE2 is used to simulate this integral effect test (see [START_REF] Geiser | The use of integral effect tests for the justification of new evaluation models based on the BEPU approach[END_REF] for the full details of the modeling). During an IBLOCA, the reactor coolant system minimum mass inventory and the Peak Cladding Temperature (PCT) are obtained shortly after the beginning of the accumulators' injection. Figure 3 shows the CATHARE2 prediction and the experimental values of the maximal cladding temperature (also called maximal heater rod temperature) obtained during the test. The conclusion of [START_REF] Geiser | The use of integral effect tests for the justification of new evaluation models based on the BEPU approach[END_REF], which also presents other results, is that the CATHARE2 modeling of the LSTF allows to reproduce the global trends of the different physical phenomena during the transient of the experimental test. In our study, the main output variable of interest will be a single scalar which is the PCT during the accident transient (as an example, see the peak in Figure 3). This quantity is derived from the physical outputs provided by CATHARE2 code.

The input parameters of the CATHARE2 code correspond to various system parameters as boundary conditions, some critical flow rates, interfacial friction coefficients, condensation coeffi- I gives more details about these uncertain inputs and their probability density functions (pdf). The nature of these uncertainties appears to be epistemic since they come from a lack of knowledge on the true value of these parameters. simulation is too important to directly perform all the statistical analysis which are required in a safety study and for which many simulations are needed. To overcome this limitation, an accurate metamodel, built from a reduced number of direct code calculations, will make it possible to develop a more complete and robust safety demonstration.

III. STEP 1: DESIGN OF EXPERIMENTS

This initial step of sampling is to define a design of n experiments for the inputs and performing the corresponding runs with the numerical model g. The obtained sample of inputs/outputs will constitute the learning sample on which the metamodel will then be fitted. The objective is therefore to investigate, most efficiently and with few simulations, the whole variation domain of the uncertain parameters in order to build a predictive metamodel which approximates as accurately as possible the output Y .

For this, we propose to use a space-filling design (SFD) of a n of experiments, this kind of design providing a full coverage of the high-dimensional input space [START_REF] Fang | Design and modeling for computer experiments[END_REF]. Among SFD types, a

Latin Hypercube Sample (LHS, [START_REF] Mckay | A comparison of three methods for selecting values of input variables in the analysis of output from a computer code[END_REF])) with optimal space-filling and good projection properties [START_REF] Woods | Design of experiments for screening[END_REF] would be well adapted. In particular, [START_REF] Fang | Design and modeling for computer experiments[END_REF][START_REF] Damblin | Numerical studies of space filling designs: Optimization of Latin hypercube samples and subprojection properties[END_REF] have shown the importance of ensuring good loworder sub-projection properties. Maximum projection designs [START_REF] Joseph | Maximum projection designs for computer experiments[END_REF] or low-centered L 2 discrepancy LHS [START_REF] Jin | An efficient algorithm for constructing optimal design of computer experiments[END_REF] are then particularly well-suited.

Mathematically, the experimental design corresponds to a n-size sample x (1) , . . . , x (n) which is performed on the model (or code) g. This yields n model output values denoted y (1) , . . . , y (n) with

y (i) = g(x (i) ). The obtained learning sample is denoted (X s , Y s ) with X s = x (1) T , . . . , x (n) T T
and Y s = y (1) , . . . , y (n) T . Then, the goal is to build an approximating metamodel of g from the

n-sample (X s , Y s ).
The number n of simulations is a compromise between the CPU time required for each simulation and the number of input parameters. For uncertainty propagation and metamodelbuilding purpose, some rules of thumb propose to choose n at least as large as 10 times the dimension d of the input vector [START_REF] Loeppky | Choosing the sample size of a computer experiment: A practical guide[END_REF][START_REF] Marrel | An efficient methodology for modeling complex computer codes with Gaussian processes[END_REF].

To build the metamodel for the IBLOCA test case, n = 500 CATHARE2 simulations are performed following a space-filling LHS built in dimension d = 27. The histogram of the obtained values for the output of interest, namely the PCT, is given by Figure 4 (temperature is in • C). A kernel density estimator [START_REF] Parzen | On Estimation of a Probability Density Function and Mode[END_REF] of the data is also added on the plot to provide an estimator of the probability density function. A bimodality seems to be present in the histogram. It underlines the existence of bifurcation or threshold effects in the code, probably caused by a phenomenon of counter current flow limitation between the bend of hot legs and the steam generator inlet plena. Remark III.1 Note that the input values are sampled following their prior distributions defined on their variation ranges. Indeed, as we are not ensured to be able to build a sufficiently accurate metamodel, we prefer sample the inputs following the probabilistic distributions in order to have at least a probabilized sample of the uncertain output, on which statistical characteristics could be estimated. Moreover, as explained in the next section, dependence measures can be directly estimated on this sample, providing first usable results of sensitivity analysis.

Finally, the bimodality that has been observed on the PCT distribution strengthens the use of advanced sensitivity indices (i.e. more general than linear ones or variance-based ones) in our subsequent analysis.

IV. STEP 2: PRELIMINARY SCREENING BASED ON DEPENDENCE MEA-SURE

From the learning sample, a screening technique is performed in order to identify the primary influential inputs (PII) on the variability of model output. It has been recently shown that screening based on dependence measures [START_REF] Veiga | Global sensitivity analysis with dependence measures[END_REF][START_REF] Lozzo | New improvements in the use of dependence measures for sensitivity analysis and screening[END_REF] or on derivative-based global sensitivity measures [START_REF] Kucherenko | Derivative-Based Global Sensitivity Measures[END_REF][START_REF] Roustant | Poincaré inequalities on intervals -application to sensitivity analysis[END_REF] are very efficient methods which can be directly applied on a SFD. Moreover, beyond the screening job, these sensitivity indices can be quantitatively interpreted and used to order the PII by decreasing influence, paving the way for a sequential building of metamodel. In the considered IBLOCA test case, the adjoint model is not available and the derivatives of the model output are therefore not computed because of their costs. The screening step will then be based only on dependence measures, more precisely on Hilbert-Schmidt independence criterion (HSIC) indices, directly estimated from the learning sample.

The dependence measures for screening purpose have been proposed by [START_REF] Veiga | Global sensitivity analysis with dependence measures[END_REF] and [START_REF] Lozzo | New improvements in the use of dependence measures for sensitivity analysis and screening[END_REF]. These sensitivity indices are not the classical ones based on the decomposition of output variance (see [START_REF] Iooss | A review on global sensitivity analysis methods[END_REF] for a global review). They consider higher order information about the output behavior in order to provide more detailed information. Among them, the Hilbert-Schmidt independence criterion (HSIC) introduced by [START_REF] Gretton | Measuring statistical dependence with Hilbert-Schmidt norms[END_REF] builds upon kernel-based approaches for detecting dependence, and more particularly on cross-covariance operators in reproducing kernel Hilbert spaces, see Appendix A for mathematical details.

From the estimated R 2 HSIC [START_REF] Veiga | Global sensitivity analysis with dependence measures[END_REF], independence tests can be performed for a screening purpose. The objective is to separate the inputs into two sub-groups, the significant ones and the non-significant ones. For a given input X k , statistical HSIC-based tests aims at testing the null hypothesis "H (k) 0 : X k and Y are independent", against its alternative "H (k)

1 : X k and Y are dependent". The significance level c of these tests is hereinafter noted α. Several HSIC-based statistical tests are available: asymptotic versions based on an approximation with a Gamma law (for large sample size), spectral extensions and permutation-based versions for non-asymptotic case (case of small sample size). All these tests are described and compared in [START_REF] Lozzo | New improvements in the use of dependence measures for sensitivity analysis and screening[END_REF]; a guidance to use them for a screening purpose is also proposed.

So, at the end of this HSIC-based screening step, the inputs are clustered into two subgroups, PII and non-influential inputs, and the PII are ordered by decreasing R 2 HSIC .This order will be used for the sequential metamodel building in step 3.

On the IBLOCA test case, from the learning sample of n = 500 simulations, R mary Influential Inputs). These inputs are given by Table II with their estimated R 2 HSIC . Ordering them by decreasing R 2 HSIC reveals:

• the large influence of the interfacial friction coefficient in the horizontal part of the hot legs (X 10 ),

• followed by the minimum stable film temperature in the core X 2 , the interfacial friction coefficient in the SG inlet plena X 12 and the wall to liquid friction (in under-saturated break flow conditions) in the break line X 22 ,

• followed by seven parameters with a lower influence: the interfacial friction coefficients in the upper plenum X 15 , the downcomer X 13 , the core X 14 and the SG outlet plena and crossover legs together X 9 , the heat transfer coefficient in the core for film boiling X 5 , the interfacial friction coefficient of the saturated break flow X 26 and the condensation coefficient in the jet during the injection X 27 .

These results clearly underline the predominance influence of the uncertainties on various interfacial friction coefficients. From the learning sample, some scatterplots of the PCT with respect to some well-chosen inputs (the three most influential ones: X 2 , X 10 and X 12 ) are displayed in Figure 5. An additional local regression using weighted linear least squares and a first degree polynomial model (moving average filter) is added on each scatterplot to extract a possible tendency. One can observe that larger values of the interfacial friction coefficient in the horizontal part of the hot legs (X 10 ) lead to larger values of the PCT. This can be explained by the increase of vapor which brings the liquid in the horizontal part of hot legs, leading to a reduction of the liquid water return from the rising part of the U-tubes of the SG to the core (through the hot branches and the upper plenum). Since the amount of liquid water available to the core cooling is reduced, higher PCT are observed.

In addition, we notice a threshold effect concerning this input: beyond a value of 2, the water 13 non-return effect seems to have been reached, and X 10 no longer appears to be influential. We also note that the minimum stable film temperature in the core (X 2 ) shows a trend: the more it increases, the lower the PCT. This is explained by the fact that in the film boiling regime in the core (i.e. when the rods are isolated from the liquid by a film of vapor), X 2 represents (with a decrease in heat flux) the temperature from which the thermal transfer returns to the nucleate boiling regime. Thus, the larger X 2 , the faster the re-wetting of the rods, the faster the cladding temperature excursion is stopped, and thus the lower the PCT. In the next steps, only the eleven PII are considered as explanatory variables in the joint metamodel and will be successively included in the building process. The other sixteen variables will be joined in a so-called uncontrollable parameter.

V. STEP 3: JOINT GP METAMODEL WITH SEQUENTIAL BUILDING PRO-CESS Among all the metamodel-based solutions (polynomials, splines, neural networks, etc.), we focus our attention on the Gaussian process (Gp) regression, which extends the kriging principles of geostatistics to computer experiments by considering the correlation between two responses of a computer code depending on the distance between input variables. The Gp-based metamodel presents some real advantages compared to other metamodels: exact interpolation property, simple analytical formulations of the predictor, availability of the mean squared error of the predictions and the proved capabilities for modeling of numerical simulators (see [START_REF] Sacks | Design and analysis of computer experiments[END_REF], [START_REF] Santner | The design and analysis of computer experiments[END_REF] or [START_REF] Marrel | An efficient methodology for modeling complex computer codes with Gaussian processes[END_REF]). The reader can refer to [START_REF] Rasmussen | Gaussian processes for machine learning[END_REF] for a detailed review on Gp metamodel.

However, for its application to complex industrial problems, developing a robust implementation methodology is required. Indeed, it often implies the estimation of several hyperparameters involved in the covariance function of the Gp (e.g. usual case of anisotropic stationary covariance function). Therefore, some difficulties can arise from the parameter estimation procedure (instability, high number of hyperparameters, see [START_REF] Marrel | An efficient methodology for modeling complex computer codes with Gaussian processes[END_REF] for example). To tackle this issue, we propose a progressive estimation procedure based on the result of the previous screening step and using a joint Gp approach [START_REF] Marrel | Global sensitivity analysis of stochastic computer models with joint metamodels[END_REF]. The interest of the previous screening step becomes twofold. First, the input space, on which each component of the joint Gp is built, can be reduced to the PII space (only the PII are explanatory inputs of Gp). Secondly, the joint Gp is built with a sequential process where the ranked PII are successively included as explanatory inputs in the metamodel. It is expected that these two uses of screening results could significantly make the joint Gp building easier and more efficient.

V.A. Sequential Gp-building process based on successive inclusion of PII as explanatory variables

At the end of the screening step, the PII are ordered by decreasing influence (decreasing R 2 HSIC ). They are successively included as explanatory inputs in the Gp metamodel while the other inputs (the remaining PII and the other non-PII inputs) are joined in a single macro-parameter which is considered as an uncontrollable parameter (i.e. a stochastic parameter, notion detailed in section V.B). Thus, at the j th iteration, a joint Gp metamodel is built with, as explanatory inputs, the j first ordered PII. The definition and building procedure of a joint Gp is fully described in [START_REF] Marrel | Global sensitivity analysis of stochastic computer models with joint metamodels[END_REF] and summarized in the next subsection.

However, building a Gp or a joint Gp involves to perform a numerical optimization in order to estimate all the parameters of the metamodel (covariance hyperparameters and variance parameter). As we usually consider in computer experiments anisotropic (stationary) covariance, the number of hyperparameters linearly increases with the number of inputs. In order to improve the robustness of the optimization process and deal with a large number of inputs, the estimated hyperparameters obtained at the (j -1) th iteration are used, as starting points for the optimization algorithm. This procedure is repeated until the inclusion of all the PII. Note that this sequential estimation process is directly adapted from the one proposed by [START_REF] Marrel | An efficient methodology for modeling complex computer codes with Gaussian processes[END_REF].

V.B. Joint Gp metamodeling

We propose to use a joint Gp metamodeling to handle the group of non-PII inputs and capture their residual effect. In the framework of stochastic computer codes, [START_REF] Zabalza | Prediction and density estimation of a horizontal well productivity index using generalized linear models[END_REF] proposed to model the mean and dispersion of the code output by two interlinked Generalized Linear Models (GLM), called "joint GLM". This approach has been extended by [START_REF] Marrel | Global sensitivity analysis of stochastic computer models with joint metamodels[END_REF] to several nonparametric models and the best results on several tests are obtained with two interlinked Gp models, called "joint Gp". In this case, the stochastic input is considered as an uncontrollable parameter denoted X ε (i.e. governed by a seed variable).

We extend this approach to a group of non-explanatory variables. More precisely, the input variables X = (X 1 , . . . , X d ) are divided in two subgroups: the explanatory ones denoted X exp and the others denoted X ε . The output is thus defined by Y = g(X exp , X ε ) and the metamodeling process will now focus on fitting the random variable Y |X exp d . Under this hypothesis, the joint metamodeling approach yields building two metamodels, one for the mean Y m and another for the

dispersion component Y d : Y m (X exp ) = E(Y |X exp ) (2)
Y d (X exp ) = Var(Y |X exp ) = E (Y -Y m (X exp )) 2 |X exp . (3)
where E[Z] is the usual notation for the expected (i.e. mean) value of a random variable Z.

To fit these mean and dispersion components, we propose to use the methodology proposed by [START_REF] Marrel | Global sensitivity analysis of stochastic computer models with joint metamodels[END_REF]. To summarize, it consists in the following steps. First, an initial Gp denoted Gp m,1 is d Y |Xexp (i.e. Y knowing Xexp) is a random variable as its value depends on the uncontrollable random variable Xε.

estimated for the mean component with homoscedastic nugget effect e . A nugget effect is required to relax the interpolation property of the Gp metamodel. Indeed, this property, which would yield zero residuals for the whole learning sample, is no longer desirable as the output Y |X exp is stochastic. Then, a second Gp, denoted Gp v,1 , is built for the dispersion component with, here also, an homoscedastic nugget effect. Gp v,1 is fitted on the squared residuals from the predictor of Gp m,1 . Its predictor is considered as an estimator of the dispersion component. The predictor of Gp v,1 provides an estimation of the dispersion at each point, which is considered as the value of the heteroscedastic nugget effect. The homoscedastic hypothesis is so removed and a new Gp, denoted Gp m,2 , is fitted on data, with the estimated heteroscedastic nugget. This nugget, as a function of X exp , accounts for a potential interaction between X exp and the uncontrollable parameter X ε . Finally, the Gp on the dispersion component is updated from Gp m,2 following the same methodology as for Gp v,1 .

Remark V.1 Note that some parametric choices are made for all the Gp metamodels: a constant trend and a Matérn stationary anisotropic covariance are chosen. All the hyperparameters (covariance parameters) and the nugget effect (when homoscedastic hypothesis is done) are estimated by maximum likelihood optimization process.

V.C. Assessment of metamodel accuracy

To evaluate the accuracy of a metamodel, we use the predictivity coefficient Q 2 :

Q 2 = 1 - ntest i=1 y (i) -ŷ(i) 2 ntest i=1 y (i) -1 ntest ntest i=1 y (i) 2 (4)
where (x (i) ) 1≤i≤ntest is a test sample, (y (i) ) 1≤i≤ntest are the corresponding observed outputs and (ŷ (i) ) 1≤i≤ntest are the metamodel predictions. Q 2 corresponds to the coefficient of determination in prediction and can be computed on a test sample independent from the learning sample or by cross-validation on the learning sample. The closer to one the Q 2 , the better the accuracy e Borrowed from geostatistics to refer to an unexpected nugget of gold found in a mining process, a constant nugget effect assumes an additive white noise effect, whose variance constitutes the nugget parameter. Most often, this variance is assumed to be constant, independent from the inputs (here Xexp), and the nugget effect is called homoscedastic. When this variance depends on the value of x (i.e. is a function of X), the nugget effect is called heteroscedastic.

of the metamodel. In the framework of joint Gp-modeling, Q 2 criterion will be used to assess the accuracy of the mean part of the joint Gp, namely Gp m,• , whether Gp m,• is a homoscedastic (Gp m,1 ) or heteroscedastic Gp (Gp m,2 ). This quantitative information can be supplemented with a plot of predicted against observed values (ŷ (i) with respect to y (i) ) or a quantile-quantile plot.

To evaluate the quality of the dispersion part of a joint metamodel (denoted Gp v,• ), we use the graphical tool introduced in [START_REF] Marrel | Global sensitivity analysis of stochastic computer models with joint metamodels[END_REF] to assess the accuracy of the confidence intervals predicted by a Gp metamodel. For a given Gp metamodel, It consists in evaluating the proportions of observations that lie within the α-theoretical confidence intervals which are built with the mean squared error of Gp predictions (the whole Gp structure is used and not only the conditional mean). These proportions (i.e. the observed confidence intervals) can be visualized against the α-theoretical confidence intervals, for different values of α.

V.D. Application on IBLOCA test case

The joint Gp metamodel is built from the learning sample of n = 500: the eleven PII identified at the end of the the screening step are considered as the explanatory variables while the sixteen others are considered as the uncontrollable parameter. Gps on mean and dispersion components are built using the sequential building process described in section V.A where PII ordered by decreasing R 2 HSIC are successively included in Gp. Q 2 coefficient of mean component Gp m is computed by cross validation at each iteration of the sequential building process. The results which are given by Table III show an increasing predictivity until its stabilization around 0.87, which illustrates the robustness of the Gp building process. The first four PII make the major contribution yielding a Q 2 around 0.8, the four following ones yield minor improvements (increase of 0.02 on average for each input) while the three last PII does not improve the Gp predictivity. Note that, in this application, these results remain unchanged whether we consider homoscedastic Gp (Gp m,1 ) or heteroscedastic Gp (Gp m,2 ), the heteroscedastic nugget effect not significantly improving the Gp predictor for the mean component. Thus, only 13% of the output variability remains here not explained by Gp m , this includes both the inaccuracy of the Gp m (part of Y m not fitted by Gp) and the total effect of the uncontrollable parameter, i.e. the group of non-selected inputs.

For this exercise, 600 remaining CATHARE2 simulations, different from the learning sample,

TABLE III

Evolution of Gp m metamodel predictivity during the sequential process building, for each new additional PII. PII X 10 X 2 X 12 X 22 X 15 X 13 X 9 X 5 X 14 X 26 X 27 Q 2 0.60 0.64 0.70 0.79 0.81 0.83 0.85 0.85 0.87 0.87 0.87 are also available. As they are not used to build the Gp metamodel, this set of simulations will constitutes a test sample. The Q 2 computed on this test sample is Q 2 = 0.90 for both Gp m,1 and Gp m,2 , which is consistent with the estimations by cross-validation. Now, to assess the quality of dispersion part of the joint Gp, the predicted confidence intervals are compared with the theoretical ones (cf. Section V.C). Figure 6 gives the results obtained by Gp m,1 and Gp m,2 on the learning sample (by cross-validation) and the test sample, since available here. It clearly illustrates both the interest of considering a heteroscedatic nugget effect and the efficiency of using a joint Gp model to fit and predict this nugget. It can be seen that the joint Gp yields the most accurate confidence intervals in prediction, especially for the test sample. Indeed, all its points are close to the theoretical y = x line (a deviation is only observed for the learning sample for the highest α), while the homoscedastic Gp tends to give too large confidence intervals. Thus, in this case, the heteroscedasticity hypothesis is justified and, consequently, the proposed joint Gp model is clearly more competitive than the simple one.

VI. STEP 4A: VARIANCE-BASED SENSITIVITY ANALYSIS

Sensitivity Analysis methods allow to answer the question "How do the input parameters variations contribute, qualitatively or quantitatively, to the variation of the output?" [START_REF] Saltelli | Global sensitivity analysis -The primer[END_REF]. These tools can detect non-significant input parameters in a screening context, determine the most significant ones, measure their respective contributions to the output or identify an interaction between several inputs which impacts strongly the model output. From this, engineers can guide the characterization of the model by reducing the output uncertainty: for instance, they can calibrate the most influential inputs and fix the non-influential ones to nominal values. Many surveys on sensitivity analysis exist in the literature, such as [START_REF] Kleijnen | Sensitivity analysis and related analyses: a review of some statistical techniques[END_REF], [START_REF] Frey | Identification and review of sensitivity analysis methods[END_REF] or [START_REF] Helton | Survey of sampling-based methods for uncertainty and sensitivity analysis[END_REF]. Sensitivity analysis can be divided into two sub-domains: the local sensitivity analysis and the global sensitivity analysis. The first one studies the effects of small input perturbations around nominal values on the model output [START_REF] Cacuci | Sensitivity theory for nonlinear systems. I. Nonlinear functional analysis approach[END_REF] Fig. 6. Proportion of observations that lie within the α-confidence interval predicted by the Gp, according to the theoretical α. Top: results for homoscedastic Gp (Gp m,1 ) on the learning sample (left) and on the test sample (right). Bottom: results for heteroscedastic Gp (Gp m,2 ) on the learning sample (left) and on the test sample (right). while the second one considers the impact of the input uncertainty on the output over the whole variation domain of uncertain inputs [START_REF] Saltelli | Global sensitivity analysis -The primer[END_REF]. We focus here on one of the most widely used global sensitivity indices, namely Sobol' indices which are based on output variance decomposition.

A classical approach in global sensitivity analysis is to compute the first-order and total Sobol' indices which are based on the output variance decomposition [START_REF] Sobol | Sensitivity estimates for non linear mathematical models[END_REF][START_REF] Homma | Importance measures in global sensitivity analysis of non linear models[END_REF], see Appendix B for mathematical details. Sobol' indices are widely used in global sensitivity analysis because they are easy to interpret and directly usable in a dimension reduction approach. However, their estimation (based on Monte-Carlo methods for example) requires a large number of model evaluations, which is intractable for time expensive computer codes. An classical alternative solution consists in using a metamodel to compute these indices. Note that a connection can be made between the estimation error of Sobol' indices when using a metamodel and the predictivity coefficient Q 2 of the metamodel. Indeed, when the Q 2 is estimated on a probabilized sample of the inputs (in other words when it is computed under the probability distribution of the inputs), it provides an estimation of the part of variance unexplained by the metamodel. This can be kept in mind when interpreting the Sobol' indices estimated with the metamodel.

VI.A. Sobol' indices with a joint Gp metamodel

In the case where a joint Gp metamodel is used to take into account an uncontrollable input X ε , [START_REF] Iooss | Global sensitivity analysis of computer models with functional inputs[END_REF] and [START_REF] Marrel | Global sensitivity analysis of stochastic computer models with joint metamodels[END_REF] have shown how to deduce Sobol' indices from this joint metamodel, see Appendix C for mathematical details.

Therefore, from a joint Gp, it is only possible to estimate Sobol' indices of any subset of X exp (equation ( 12)) and the total Sobol' index of X ε (equation ( 13)). The latter is interpreted as the total sensitivity index of the uncontrollable process. The individual index of X ε or any interaction index involving X ε are not directly reachable from joint Gp; their contributions in S T ε can not be distinguished. This constitutes a limitation of this approach. However, the potential interactions between X ε and inputs of X exp could be pointed out, considering all the primary and total effects of all the other parameters. The sensitivity analysis of Y d can also be a relevant indicator: if a subset X u of X exp is not influential on Y d , we can deduce that S uε is equal to zero. Note that in practice, Var(Y ) which appears in both equations ( 12) and ( 13) can be estimated directly from the learning sample (empirical estimator) or from the fitted joint Gp, using equation [START_REF] Petruzzi | Approaches, relevant topics, and internal method for uncertainty evaluation in predictions of thermal-hydraulic system codes[END_REF].

VI.B. Results on IBLOCA test case

From the joint Gp built in section V.D, Sobol' indices of PII are estimated from Gp m by using equation [START_REF] Fang | Design and modeling for computer experiments[END_REF], Var(Y ) being estimated with Gp m and Gp d using equation [START_REF] Petruzzi | Approaches, relevant topics, and internal method for uncertainty evaluation in predictions of thermal-hydraulic system codes[END_REF]. For this, intensive Monte Carlo methods are used (see e.g. the pick-and-freeze estimator of [START_REF] Gamboa | Statistical inference for Sobol pick freeze Monte Carlo methods[END_REF]): tens of thousands simulations of Gp are done. Remind that predictions of Gp are very inexpensive (few seconds for several thousand simulations), especially with respect to the thermal-hydraulic simulator. The obtained first Sobol' indices of PII are given by Table IV and represent 85% of the total variance of the output. Qualitative results of HSIC indices are confirmed and refined: X 10 remains the major influential input with 59% of explained variance, followed to a lesser extend by X 12 and X 22 with for each of them 8% of variance. The partial total Sobol' indices involving only PII and derived from Gp m show that additional 4% of variance is due to interaction between X 10 , X 12 and X 22 .

The related second order Sobol' indices are estimated and the significant ones are given in Table IV. The other PII have negligible influence. In short, the set of PII explain a total 89% of the output variance, of which 79% is only due to X 10 , X 12 and X 22 .

TABLE IV

First and second Sobol' indices of PII (in %), estimated with Gp m of the joint Gp metamodel.

PII X 10 X 2 X 12 X 22 X 15 X 13 X 9 X 5 X 14 X 26 X 27 1 st -order Sobol' indices 59 3 8 8 2 1 2 0 2 0 0
Interaction between PII X 10 × X 22 X 10 × X 12 2 nd -order Sobol' indices 3 1

For the physical interpretation, these results confirm those revealed in Section IV, with a rigorous quantification of inputs' importance: the interfacial friction coefficient in the horizontal part of the hot legs (X 10 ) is the main contributor to the uncertainty of the PCT. Moreover, some results have not been revealed by the HSIC-based screening analysis of Table II. At present, Sobol'

indices clearly indicate that the interfacial friction coefficient in the SG inlet plena X 12 and the wall to liquid friction (in under-saturated break flow conditions) in the break line X 22 are more influential than the minimum stable film temperature in the core X 2 . X 22 has a significant influence on the PCT because its low values lead to higher break flow rates, resulting in a loss of the primary water mass inventory at the higher break, and thus a more significant core uncovering (then higher PCT). For X 12 , its higher values lead to a greater supply (by the vapor) of liquid possibly stored in the water plena to the rest of the primary loops (then lower PCT). Table IV also shows that there are some small interaction effects (possibly antagonistic) between X 10 and X 12 , as well as between X 10 and X 22 . Let us remark that deepening this question (which is outside the scope of this paper) would be possible via plotting, from the Gp metamodel, the conditional expectations of the PCT as a function of the interaction variables.

At present, from Gp d and using equation ( 13), the total effect of the group of the sixteen non-PII inputs (i.e. X ε ) is estimated to 9.7%. This includes its effect alone and in interaction with the PII. To further investigate these interactions, Sobol indices of Y d are estimated and HSIC-based statistical dependence tests are applied on Y d . They reveal that only X 10 , X 14 , X 2 , X 22 and X 4

potentially interact with the group of non-PII inputs. At this stage of our analysis, this result cannot be physically interpreted.

VII. STEP 4B: QUANTILE ESTIMATION

As already said in the introduction, once a predictive Gp metamodel has been built, it can be used to perform uncertainty propagation and in particular, estimate probabilities or, as here, quantiles.

VII.A. Gp-based quantile estimation

The most trivial and intuitive approach to estimate a quantile with a Gp metamodel is to apply the quantile definition to the predictor of the metamodel. This direct approach yields a so called plug-in estimator. More precisely, with a Gp metamodel, this approach consists in using only the predictor of the Gp metamodel (i.e. the expectation conditional to the learning sample) in order to estimate the quantile. As the expectation of the Gp mean is a deterministic function of the inputs, this provides a deterministic expression of the quantile but no confidence intervals are available. Moreover, for high (resp. low) quantiles, this methods tends to substantially underestimate (resp. overestimate) the true quantile because the metamodel predictor is usually constructed by smoothing the computer model output values (see an analysis of this phenomenon in [START_REF] Cannamela | Controlled stratification for quantile estimation[END_REF]).

To overcome this problem, [START_REF] Oakley | Estimating percentiles of uncertain computer code outputs[END_REF] has proposed to take advantage of the probabilistic-type Gp metamodel by using its entire structure: not only the mean of the conditional Gp metamodel but also its covariance structure are taken into account. In this full-Gp based approach also called modular Bayesian approach, the quantile definition is therefore applied to the global Gp metamodel and leads to a random variable. The expectation of this random variable can be then considered as a quantile estimator. Its variance and, more generally, all its distribution can then be used as an indicator of the accuracy of the quantile estimate. Confidence intervals can be deduced, which constitutes a great advantage of this full-Gp based approach. Moreover, the efficiency of this approach for high quantile (of the order of 95%) has been illustrated by [START_REF] Rutherford | A response-modeling alternative to surrogate models for support in computational analyses[END_REF].

In practice, the estimation of quantile with the full-Gp approach is based on stochastic simulations (conditionally to the learning sample) of the Gp metamodel, by using the standard method of conditional simulations [START_REF] Chilès | Geostatistics: Modeling spatial uncertainty[END_REF]. We recall just that, to generate conditional Gaussian simulations on a given interval, a discretization of the interval is first considered, then the vector of the conditional expectation and the matrix of the conditional covariance on the discretized interval are required. Note that [START_REF] Gratiet | A Bayesian approach for global sensitivity analysis of (multifidelity) computer codes[END_REF] uses this approach with Gp conditional simulations for the estimation of Sobol' indices and their associated confidence intervals.

In this paper, from our joint Gp model, we will compare the full-Gp approach applied to either the homoscedastic Gp or the heteroscedastic Gp and will proceed as follows:

• For the homoscedastic Gp, the standard technique of conditional simulations is applied to Gp m,1 (built to estimate Y m , with a constant nugget effect).

• For the heteroscedastic Gp, we propose a new technique to simulate the conditional Gp trajectories:

1. The heteroscedastic Gp built for Y m (namely Gp m,2 ) provides the conditional expectation vector and a preliminary conditional covariance matrix.

2. The Gp built for Y d (namely Gp d,2 ) is predicted and provides the heteroscedastic nugget effect which is added to the diagonal of the conditional covariance matrix of the previous step.

3. Conditional simulations are then done using the standard method [START_REF] Chilès | Geostatistics: Modeling spatial uncertainty[END_REF].

VII.B. Results on IBLOCA test case

In this section, we focus on the estimation of the 95%-quantile of the PCT for the IBLOCA application. From the learning sample of size n = 500 and the joint Gp, we compare here the following approaches to estimate the PCT quantile:

• The classical empirical quantile estimator, denoted qemp 95 . A bootstrap method (see for example [START_REF] Zio | Bootstrap and Order Statistics for Quantifying Thermal-Hydraulic Code Uncertainties in the Estimation of Safety Margins[END_REF]) makes it possible to obtain in addition a 90%-confidence interval for this empirical quantile.

• The plug-in estimators from the homoscedastic or the heteroscedastic Gp, denoted qGpm,1-pi 95 and qGpm,2-pi

95

. No confidence intervals are obtained using this estimation method.

• The full-Gp estimators from the homoscedastic or the heteroscedastic Gp, denoted qGpm,1-full 95 and qGpm,2-full 95 respectively. As explained in the previous section, confidence intervals can be deduced with this full-Gp approach.

Table V synthesizes all the results obtained for the PCT 95%-quantile with these different approaches given above. In addition, 90%-confidence intervals are given when they are available.

As explained in Section V.D, we also have 600 other CATHARE2 simulations, for a total of 1100 PCT simulations (learning sample plus test sample). A reference value of the 95%-quantile is obtained from this full sample with the classical empirical quantile estimator:

qref 95 = 742.28 • C .
The empirical estimator based on the learning sample is imprecise but its bootstrap-based confidence interval is consistent with the reference value. As expected, the plug-in approaches provide poor estimations of the quantile, which is here strongly underestimated. The full-Gp approach based on the homoscedastic assumption overestimates the quantile; this is consistent with the analysis of Gp confidence intervals in Figure 6 (too large confidence intervals provided by homoscedastic Gp). Finally, the best result is obtained with the conditional simulation method based on the heteroscedastic Gp metamodel, built with the joint Gp method. This approach yields a more accurate prediction than the usual homoscedastic Gp and outperforms the empirical estimator in terms of confidence interval. Once again, the heteroscedasticity hypothesis is clearly relevant in this case. To cope with this limitation, we propose a methodology mainly based on a predictive joint Gp metamodel, built with an efficient sequential algorithm. First, an initial screening step based on advanced dependence measures and associated statistical tests enabled to identify a group of significant inputs, allowing a reduction of the dimension. The efforts of optimization when fitting the metamodel can then be concentrated on the main influential inputs. The robustness of metamodeling is thus increased. Moreover, thanks to the joint metamodel approach, the nonselected inputs are not completely removed: the residual uncertainty due to dimension reduction is integrated in the metamodel and the global influence of non-selected inputs is so controlled. From this joint Gp metamodel, accurate uncertainty propagation and quantitative sensitivity analysis, not feasible with the numerical model because of its computational cost, become accessible. Thus, the uncertainties of model inputs are propagated inside the joint Gp to estimate Sobol' indices, failure probabilities and/or quantiles, without additional code simulations.

Thus, on the IBLOCA application, a predictive Gp metamodel is built with only a few hundred of code simulations (500 code simulations for 27 uncertain inputs). From this joint Gp, a quantitative sensitivity analysis based on variance decomposition is performed without additional code simulation: Sobol' indices are computed and reveal that the output is mainly explained by four uncertain inputs. One of them, namely the interfacial friction coefficient in the hot legs, is strongly influential with around 60% of output variance explained, the three others being of minor influence. The quite less influence of all the other inputs is also confirmed. Note that a direct and accurate computation of Sobol' indices with the thermal-hydraulic code would have required tens of thousands of code simulations.

The physical interpretation of the results obtained with the screening step and the variancebased sensitivity analysis step are known to be useful for modelers and engineers. This study has demonstrated this once again, in the particular case of an IBLOCA safety analysis, by revealing some physical effects on the PCT of influential inputs which cannot be understood without a global statistical approach (e.g. the threshold effect due to the interfacial friction coefficient in the horizontal part of the hot legs). [START_REF] Geiser | Sensitivity analysis applied to LOCA integral effect tests for the justification of the BEPU approach[END_REF] has also underlined the importance of sensitivity analysis in a validation methodology in order to identify relevant physical model uncertainties on which safety engineers must focus their efforts. Counter-intuitive effects are also present during an IBLOCA transient and only a global understanding of physical phenomena can help. As a perspective of the present work, extending the screening and sensitivity analysis tools to a joint analysis of several relevant output variables of interest (as the PCT time, the primary pressure and the core differential pressure) would be essential.

In the IBLOCA test case, we are particularly interested by the estimation of the 95%quantile of the model output temperature. In nuclear safety, as in other engineering domains, methods of conservative computation of quantiles [START_REF] Nutt | Evaluation of nuclear safety from the outputs of computer codes in the presence of uncertainties[END_REF][START_REF] Hessling | Robustness of Wilks' conservative estimate of confidence intervals[END_REF] have been largely studied. We have shown in the present work how to use and simulate the joint Gp metamodel to reach the same objective: the uncertainty of the influential inputs are directly and accurately propagated through the mean component of the joint metamodel while a confidence bound is derived from the dispersion component in order to take into account the residual uncertainty of the other inputs. Results on the IBLOCA test case show that joint Gp provides a more accurate estimation of the 95%-quantile than the empirical approach, at equal calculation budget. Besides, this estimation is very close from the reference value obtained with 600 additional code simulations. Furthermore, the interest of the heteroscedastic approach in joint Gp is also emphasized: the estimated quantile and the associated confidence interval are much better than those of the homoscedastic approach.

As a future application of modern statistical methods on IBLOCA safety issues, one should mention the use of Gp metamodels to identify penalizing thermal-hydraulic transients, with respect to some particular scenario inputs. As in the present work, strong difficulties are raised by the cpu time cost of the computer code and the large number of uncertain (and uncontrollable) inputs.

1 .

 1 Step 1: Design of experiments. Knowing the variation domain of the input variables, a design of n numerical experiments is firstly performed and yields n model output values.

Fig. 1 .

 1 Fig. 1. General workflow of the statitical methodology.

Fig. 2 .

 2 Fig. 2. Large Scale Test Facility showing the main components and the hot and cold legs.

Fig. 3 .

 3 Fig. 3. Experimental values and physical simulation output of the CATHARE2 model: maximal rod cladding temperature during the transient.

Fig. 4 .

 4 Fig. 4. Histogram of the PCT from the learning sample of n = 500 simulations.

Fig. 5 .

 5 Fig. 5. Scatterplots with local polynomial regression of PCT according to several inputs, from the learning sample of n = 500 simulations.

TABLE I

 I List of the 27 uncertain input parameters and associated physical models in CATHARE2 code.

	Type of inputs	Inputs pdf a Physical models
	Heat transfer	X 1	N	Departure from nucleate boiling
	in the core	X 2	U	Minimum film stable temperature
		X 3	LN	HTC b for steam convection
		X 4	LN	Wall-fluid HTC
		X 5	N	HTC for film boiling
	Heat transfer in the steam	X 6	LU	HTC forced wall-steam convection
	generators (SG) U-tube	X 7	N	Liquid-interface HTC for film condensation
	Wall-steam friction in core	X 8	LU
	Interfacial friction	X 9	LN	SG outlet plena and crossover legs together
		X 10	LN	Hot legs (horizontal part)
		X 11	LN	Bend of the hot legs
		X 12	LN	SG inlet plena
		X 13	LN	Downcomer
		X 14	LN	Core
		X 15	LN	Upper plenum
		X 16	LN	Lower plenum
		X 17	LN	Upper head
	Condensation	X 18	LN	Downcomer
		X 19	U	Cold leg (intact)
		X 20	U	Cold leg (broken)
		X 27	U	Jet
	Break flow	X 21	LN	Flashing (undersaturated)
		X 22	N	Wall-liquid friction (undersaturated)
		X 23	N	Flashing delay (undersaturated)
		X 24	LN	Flashing (saturated)
		X 25	N	Wall-liquid friction (saturated)
		X 26	LN	Global interfacial friction (saturated)
	Our objective with this use-case is to provide a good metamodel for sensitivity analysis,
	uncertainty propagation and, more generally, safety studies. Indeed, the cpu-time cost of one
	a U, LU, N and LN respectively stands for Uniform, Log-Uniform, Normal and Log-Normal probability distri-
	butions.

b Heat Transfer Coefficient.

TABLE V

 V In the framework of the estimation of safety margins in nuclear accident analysis, it is essential to quantitatively assess the uncertainties tainting the results of Best-estimate codes. In this context, this paper has been focused on an advanced statistical methodology for Best Estimate Plus Uncertainty (BEPU) analysis, illustrated by a high dimensional thermal-hydraulic test case simulating accidental scenario in a Pressurized Water Reactor (IBLOCA test case). Some statistical analyses such as the estimation of high-level quantiles or quantitative sensitivity analysis (e.g., estimation of Sobol' indices based on variance decomposition) may call in practice for several thousand of code simulations. Complex computer codes, as the ones used in thermal-hydraulic accident scenario simulations, are often too cpu-time expensive to be directly used to perform these studies.

	Results for the 95%-quantile estimates of the PCT (in • C) with its 90%-confidence interval (CI)
	when available.
		Reference f	Empirical g	Plug-in g,h homo-Gp hetero-Gp	Full-Gp g,h (conditional simulations) homo-Gp hetero-Gp
		qref 95	qemp 95	qGpm,1-pi 95	qGpm,2-pi 95	qGpm,1-full 95	qGpm,2-full 95
	Mean	742.28	746.80	736.26	735.83	747.11	741.46
	90%-CI	-	[736.7; 747.41]	-	-	[742.93; 751.32]	[738.76; 744.17]

f The reference method uses the full sample of size n = 1100. g Empirical and Gp-based methods are applied from the learning sample of size n = 500. h The plug-in and full-Gp estimators are based on either the homoscedastic or heteroscedastic Gp metamodel.

VIII. CONCLUSION

ACKNOWLEDGMENTS

We are grateful to Henri Geiser and Thibault Delage who performed the computations of the CATHARE2 code.

A. APPENDIX: HSIC DEPENDENCE MEASURES

If we consider two reproducing kernel Hilbert spaces F k and G of functions X k → R and Y → R respectively, the crossed-covariance C X k ,Y operator associated to the joint probabilistic distribution of (X k , Y ) is the linear operator defined for every f X k ∈ F k and g Y ∈ G by:

C X k ,Y generalizes the covariance matrix by representing higher order correlations between X k and Y through nonlinear kernels. The HSIC criterion is then defined by the Hilbert-Schmidt norm of the cross-covariance operator:

From this, [START_REF] Veiga | Global sensitivity analysis with dependence measures[END_REF] introduces a normalized version of the HSIC which provides a sensitivity index of

The closer to one the R 2 HSIC,k , the stronger the dependence between X k and Y . In practice, [START_REF] Gretton | Measuring statistical dependence with Hilbert-Schmidt norms[END_REF] propose a Monte Carlo estimator of HSIC(X k , Y ) and a plug-in estimator can be deduced for R 2 HSIC,k . Note that Gaussian kernel functions with empirical estimations of the variance parameter are used in our application (see [START_REF] Gretton | Measuring statistical dependence with Hilbert-Schmidt norms[END_REF] for more details).

B. APPENDIX: SOBOL' INDICES

+∞, we can apply the Hoeffding decomposition to the random variable g(X) [START_REF] Hoeffding | A class of statistics with asymptotically normal distributions[END_REF]:

where

, for all u ⊂ {1, . . . , d}. All the 2 d terms in [START_REF] De Crécy | Uncertainty and sensitivity analysis of the LOFT L2-5 test: Results of the BEMUSE programme[END_REF] have zero mean and are mutually uncorrelated with each other. This decomposition is unique and leads to the Sobol' indices. These are the elements of the g(X) variance decomposition according to the different groups of input parameter interactions in [START_REF] De Crécy | Uncertainty and sensitivity analysis of the LOFT L2-5 test: Results of the BEMUSE programme[END_REF]. More precisely, for each u ⊂ {1, . . . , d}, the first-order and total Sobol' indices of X u are defined by:

S u represents the part of the output variance explained by X u , independently from the other inputs, and S T u is the part of the output variance explained by X u considered separately and in interaction with the other input parameters.

In practice, we are usually interested in the first-order sensitivity indices S 1 , . . . , S d , the total ones S T 1 , . . . , S T d and sometimes in the second-order ones S ij , 1

C. APPENDIX: SOBOL' INDICES FROM A JOINT METAMODEL

By adopting the same notations as in Section V.B where X exp is the vector of all the controllable (or explanatory) inputs and X ε denotes the uncontrollable input (group of non explanatory inputs), we can show that the variance of the output variable Y (X exp , X ε ) can be rewritten and deduced from the two components Y m and Y d (equations ( 2) and ( 3)):

where E X (resp. Var X ) denotes the mean (resp. variance) operator with respect to the probability density function of X. Furthermore, the variance of Y is the sum of the contributions of both all the explanatory inputs in X exp and X ε :

where

Variance of the mean component Y m (X exp ) denoted hereafter Y m can be also decomposed:

As

Similarly, the total sensitivity index of X ε is given by:

Note that, as Y d (X exp ) is a positive random variable, the positivity of S T ε is guaranteed.


Xi Liu

Badong Chen
email: chenbd@mail.xjtu.edu.cn

Bin Xu b Zongze Wu

Paul Honeine

Maximum Correntropy Unscented Filter

Keywords: Unscented Kalman Filter (UKF), Unscented Transformation (UT), Maximum Correntropy Criterion (MCC)

The unscented transformation (UT) is an efficient method to solve the state estimation problem for a non-linear dynamic system, utilizing a derivative-free higher-order approximation by approximating a Gaussian distribution rather than approximating a non-linear function. Applying the UT to a Kalman filter type estimator leads to the well-known unscented Kalman filter (UKF). Although the UKF works very well in Gaussian noises, its performance may deteriorate significantly when the noises are non-Gaussian, especially when the system is disturbed by some heavy-tailed impulsive noises. To improve the robustness of the UKF against impulsive noises, a new filter for nonlinear systems is proposed in this work, namely the maximum correntropy unscented filter (MCUF). In MCUF, the UT is applied to obtain the prior estimates of the state and covariance matrix, and a robust statistical linearization regression based on the maximum correntropy criterion (MCC) is then used to obtain the posterior estimates of the state and covariance matrix. The satisfying performance of the new algorithm is confirmed by two illustrative examples.

Introduction

Estimation problem plays a key role in many fields, including communication, navigation, signal processing, optimal control and so on [START_REF] Dash | Adaptive complex unscented Kalman filter for frequency estimation of time-varying signals[END_REF][START_REF] Li | Fault tolerant navigation method for satellite based on information fusion and unscented Kalman filter[END_REF][START_REF] Li | Unscented Kalman filter over unreliable communication networks with markovian packet dropouts[END_REF][START_REF] Partovibakhsh | Adaptive unscented Kalman filter-based online slip ratio control of wheeled-mobile robot[END_REF]. The Kalman filter (KF) assists in obtaining accurate state estimation for a linear dynamic system, which provides an optimal recursive solution under minimum mean square error (MMSE) criterion [START_REF] Bryson | Applied optimal control: Optimization, estimation and control[END_REF][START_REF] Kalman | A new approach to linear filtering and prediction problems[END_REF][START_REF] Nahi | Estimation theory and applications[END_REF]. Nevertheless, most practical systems are inherently nonlinear, and it is not easy to implement an optimal filter for nonlinear systems. To solve the nonlinear filtering problem, so far many sub-optimal nonlinear extensions of the KF have been developed by using some approximations, among which the extended Kalman filter (EKF) [START_REF] Anderson | Optimal filtering[END_REF] and unscented Kalman filter (UKF) [START_REF] Julier | A new method for the nonlinear transformation of means and covariances in filters and estimators[END_REF] are two widely used ones. As a popular nonlinear extension of KF, the EKF approximates the nonlinear system by its first order linearization and uses the original KF on this approximation. However, the crude approximation may lead to divergence of the filter when the function is highly non-linear. Moreover, the cumbersome derivation of the Jacobian matrices often leads to the implementation difficulties. The UKF is an alternative to the EKF, which approximates the probability distribution of the state by a set of deterministically chosen sigma points and propagates the distribution though the non-linear equations. The UKF does not need to calculate the Jacobian matrices and can obtain a better performance than the EKF. However, the UKF may performs poorly when the system is disturbed by some heavy-tailed non-Gaussian noises, which occur frequently in many real-world applications of engineering. The main reason for this is that the UKF is based on the MMSE criterion and thus exhibits sensitivity to heavy-tailed noises [START_REF] Schick | Robust recursive estimation in the presence of heavy-tailed observation noise[END_REF]. Some methods have been proposed to cope with this problem in the literatures. In particular, the Huber's generalized maximum likelihood methodology is an important one [START_REF] El-Hawary | Robust regression-based EKF for tracking underwater targets[END_REF][START_REF] Huber | Robust estimation of a location parameter[END_REF][START_REF] Karlgaard | Huber-based divided difference filtering[END_REF]; X. [START_REF] Wang | Huber-based unscented filtering and its application to vision-based relative navigation[END_REF], which is a combined minimum 1 and 2 norm estimation method. Furthermore, the statistical linear method instead of the first order linearization was introduced in [START_REF] Karlgaard | Huber-based divided difference filtering[END_REF]. Besides Huber's robust statistics, information theoretic quantities (e.g. entropy, mutual information, divergence, etc.) can also be used as a robust cost for estimation problems (B. [START_REF] Chen | System parameter identification: Information criteria and algorithms[END_REF][START_REF] Principe | Information theoretic learning: Renyis entropy and kernel perspectives[END_REF]. As a localized similarity measure in information theoretic learning (ITL), the correntropy has recently been successfully applied in robust machine learning and non-Gaussian signal processing [START_REF] Bessa | Entropy, and correntropy against minimum square error in offline, and online three-day ahead wind power forecasting[END_REF][START_REF] Chen | Maximum correntropy estimation is a smoothed MAP estimation[END_REF][START_REF] Chen | Steady-state mean-square error analysis for adaptive filtering under the maximum correntropy criterion[END_REF][START_REF] Chen | Generalized correntropy for robust adaptive filtering[END_REF][START_REF] Chen | Recursive robust least squares support vector regression based on maximum correntropy criterion[END_REF][START_REF] He | Robust recognition via information theoretic learning[END_REF][START_REF] He | Robust principal component analysis based on maximum correntropy criterion[END_REF]He, Zheng, & Hu, 2011;[START_REF] Liu | Correntropy: Properties, and applications in non-gaussian signal processing[END_REF][START_REF] Shi | Convex combination of adaptive filters under the maximum correntropy criterion in impulsive interference[END_REF][START_REF] Singh | Using correntropy as a cost function in linear adaptive filters[END_REF][START_REF] Wang | Robust hyperspectral unmixing with correntropy-based metric[END_REF][START_REF] Xu | A pitch detector based on a generalized correlation function[END_REF][START_REF] Zhao | Kernel adaptive filtering with maximum correntropy criterion[END_REF]. The adaptive filtering algorithms under the maximum correntropy criterion (MCC) can achieve excellent performance in heavy-tailed non-Gaussian noises (B. [START_REF] Chen | Steady-state mean-square error analysis for adaptive filtering under the maximum correntropy criterion[END_REF][START_REF] Chen | Generalized correntropy for robust adaptive filtering[END_REF][START_REF] Principe | Information theoretic learning: Renyis entropy and kernel perspectives[END_REF][START_REF] Shi | Convex combination of adaptive filters under the maximum correntropy criterion in impulsive interference[END_REF][START_REF] Singh | Using correntropy as a cost function in linear adaptive filters[END_REF][START_REF] Wang | Robust hyperspectral unmixing with correntropy-based metric[END_REF][START_REF] Zhao | Kernel adaptive filtering with maximum correntropy criterion[END_REF]. In particular, in a recent work [START_REF] Chen | Maximum correntropy Kalman filter[END_REF], a linear Kalman type filter has been developed under MCC, which can outperform the original KF significantly especially when the disturbance noises are impulsive.

The goal of this paper is to develop an unscented non-linear Kalman type filter based on the MCC, called the maximum correntropy unscented filter (MCUF). In the MCUF, the unscented transformation (UT) is applied to get a prior estimation of the state and covariance matrix, and a statistical linearization regression model based on MCC is used to obtain the posterior state and covariance. The new filter adopts the UT and statistical linear approximation instead of the first order approximation as in EKF to approximate the nonlinearity, and uses the MCC instead of the MMSE to cope with the non-Gaussianity, and hence can achieve desirable performance in highly nonlinear and non-Gaussian conditions. Moreover, the proposed MCUF is suitable for online implementation on account of the retained recursive structure.

The rest of the paper is organized as follows. In Section 2, we briefly introduce the MCC. In Section 3, we derive the MCUF algorithm. In Section 4, we present two illustrative examples and show the desirable performance of the proposed algorithm. Finally, Section 5 concludes this paper.

Maximum correntropy criterion

Correntropy is a generalized similarity measure between two random variables. Given two random variables X, Y ∈ R with joint distribution function F XY (x, y), the correntropy is defined by

V (X, Y ) = E [κ(X, Y )] = κ(x, y)dF XY (x, y) (1)
where E denotes the expectation operator, and κ(•, •) is a shift-invariant Mercer kernel. In this work, without mentioned otherwise, the used kernel function is the Gaussian kernel:

κ(x, y) = G σ (e) = exp - e 2 2σ 2 (2)
where e = x -y, and σ > 0 stands for the kernel bandwidth. In many practical situations, we have only a set of finite data and the joint distribution F XY is unknown. In these cases, one can estimate the correntropy using a sample mean estimator:

V (X, Y ) = 1 N N i=1 G σ (e(i)) (3)
where e(i) = x(i) -y(i), with {x(i), y(i)} N i=1 being N samples drawn from F XY . Using the Taylor series expansion for the Gaussian kernel yields

V (X, Y ) = ∞ n=0 (-1) n 2 n σ 2n n! E (X -Y ) 2n (4)
Thus, the correntropy is a weighted sum of all even order moments of the error variable X -Y . The kernel bandwidth appears as a parameter to weight the second order and higher order moments. With a very large kernel bandwidth (compared to the dynamic range of the data), the correntropy will be dominated by the second order moment. Suppose our goal is to learn a parameter vector W of an adaptive model, and let x(i) and y(i) denote, respectively, the model output and the desired response. The MCC based learning can be formulated as solving the following optimization problem:

W = arg max W ∈Ω 1 N N i=1 G σ (e(i)) (5)
where W denotes the optimal solution, and Ω denotes a feasible set of the parameter.

Maximum correntropy unscented filter

In this section, we propose to combine the MCC and a statistical linear regression model together to derive a novel non-linear filter, which can perform very well in non-Gaussian noises, since correntropy embraces second and higher order moments of the error. Let's consider a nonlinear system described by the following equations:

x(k) = f (k -1, x(k -1)) + q(k -1), (6)
y(k) = h (k, x(k)) + r(k). (7)
where x(k) ∈ R n denotes a n-dimensional state vector at time step k, y(k) ∈ R m represents an mdimensional measurement vector, f is a nonlinear system function, and h is a nonlinear measurement function and both are assumed to be continuously differentiable. The process noise q(k -1) and measurement noise r(k) are generally assumed to meet the independence assumption with zero mean and covariance matrices

E q(k -1)q T (k -1) = Q(k -1), E r(k)r T (k) = R(k) (8)
Similar to other Kalman type filters, the MCUF also includes two steps, namely the time update and measurement update:

Time update

A set of 2n+1 samples, also called sigma points, are generated from the estimated state x(k-1|k-1) and covariance matrix P(k -1|k -1) at the last time step k -1:

χ 0 (k -1|k -1) = x(k -1|k -1), χ i (k -1|k -1) = x(k -1|k -1) + (n + λ)P(k -1|k -1) i , for i = 1 . . . n, χ i (k -1|k -1) = x(k -1|k -1) - (n + λ)P(k -1|k -1) i-n , for i = n + 1 . . . 2n. ( 9
)
where

(n + λ)P(k -1|k -1)
i is the i-th column of the matrix square root of (n + λ)P(k -1|k -1), with n being the state dimension and λ being a composite scaling factor, given by

λ = α 2 (n + φ) -n (10)
where α determines the spread of the sigma points, usually selected as a small positive number, and φ is a parameter that is often set to 3 -n.

The transformed points are then given through the process equation:

χ i * (k|k -1) = f k -1, χ i (k -1|k -1) , for i = 0 . . . 2n (11)
The prior state mean and covariance matrix are thus estimated by

x(k|k -1) = 2n i=0 w i m χ i * (k|k -1), (12)
P(k|k -1) = 2n i=0 w i c χ i * (k|k -1) -x(k|k -1) × χ i * (k|k -1) -x(k|k -1) T + Q(k -1). ( 13
)
in which the corresponding weights of the state and covariance matrix are

w 0 m = λ (n + λ) , w 0 c = λ (n + λ) + (1 -α 2 + β), w i m = w i c = 1 2(n + λ) , for i = 1 . . . 2n. ( 14
)
where β is a parameter related to the prior knowledge of the distribution of x(k) and is set to 2 in the case of the Gaussian distribution.

Measurement update

Similarly, a set of 2n + 1 sigma points are generated from the prior state mean and covariance matrix

χ 0 (k|k -1) = x(k|k -1), χ i (k|k -1) = x(k|k -1) + (n + λ)P(k|k -1) i , for i = 1 . . . n, χ i (k|k -1) = x(k|k -1) - (n + λ)P(k|k -1) i-n , for i = n + 1 . . . 2n. ( 15
)
These points are transformed through the process equation as

γ i (k) = h(k, χ i (k|k -1)), for i = 0 . . . 2n (16)
The prior measurement mean can then be obtained as

y(k) = 2n i=0 w i m γ i (k), (17)
Further, the state-measurement cross-covariance matrix is given by

P xy (k) = 2n i=0 w i c χ i (k|k -1) -x(k|k -1) γ i (k) -y(k) T . (18)
Next, we apply a statistical linear regression model based on the MCC to accomplish the measurement update. First, we formulate the regression model. We denote the prior estimation error of the state by

η(x(k)) = x(k) -x(k|k -1) (19)
and define the measurement slope matrix as

H(k) = P -1 (k|k -1)P xy (k) T ( 20
)
Then the measurement equation ( 7) can be approximated by (X. [START_REF] Wang | Huber-based unscented filtering and its application to vision-based relative navigation[END_REF])

y(k) ≈ y(k) + H(k)(x(k) -x(k|k -1)) + r(k) (21)
Combining (12) ( 17) and ( 21), we obtain the following statistical linear regression model:

x(k|k -1) y(k) -y(k) + H(k) x(k|k -1) = I H(k) x(k) + ξ(k) (22)
where ξ(k) is

ξ(k) = η(x(k)) r(k)
with

Ξ(k) = E ξ(k)ξ T (k) = P(k|k -1) 0 0 R(k) = S p (k|k -1)S T p (k|k -1) 0 0 S r (k)S T r (k) = S(k)S T (k) (23)
Here, S(k) can be obtained by the Cholesky decomposition of Ξ(k). Left multiplying both sides of ( 22) by S -1 (k), the statistical regression model is transformed to

D(k) = W(k)x(k) + e(k) (24)
where

D(k) = S -1 (k) x(k|k -1) y(k) -y(k) + H(k) x(k|k -1) , W(k) = S -1 (k) I H(k) , e(k) = S -1 (k)ξ(k).
It can be seen that E e(k)e T (k) = I.

We are now in a position to define a cost function based on the MCC:

J L (x(k)) = L i=1 G σ (d i (k) -w i (k)x(k)) ( 25
)
where

d i (k) is the i-th element of D(k), w i (k) is the i-th row of W(k), and L = n + m is the dimension of D(k).
Under the MCC, the optimal estimate of x(k) arises from the following optimization:

x(k) = arg max x(k) J L (x(k)) = arg max x(k) L i=1 G σ (e i (k)) ( 26
)
where e i (k) is the i-th element of e(k) given by

e i (k) = d i (k) -w i (k)x(k) (27)
The optimal solution of x(k) can be solved through

∂J L (x(k)) ∂x(k) = 0 (28)
It follows easily that

x(k) = L i=1 G σ (e i (k)) w T i (k)w i (k) -1 × L i=1 G σ (e i (k)) w T i (k)d i (k) (29) Since e i (k) = d i (k) -w i (k)x(k)
, the equation ( 29) is actually a fixed-point equation with respect to x(k) and can be rewritten as

x(k) = g (x(k)) (30)
Hence, a fixed-point iterative algorithm can be obtained as (

x(k) t+1 = g ( x(k) t ) (31)
with x(k) t being the estimated state x(k) at the t-th fixed-point iteration .

The fixed-point equation ( 29) in matrix form can also be expressed as

x(k) = W T (k)C(k)W(k) -1 W T (k)C(k)D(k) ( 32
)
where

C(k) = C x (k) 0 0 C y (k) , with C x (k) = diag (G σ (e 1 (k)) , ..., G σ (e n (k))), C y (k) = diag (G σ (e n+1 (k)) , ..., G n+m (e n+m (k))),
which can be further written as (see the Appendix A for a detailed derivation):

x(k) = x(k|k -1) + K(k) (y(k) -y(k)) (33)
where

     K(k) = P(k|k -1)H T (k) H(k)P(k|k -1)H T (k) + R(k) -1 P(k|k -1) = S p (k|k -1)C -1 x (k)S T p (k|k -1) R(k) = S r (k)C -1 y (k)S T r (k) (34)
Meanwhile, the corresponding covariance matrix is updated by

P(k|k) = I -K(k)H(k) P(k|k -1) I -K(k)H(k) T + K(k)R(k)K T (k) (35)
Remark 1: Since K(k) relies on P(k|k -1) and R(k), both related to x(k) through C x (k) and C y (k), respectively, the equation ( 33) is a fixed-point equation of x(k). One can solve (33) by a fixed-point iterative method. The initial value of the fixed-point iteration can be set to x(k|k) 0 = x(k|k -1) or chosen as the least-squares solution x(k|k

) 0 = W T (k)W(k) -1 W T (k)D(k).
The value after convergence is the posterior estimate of the state x(k|k).

A detailed description of the proposed MCUF algorithm is as follows:

1) Choose a proper kernel bandwidth σ and a small positive ε; Set an initial estimate x(0|0) and corresponding covariance matrix P(0|0); Let k = 1; 2) Use equations ( 9)∼( 14) to obtain prior estimate x(k|k -1) and covariance P(k|k -1), and calculate S p (k|k -1) by Cholesky decomposition ; 3) Use ( 14)∼( 17) to compute the prior measurement y(k) and use (13) ( 18) and ( 20) to acquire the measurement slope matrix H(k), and construct the statistical linear regression model ( 22); 4) Transform ( 22) into (24), and let t = 1 and x(k|k

) 0 = W T (k)W(k) -1 W T (k)D(k);
5) Use ( 36)∼(42) to compute x(k|k) t ;

x(k|k

) t = x(k|k -1) + K(k) (y(k) -y(k)) (36)
with

K(k) = P(k|k -1)H T (k) H(k) P(k|k -1)H T (k) + R(k) -1 , (37)
P(k|k -1) = S p (k|k -1) C -1 x (k)S T p (k|k -1), ( 38)

R(k) = S r (k) C -1 y (k)S T r (k). ( 39
) C x (k) = diag (G σ ( e 1 (k)) , ..., G σ ( e n (k))) (40) C y (k) = diag (G σ ( e n+1 (k)) , ..., G σ ( e n+m (k))) (41) e i (k) = d i (k) -w i (k) x(k|k) t-1 (42)
6) Compare the estimation at the current step and the estimation at the last step. If (43) holds, set x(k|k) = x(k|k) t and continue to step 7); Otherwise, t + 1 → t, and go back to step 5).

x(k|k) t -x(k|k) t-1

x(k|k) t-1 ≤ ε (43) 7) Update the posterior covariance matrix by ( 44), k + 1 → k and go back to step 2).

P(k|k) = I -K(k)H(k) P(k|k -1) I -K(k)H(k) T + K(k)R(k) K T (k) (44)
Remark 2: As one can see, (9)∼( 18) are the unscented transformation (UT). In the MCUF, we use a statistical linear regression model and the MCC to obtain the posterior estimates of the state and covariance. With the UT and statistical linear approximation, the proposed filter can achieve a more accurate solution than the first order linearization based filters. Moreover, the usage of MCC will improve the robustness of the filter against large outliers. Usually the convergence of the fixed-point iteration to the optimal solution is very fast (see Section 4). Thus, the computational complexity of MCUF is not high. The kernel bandwidth σ, which can be set manually or optimized by trial and error methods in practical applications, is a key parameter in MCUF. In general, a smaller kernel bandwidth makes the algorithm more robust (with respect to outliers), but a too small kernel bandwidth may lead to slow convergence or even divergence of the algorithm. From (B. [START_REF] Chen | Convergence of a fixed-point algorithm under maximum correntropy criterion[END_REF], we know that if the kernel bandwidth is larger than a certain value, the fixed-point equation ( 29) will surely converge to a unique fixed point. When σ → ∞, the MCUF will obtain the least-squares solution

W T (k)W(k) -1 W T (k)D(k).

Illustrative examples

In this section, we present two illustrative examples to demonstrate the performance of the proposed MCUF algorithm. Moreover, the performance is measured using the following benchmarks:

MSE 1 (k) = 1 M M m=1 (x(k) -x(k|k)) 2 , for k = 1 . . . K (45) MSE 2 (m) = 1 K K k=1 (x(k) -x(k|k)) 2 , for m = 1 . . . M (46) MSE = 1 M M m=1 MSE 2 (m) = 1 K K k=1 MSE 1 (k) ( 47
)
where K is the total time steps in every Monte Carlo run and M represents the total number of Monte Carlo runs.

Example 1

Consider the univariate nonstationary growth model (UNGM), which is often used as a benchmark example for nonlinear filtering. The state and measurement equations are given by 48)

x(k) = 0.5x(k -1) + 25 x(k -1) 1 + x(k -1) 2 + 8 cos (1.2(k -1)) + q(k -1), (
y(k) = x(k) 2 20 + r(k). (49)
First, we consider the case in which the noises are all Gaussian, that is,

q(k -1) ∼ N (0, 1) r(k) ∼ N (0, 1)
Table 1 lists the MSEs of x, defined in (47), and the average fixed-point iteration numbers. In the simulation, the parameters are set as K = 500, M = 100. Since all the noises are Gaussian, the UKF achieves the smallest MSE among all the filters. In this example, we should choose a larger kernel bandwidth in MCUF to have a good performance. One can also observe that the average iteration numbers of the MCUF are relatively small especially when the kernel bandwidth is large. Second, we consider the case in which the process noise is still Gaussian but the measurement noise is a heavy-tailed (impulsive) non-Gaussian noise, with a mixed-Gaussian distribution, that is, q(k -1) ∼ N (0, 1) r(k) ∼ 0.8N (0, 1) + 0.2N (0, 400) Table 2 illustrates the corresponding MSEs of x and average iteration numbers. As one can see, in impulsive noises, when kernel bandwidth is too small or too large, the performance of MCUF will be not good. However, with a proper kernel bandwidth (say σ = 2.0), the MCUF can outperform all the filters, achieving the smallest MSE. In addition, it is evident that the larger the kernel bandwidth, the faster the convergence speed. In general, the fixed-point algorithm in MCUF will converge to the optimal solution in only few iterations. We also investigate the influence of the threshold ε on the performance. The MSEs of x and average iteration numbers with different ε (The kernel bandwidth is set at σ = 2.0) are given in Table 3. Usually, a smaller ε results in a slightly lower MSE but a larger iteration number for convergence. Without mentioned otherwise, we choose ε = 10 -6 in this work. Further, we consider the situation where the process and measurement noises are all non-Gaussian with mixed-Gaussian distributions: q(k -1) ∼ 0.8N (0, 0.1) + 0.2N (0, 10) r(k) ∼ 0.8N (0, 1) + 0.2N (0, 400)

With the same parameters setting as before, the results are presented in Table 4. As expected, with a proper kernel bandwidth the MCUF can achieve the best performance.

Example 2

In this example, we consider a practical model [START_REF] Simon | Optimal state estimation: Kalman, H ∞ and nonlinear approaches[END_REF]. and the performance of EKF [START_REF] Anderson | Optimal filtering[END_REF], Huber-EKF (HEKF) (El-Hawary & Jing, 1995), UKF [START_REF] Julier | A new method for the nonlinear transformation of means and covariances in filters and estimators[END_REF] and HUKF (X. [START_REF] Wang | Huber-based unscented filtering and its application to vision-based relative navigation[END_REF] are also presented for comparison purpose. The goal is to estimate the position, velocity and ballistic coefficient of a vertically falling body at a very high altitude. The measurements are taken by a radar system each 0.1s. By rectangle integral in the discrete time with sampling period ∆T , we obtain the following model:

x 1 (k 1 ) = x 1 (k 1 -1) + ∆T x 2 (k 1 -1) + q 1 (k 1 -1) x 2 (k 1 ) = x 2 (k 1 -1) + ∆T ρ 0 exp (-x 1 (k 1 -1)/a) ×x 2 (k 1 -1) 2 x 3 (k 1 -1)/2 -∆T g + q 2 (k 1 -1) x 3 (k 1 ) = x 3 (k 1 -1) + q 3 (k 1 -1) (50) y(k) = b 2 + (x 1 (k) -H) 2 + r(k) (51)
where the sampling time is ∆T = 0.001s, that is k = 100k 1 , the constant ρ 0 is ρ 0 = 2, the constant a that relates the air density with altitude is a = 20000, the acceleration of gravity is g = 32.2f t/s 2 , the located altitude of radar is H = 100000f t, and the horizontal range between the body and the radar is b = 100000f t.

The state vector x(k) = x 1 (k) x 2 (k) x 3 (k) T contains the position, velocity and ballistic coefficient. Similar to [START_REF] Simon | Optimal state estimation: Kalman, H ∞ and nonlinear approaches[END_REF], we do not introduce any process noise in this model. The initial state is assumed to be x(0) = 300000 -20000 1/1000 T , and the initial estimate is x(0|0) = 300000 -20000 0.0009 T with covariance matrix P(0|0) = diag([1000000, 4000000, 1/1000000]). First, we assume that the measurement noise is Gaussian, that is,

q 1 (k 1 ) ∼ 0 q 2 (k 1 ) ∼ 0 q 3 (k 1 ) ∼ 0 r(k) ∼ N (0, 10000)
In the simulation, we consider the motion during the first 50s, and make 100 independent Monte Carlo runs, that is, K = 500, M = 100. Fig. 1 ∼ Fig. 3 show the MSE 1 (as defined in (45)) of x 1 , x 2 and x 3 for different filters in Gaussian noise. The corresponding MSEs and average fixed-point iteration numbers are summarized in Table 5 andTable 6. As one can see, in this case, since the noise is Gaussian, the UKF performs the best. However, the proposed MCUF with large kernel bandwidth can outperform the other two robust Kalman type filters, namely HEKF and HUKF. One can also observe that the average iteration numbers of the MCUF are very small especially when the kernel bandwidth is large. 1.0521

Second, we consider the case in which the measurement noise is a heavy-tailed (impulsive) non-Gaussian noise, with a mixed-Gaussian distribution, that is, q 1 (k 1 ) ∼ 0 q 2 (k 1 ) ∼ 0 q 3 (k 1 ) ∼ 0 r(k) ∼ 0.7N (0, 1000) + 0.3N (0, 100000) Fig. 4 ∼ Fig. 6 demonstrate the MSE 1 of x 1 , x 2 and x 3 for different filters in non-Gaussian noise, and Table 7 andTable 8 summarize the corresponding MSEs and average iteration numbers respectively. We can see clearly that the three robust Kalman type filters (HEKF, HUKF, MCUF) are superior to their non-robust counterparts (EKF, UKF). When the kernel bandwidth is very large, the MCUF achieves almost the same performance as that of UKF. In contrast, with a smaller kernel bandwidth, the MCUF can outperform the UKF significantly. Especially, when σ = 2.0, the MCUF exhibits the smallest MSE among all the algorithms. Again, the fixed-point algorithm in MCUF will converge to the optimal solution in very few iterations. 1.0917

Conclusion

In this work, we propose a novel nonlinear Kalman type filter, namely the maximum correntropy unscented filter (MCUF), by using the unscented transformation (UT) to get the prior estimates of the state and covariance matrix and applying a statistical linearization regression model based on the maximum correntropy criterion (MCC) to obtain the posterior estimates (solved by a fixedpoint iteration) of the state and covariance. Simulation results demonstrate that with a proper kernel bandwidth, the MCUF can achieve better performance than some existing algorithms including EKF, HEKF, UKF and HUKF particularly when the underlying system is disturbed by some impulsive noises.

D (k) = S -1 (k)

x (k|k -1) y (k) -h (k, x(k|k -1)) + H(k) x(k|k -1) = S -1 p (k|k -1) x (k|k -1) S -1 r (y(k) -h (k, x(k|k -1)) + H(k) x(k|k -1)) (A3) By (A1) and (A2), we have

W T (k) C (k) W (k) -1 = S -1 p T C x S -1 p + H T S r -1 T C y S -1 r H -1 (A4)
where we denote S p (k|k -1) by S p , S r (k) by S r , C x (k) by C x and C y (k) by C y for simplicity.

Using the matrix inversion lemma with the identification:

S -1 p T C x S -1 p → A, H T → B, H → C, S -1 r T C y S -1 r → D.
We arrive at

W T (k) C (k) W (k) -1 = S p C -1 x S T p -S p C -1 x S T p H T (S r C -1 y S T r + HS p C -1 x S T p H T ) -1 HS p C -1 x S T

Filter

Filter

Figure 1 .Filter

 1 Figure 1. MSE 1 of x 1 in Gaussian noise

Figure 4 .Filter

 4 Figure 4. MSE 1 of x 1 in non-Gaussian noise

  y S -1 r (y(k) -h (k, x(k|k -1)) + H(k) x(k|k -1)) (A6)Combining (32), (A5) and (A6), we have (33).

Table 1 .

 1 MSEs of x and Average Iteration Numbers in Gaussian Noises.

	Filter	MSE of x Average iteration number
	UKF	68.9766	-
	MCUF σ = 2.0, ε = 10 -6	108.9796	5.0624
	MCUF σ = 3.0, ε = 10 -6	94.5856	4.5431
	MCUF σ = 5.0, ε = 10 -6	83.7554	3.7323
	MCUF σ = 8.0, ε = 10 -6	86.1612	3.0431
	MCUF σ = 10, ε = 10 -6	85.4109	2.7919

Table 2 .

 2 MSEs of x and Average Iteration Numbers in Gaussian ProcessNoise and Non-Gaussian Measurement Noise.

Table 3 .

 3 MSEs of x and Average Iteration Numbers with Different ε.

	Filter	MSE of x Average iteration number
	MCUF σ = 2.0, ε = 10 -1	69.7465	1.1142
	MCUF σ = 2.0, ε = 10 -2	69.5853	1.3777
	MCUF σ = 2.0, ε = 10 -4	69.0197	2.1753
	MCUF σ = 2.0, ε = 10 -6	68.8073	3.0282
	MCUF σ = 2.0, ε = 10 -8	68.8148	3.8770

Table 4 .

 4 MSEs of x and Average Iteration Numbers in Non-Gaussian Process Noise and Measurement Noise.

Funding

This work was supported by 973 Program (No. 2015CB351703) and the National Natural Science Foundation of China (No. 61372152).

Appendix A. Derivation of (33)