GLOBAL SENSITIVITY ANALYSIS OF COMPUTER MODELS WITH FUNCTIONAL INPUTS

B. Iooss1*, M. Ribatet2
1CEA Cadarache, France, 2CEMAGREF Lyon, France
bertrand.iooss@cea.fr

Introduction

This work addresses the problem of global sensitivity analyses of numerical models when some input parameters e are not scalar but functional: e(x) is a one or multi-dimensional stochastic function where x can be spatial coordinates, time scale or any other physical parameter. This work focuses with models which depend on scalar parameter vector p and use some stochastic processes simulations or random fields e(x) during their running. The output value Y of the computer code depends on the realizations of these random functions. For example, the oil reservoirs production evaluations depend on geostatistical simulations of the heterogeneous geological media (Zabalza-Mezghani et al. [7], Ruffo et al. [4]), the environmental assessment on soil models use spatially distributed maps affected by random errors (Tarantola et al. [6]), the thermo-mechanical behaviour of the fuel rod under nuclear irradiation depends on temporal parameters modelized as stochastic processes (Iooss and Ribatet [1]), … These models are typically non linear with strong interactions between input parameters. Therefore, we concentrate our methodology on the estimation of the variance based sensitivity indices; that is, the so-called Sobol indices (Sobol [5]).
To deal with this situation, a first natural approach consists in the discretization of the input functional parameter e(x) or its decomposition in an appropriate basis of orthogonal functions. Then, for all the new scalar parameters which represent e(x), sensitivity indices are computed. However, in the case of complex functional parameters, this approach seems to be rapidly impracticable as these parameters cannot be represented by a small number of scalar parameters (Tarantola et al. [6]). Moreover, when dealing with non physical parameters, sensitivity indices interpretation would be quite difficult. Indeed, most often, physicists would prefer to obtain one global sensitivity index concerning e(x). 
In this work, four strategies to resolve this problem are presented. The advantages and drawbacks of each approach are discussed and some applications will be shown during the conference presentation.

1.	Sobol indices by Monte Carlo calculations.

To resolve the problem of correlated input parameters in the Sobol indices calculations, Jacques et al. [2] have proposed to use the multi-dimensional sensitivity indices (Sobol [5]): each group of correlated parameters is considered as a multi-dimensional parameter. The simple Monte-Carlo algorithms (Sobol [5]) can then be used to calculate the different Sobol indices (first order, second order, …, total). In our case, this solution is fully applicable by considering the input functional parameter e(x) as an unique input multi-dimensional parameter. However, it is supposed that e(x) is not correlated with any other input parameters. Moreover, the calculation of Sobol indices by the simple Monte-Carlo method requires more than thousand model evaluations for one input parameter. In complex industrial applications, it is often unrealisable.
 
2.	Using an intermediate “trigger” parameter governing the random function simulation.

Tarantola et al. [6] have proposed another solution by introducing a scalar input parameter x ~ U[0,1] governing the simulation of the random function.  For each sample simulation, if x < 0.5, the functional parameter e(x) is fixed to a nominal value e0(x) (for example the mean <e(x)>); if x > 0.5, the functional parameter e(x) is simulated. Therefore, the sensitivity index of x quantifies the influence of the random function on the model output variable.  Contrary to the previous case, various sensitivity analysis methods can be applied (regression, FAST, Sobol, …). However, as the previous method, it also requires the use of the computer model to perform the sensitivity analysis. Moreover, x reflects only the presence or the absence of the stochastic errors on e0(x), and the term Var[E(Y|x)] does not quantify properly the contribution of e(x) to Var(Y).

3.	Using a metamodel.

To perform a variance-based sensitivity analysis on long running time model, some work propose to replace the computer code by an approximate mathematical function, called a metamodel (Ruffo et al. [4]). For metamodels with sufficient prediction capabilities, the bias due to the use of the metamodel instead of the true model is negligible. Several choices of metamodel can be found in the literature: polynomials, splines, Gaussian processes, neural networks, … Therefore, a first solution would be to fit a metamodel including the multi-dimensional scalar parameters representing e(x) (i.e. its discretization or its decomposition in an appropriate basis). However, this seems to be impracticable due to the potential large number of scalar parameters. A second solution consists in replacing the continuous parameter x by a discrete parameter x’ =1x > 0.5 , which can correspond to the scenario parameter of Ruffo et al. [4] (where the number of geostatistical realizations is finite and fixed, and where each different value of the discrete parameter corresponds to a different realization). However, the situation differs in our context; as the previous method, the term Var[E(Y|x’)] does not quantify properly the contribution of e(x) to Var(Y). 
The last solution considers e(x) as an uncontrollable parameter and a metamodel is fitted in function of the other scalar parameters p: Ym = E(Y|p). Therefore, using the relation Var(Y) = Var[E(Y|p)] + E[Var(Y|p)], it can be seen that the total sensitivity index of e(x) corresponds to the expectation of the unexplained part of the metamodel E[Var(Y|p)]. However, classical parametric metamodels are not adapted to correctly estimate E(Y|p) because of the heteroscedasticity of this situation (some examples are shown in Iooss and Ribatet [1]). This can cause some misspecifications of the sensitivity indices.

4.	The joint modeling approach

Zabalza-Mezghani et al. [7] have proposed a solution by using a theory developed for experimental data: the simultaneous fitting of the mean and the dispersion by two interlinked generalized linear models (McCullagh and nelder [3]).  This approach, called the joint modeling, has been recently extended by Iooss and Ribatet [1] to non parametric models. This generalization allows more complexity and flexibility while fitting the data. They propose to use generalized additive models (GAMs) based on penalized smoothing splines. The GAM allows model and variable selections via quasi-likelihood function, statistical tests on coefficients, graphical display, … Its restriction, compared to other complex metamodels stands in the additive effects hypothesis.
Therefore, two metamodels are obtained: one for the mean component Ym = E(Y|p); and the other one for the dispersion component Yd = Var(Y|p). The sensitivity indices of p are computed using Ym with the standard procedure, while the total sensitivity indice of e(x) is computed from E(Yd). Using the explicit formula on Yd and the associated regression diagnostics, qualitative sensitivity indices for the interactions between e(x) and the scalar parameters of p can be deduced.

Acknowledgement:

	We thank S. Cambier for useful remarks about the Sobol indices estimation on this kind of models.

References

[1] Iooss B, Ribatet M: Global sensitivity analysis of stochastic computer models with generalized additive models, Technometrics, submitted (2006)

[2] Jacques J, Lavergne C, Devictor N: Sensitivity analysis in presence of model uncertainty and correlated inputs, Reliability Engineering and System Safety 91, 1126-1134 (2006)

[3] McCullagh P, Nelder JA: Generalized linear models, Chapman & Hall, 1989.

[4] Ruffo P, Bazzana L, Consonni A, Corradi A, Saltelli A, Tarantola S: Hydrocarbon exploration risk evaluation through uncertainty and sensitivity analyses techniques, Reliability Engineering and System Safety 91, 1155-1162 (2006)

[5] Sobol IM: Global sensitivity indices for non linear mathematical models and their Monte Carlo estimates, Mathematics and Computers in Simulation 55, 271-280 (2001)

[6] Tarantola S, Giglioli N, Jesinghaus N, Saltelli A: Can global sensitivity analysis steer the implementation of models for environmental assesments and decision-making?, Stochastic Environmental Research and Risk Assesment 16, 63-76(2002)

[7] Zabalza-Mezghani I, Manceau E, Feraille M, Jourdan A: Uncertainty management: From geological scenarios to production scheme optimization, Journal of petroleum Science and Engineering 44, 11-25 (2004)