ANALYTICAL CALCULATIONS OF SOBOL INDICES FOR THE GAUSSIAN PROCESS METAMODEL

A. Marrel1*, B. Iooss1, O. Roustant2
1CEA Cadarache, France; 2Ecole des Mines de Saint Etienne, France
amandine.marrel@cea.fr

Introduction

This paper deals with the analytical formulations for the Gaussian process model of the first order Sobol indices (1) where X=(X1,…, Xd) and Y is the model input and output variables. The Gaussian process model (Sacks et al. [3]), also called kriging, is often used as a metamodel to approximate complex and expensive computer codes. Because of their lower computational cost, metamodels are useful in uncertainty, sensitivity and optimisation studies, which require a large number of function evaluations. The estimation of one Sobol index is usually done by sampling adequately the input parameters space via Monte Carlo or other algorithms, which can require thousands model evaluations (Saltelli et al. [4]). In these algorithms, the full model (computer code) can be directly replaced by a metamodel. Some recent works have addressed the possibility to take advantage of the known analytical formula of some metamodels to derive analytical expressions of Sobol indices: Chen et al. [1] which use tensor-product formulation and Oakley and O’Hagan [2] which consider the Bayesian formalism of Gaussian processes. Our work proposes to develop the analytical formulation of Sobol indices for the Gaussian process model viewed as a global stochastic model and not as a simple predictor (as in Chen et al. [1]). On a simple example, comparisons are made between the Sobol indices based on the predictor and those based on the global model, in terms of the metamodel predictivity.

The Gaussian process model (Gp model)

Gp modeling considers the deterministic computer code output ycode(x) as the realization of a random function Y(x,ω) defined on Ω x Rd, including a regression part and a Gaussian stochastic process Z(x,ω):
      (2)
with b a regression parameter vector and F(x)=[f0(x),…, fk(x)] a vector of elementary functions. We take Z(x,ω) centred, second order stationary with variance s2 and a generalized exponential correlation function:  with qi>0 and 0<pi2. Gp model parameters (b, s2, q, p) are estimated by maximum likelihood on a learning sample (Xs, Ys) made by n simulations of computer code: Xs=(x(1),…,x(n)), Ys=ycode(Xs). Conditioning the joint probability distribution by the observed data (Xs, Ys), the conditional Gaussian process YGp/(Xs, Ys) is defined on Ω x Rd and characterized by:
                  (3) 
where  , Fs=F(Xs) and is the covariance matrix of Xs. The expression of conditional Gp model mean is used as a predictor and often denoted . The variance expression corresponds to the mean square error of the predictor and gives a local indicator of prediction accuracy. These analytical expressions are important advantages in uncertainty propagation, optimization and sensitivity analysis studies, contrary to other metamodels which give only the predictor. To validate the predictor, we use a test sample, made by new computer code simulations, to estimate the predictive residual sum of squares (PRESS) and  its associated predictivity coefficient Q2 (corresponding to R² for prediction points). 

Sensitivity analysis

From the expression of conditional Gp model, two approaches to compute Sobol indices are possible: 
-	only the predictor expression is used, in other words the mean of Gp process;
-	the whole global conditional Gp process YGp/(Xs, Ys) is used.
The second approach takes into account not only the mean of conditional Gp model but also its covariance structure. In this case, following the Sobol approach (1), we define a new sensitivity measure  which is a random variable. Its expectation can be considered as a sensitivity indice (denoted) and its variance (denoted ²) as an indicator of sensitivity indice accuracy. The two approaches can be defined as follows:
         (4) 
For the two approaches, the analytical formulas above lead to numerical integrals, and more precisely to one-dimensional integrals in the case of independent inputs, which is another advantage of Gp metamodel.

To compare and study the two indices, we consider a simple example based on the Ishigami function where Xi ~U[- π ; π ]. The theoretical values of Sobol indices are known: S1 = 0.42, S2 = 0.58 and S3 = 0. To study the convergence of and relatively to the accuracy of Gp modeling represented by the coefficient Q2, we simulate different learning samples with varying size n.  For each size n, a Latin Hypercube sample is simulated, the Ishigami function is evaluated on the n data points, the conditional Gp model is built on this learning sample and the predictivity coefficient Q2 is estimated on a new test sample of size 10000. The Sobol indices are computed for the two approaches (4). To observe the results, empirical mean and standard deviation of each indice are computed for same values of Q2 (or same classes of Q2). Figure (1) illustrates the convergence of and  and figure (2) shows the convergence of the ratio of S1 and S2. Sobol indices computed with the whole model is more robust and less dispersed, particularly for low accuracy values of Gp model (Q² < 0.8). The convergence to the true value is faster with the whole Gp model. However, for high values of accuracy of Gp model (Q² > 0.9), the two approaches give the same values and the first one (with only the predictor) remains easier to compute. In all cases, the second approach taking into account the covariance structure is more robust and allows to estimate the uncertainty of the sensitivity coefficient with . It justifies the use of conditional Gp model as metamodel.                  

Figure 2: Convergence of the Sobol indice ratio S1/S2

Figure 1: Convergence of Sobol indice S3

Acknowledgement:

This work was supported by CEA/Nuclear Development and Innovation Division/Risk Control Project.

References
[1] Chen W, Jin R, Sudjianto A: Analytical metamodel-based global sensitivity analysis and uncertainty propagation for robust design. Journal of Mechanical Design 127, 875-886 (2005).
[2] Oakley JE, O’Hagan A: Probabilistic sensitivity analysis of complex models: A Bayesian approach, Journal of the Royal Statistical Society, Series B 66, 751-769 (2004)
[3] Sacks J, Welch WJ, Mitchell TJ, Wynn HP: Design and analysis of computer experiments, Statistical Science 4, 409-435 (1989)
[4] Saltelli A, Chan K, Scott EM, editors: Sensitivity Analysis, Wiley, 2000.