Ivano Azzini

Thierry A Mara
email: thierry.mara@ec.europa.eu

Rossana Rosati

Comparison of two sets of Monte Carlo estimators of Sobol' indices

Keywords: global sensitivity analysis, variance-based sensitivity indices, first-order Sobol' index, total-order Sobol' index, asymptotic normality, radiative forcing model

This study compares the performances of two sampling-based strategies for the simultaneous estimation of the first-and total-order variance-based sensitivity indices (a.k.a. Sobol' indices). The first strategy corresponds to the current approach employed by practitioners and recommended in the literature. The second one was only recently introduced by the first and last authors of the present article. Both strategies rely on different estimators of first-and total-order Sobol' indices. The asymptotic normal variances of the two sets of estimators are established and their accuracies are compared theoretically and numerically. The results show that the new strategy outperforms the current one. The global sensitivity analysis of the radiative forcing model of sulfur aerosols is performed with the new strategy. The results confirm that in this model interactions are important and only one input variable is irrelevant.

Background

Uncertainty and sensitivity analyses are essential ingredients of modelling [START_REF] Saltelli | Sensitivity analysis in practice[END_REF]. They allow to point out the key uncertain assumptions (input factors that can be random variables or random fields) responsible for the uncertainty into the model outcome of interest. This is particularly relevant when models are used for decision-making.

Assessing model output uncertainty requires several runs of the model.

Monte Carlo simulations allow to carry out this task by sampling the input factors from their joint probability distribution and propagating the sample through the model response of interest (i.e. running the model). Sensitivity analysis (SA) can then be undertaken to identify, for instance, the input factors mostly responsible for the uncertainty in the model response. Depending on the method used, SA can be conducted directly from the Monte Carlo sample at hand (i.e. the one generated to assess model output uncertainty) or can require extra Monte Carlo simulations by following an appropriate sampling design.

The method to be used depends on the sensitivity indices (also called importance measures) that the analyst wants to compute. As recommended in [START_REF] Saltelli | Sensitivity analysis in practice[END_REF] (see also [START_REF] Saltelli | On the relative importance of input factors in mathematical models: Safety assessment for nuclear waste disposal[END_REF]), the sensitivity indices to assess should be related to the question that SA is called to answer to. The authors enumerate several questions (called SA settings) that can be addressed with the so-called variance-based sensitivity indices. In the sequel, we focus on the estimation of variance-based sensitivity indices, also called Sobol' indices ( [START_REF] Sobol | Sensitivity estimates for nonlinear mathematical models[END_REF]).

As eluded previously, a Monte Carlo sample is required to carry out uncertainty analysis, that is, assessing the predictive uncertainty of the model output of interest. We assume that there is only one scalar output denoted y = f (x). The input factors are represented by a random vector of scalar variables x = (x 1 , . . . , x D ) possibly grouped into two complementary vectors (u, v). They are assumed independent of each other (for the case of dependent inputs, see for instance [START_REF] Li | Global sensitivity analysis for systems with independent and/or correlated inputs[END_REF][START_REF] Kucherenko | Estimation of global sensitivity indices for models with dependent variables[END_REF][START_REF] Mara | Variance-based sensitivity indices for models with dependent inputs[END_REF][START_REF] Mara | Non-parametric methods for global sensitivity analysis of model output with dependent inputs[END_REF][START_REF] Li | Relationship between sensitivity indices defined by variance-and covariance-based methods[END_REF][START_REF] Tarantola | Variance-based sensitivity indices of computer models with dependent inputs: the Fourier amplitude sensitivity test[END_REF]).

There exist several Sobol' indices called, first-order, (closed) second-order, and so forth. A closed second-order Sobol' index (and more generally a closed d-th order Sobol' index) can be defined as the first-order Sobol' index of a group of two (resp. d) inputs. In the sequel, we will use the term first-and total-order Sobol' indices whether they refer to an individual variable (say x i ) or a group of variables (e.g. u).

First-and total-order Sobol' indices are respectively defined as follows, 1) is the first-order Sobol' index of the group of inputs u while Eq. ( 2) is the total-order Sobol' index of u. The higher the Sobol' indices, the more the group of inputs u is important for the model response.

S u = V [E [y|u]] V [y] (1)
ST u = E [V [y|v]] V [y] (2)
ST u ≥ S u . Eq. (
The difference between S u and ST u is that the latter not only accounts for the amount of variance of y explained by the input variables within u (like S u does) but it also contains cooperative contributions due to the interactions between the variables in u with those in v. Therefore, a noticeable result is

that S u + ST v = 1.
Let d ∈ [1, D) be the number of elements in u. S u represents in percentage, the expected reduction in V [y] if the variables in u where fixed to their true value. That is why the first-order sensitivity indices of individual inputs (i.e d = 1) are to be estimated if the goal of the SA is to identify the input variable that would induce the largest reduction in variance if its value was known accurately. This SA setting is called factors prioritization. Instead, if the goal is to identify the irrelevant inputs (called screening analysis or factors fixing setting) then the individual total-order Sobol' indices are to be estimated. Indeed, we note that if ST u = 0, the variables in u do not contribute at all to the variance of y. More SA settings are discussed in [START_REF] Saltelli | On the relative importance of input factors in mathematical models: Safety assessment for nuclear waste disposal[END_REF].

There are typically two direct methods to estimate the first-and totalorder Sobol' indices. The first one uses Monte Carlo methods (e.g. [START_REF] Sobol | Sensitivity estimates for nonlinear mathematical models[END_REF][START_REF] Saltelli | Making best use of model evaluations to compute sensitivity indices[END_REF]).

The second one casts the total variance onto orthogonal functions like the Fourier expansion (a.k.a. the Fourier amplitude sensitivity test, [START_REF] Cukier | Nonlinear sensitivity analysis of multiparameter model systems[END_REF][START_REF] Saltelli | A quantitative model independent method for global sensitivity analysis of model output[END_REF][START_REF] Mara | Extension of the rbd-fast method to the computation of global sensitivity indices[END_REF]) or the polynomial chaos expansion [START_REF] Sudret | Global sensitivity analysis using polynomial chaos expansions[END_REF][START_REF] Blatman | Adaptive sparse polynomial chaos expansion based on least angle regression[END_REF][START_REF] Shao | Bayesian sparse polynomial chaos expansion for global sensitivity analysis[END_REF]. Indirect methods employ surrogate models first (also called metamodels) to mimic the input-output relationship, and then often use one of the aforementioned direct methods to compute the Sobol' indices (e.g. [START_REF] Marseguerra | Variance decomposition-based sensitivity analysis via neural networks[END_REF]). Estimating the Sobol' indices with Monte Carlo estimators is rather computationally expensive, but it does not require any assumption except that the variance of f (x) be numerically tractable. In the present work, we study the performances of two Monte Carlo estimators of Eq. ( 1) and Eq. ( 2) respectively that rely on two different sampling designs.

The paper is organised as follows: in Section 2 we introduce the two sampling strategies as well as their associated Monte Carlo estimators to compute both the first-and total-order Sobol' indices. Their asymptotic normal variances, derived in the appendices, are also compared to each other.

In Section 3, the performances of the two sets of estimators are compared through numerical exercises on well-known benchmark functions. The new set of estimators is applied to the radiative forcing model of sulfur aerosols in Section 4. Finally, the key results of our work are summarized in Section 5.

Monte Carlo estimators

Integral approximation

When Ilya M. Sobol' introduced the variance-based sensitivity indices in [START_REF] Sobol | Sensitivity estimates for nonlinear mathematical models[END_REF], he also proposed their Monte Carlo (MC) estimators. Indeed, the Sobol' indices defined in Eqs. (1-2) are nothing but ratio of integrals. Approximating these integrals numerically provides estimates of the Sobol' indices. Monte Carlo (MC) estimators rely on the fact that multidimensional integrals can be approximated via MC samples as follows,

R D f q (x)p x (x)dx ≈ 1 N N n=1 f q (x n,1 , • • • , x n,D ) (3)
where x ∼ p x , meaning that p x is the joint probability density of x and

x n = (x n,1 , • • • , x n,D
) is the n-th (out of N ) MC draw of the input factors sampled w.r.t. p x . In the sequel, we assume that x is a vector of independent input variables, that is, p x = D i=1 p x i (x i ), where p x i (x i ) is the marginal distribution of x i . Let (y A , y B , y Au , y Bu ) be four distinct model output samples whose n-th element for each of them is respectively defined as follows,

y A n = f (u A n , v A n ) = f (x A n ) y B n = f (u B n , v B n ) = f (x B n ) y Au n = f (u A n , v B n ) = f (x Au n ) y Bu n = f (u B n , v A n ) = f (x Bu n )
where

Current estimators

The most popular sampling design to compute simultaneously the firstand total-order sensitivity indices as recommended by Saltelli et al. [START_REF] Saltelli | Update 1 of: Sensitivity analysis for chemical models[END_REF] requires three samples, namely (y A , y B , y Au ) (or equivalently (y A , y B , y Bu )),

to compute the sensitivity indices of u. Their estimators are respectively defined as follows,

S SS u = 1 N N n=1 y A n y Au n -y B n 1 2N N n=1 (y A n -y B n ) 2 (4)
ST SJ u = 1 2N N n=1 y Au n -y B n 2 1 2N N n=1 (y A n -y B n ) 2 .
(

) 5
Note that we do not simplify these equations (e.g., the 2N at the numerator and denominator cancel each other) for the purpose of the discussion that just follows, but latter on, we will. The superscript SS stands for Sobol-Saltelli as the former derived the integral formulation of the numerator in [START_REF] Sobol | Estimating the approximation error when fixing unessential factors in global sensitivity analysis[END_REF] while

Saltelli proposed an estimator similar to the numerator of Eq. ( 4) in [START_REF] Saltelli | Making best use of model evaluations to compute sensitivity indices[END_REF].

The superscript SJ often refers to Sobol-Jansen although one can date back the numerator of Eq. ( 5) to Šaltenis and Dzemyda [START_REF] Šaltenis | Structure analysis of extremal[END_REF] and Jansen et al. [START_REF] Jansen | Monte Carlo estimation of uncertainty contributions from several independent multivariate sources[END_REF] (see [START_REF] Saltelli | Sensitivity analysis[END_REF] page 177). Therefore, SJ can also be read Šaltenis-Jansen.

The denominators of the previous formulas are identical but they differ from the one proposed in [START_REF] Saltelli | Making best use of model evaluations to compute sensitivity indices[END_REF][START_REF] Saltelli | Variance based sensitivity analysis of model output. design and estimator for the total sensitivity index[END_REF]. As defined, the denominator of Eqs.

(4-5) is an MC estimator of V [y]. We find it convenient to formulate the denominator in this way because it highlights the symmetry between (y A , y B ) in the denominator. Eq. ( 4) is known to provide an accurate estimate of small first-order sensitivity indices [START_REF] Sobol | Estimating the approximation error when fixing unessential factors in global sensitivity analysis[END_REF] as Eq. ( 5) does for the total-order sensitivity indices.

Importantly, although in theory ST u ≥ S u , the previous estimators do not satisfy this criterion. Indeed, by noticing that

V y = 1 2N N n=1 y A n -y B n 2
at the denominator of Eqs. (4-5) is a positive scalar, we find that,

V y ST SJ u -S SS u = 1 N N n=1 y Au n -y B n 2 -2y A n y Au n -y B n (6)
which, because -2y A n y Au n -y B n can be either positive or negative, does not ensure that ST SJ u ≥ S SS u . These observations advocate for a more symmetrical and coherent estimator for the first-order sensitivity index. This is the subject of the next subsection.

New estimators

As previously mentioned, the denominator of Eq. ( 5) converges towards

V [y], that is, lim N →∞ 1 2N N n=1 y A n -y B n 2 = lim N →∞ 1 2N N n=1 y Au n -y Bu n 2 = V [y]
while the numerator is such that,

lim N →∞ 1 2N N n=1 y B n -y Au n 2 = lim N →∞ 1 2N N n=1 y A n -y Bu n 2 = E [V [y|v]] .
Hence, the following symmetrical estimator for the total-order sensitivity index can be derived,

ST IA u = N n=1 y B n -y Au n 2 + y A n -y Bu n 2 N n=1 (y A n -y B n ) 2 + (y Au n -y Bu n ) 2 . ( 7
)
This is because, as already mentioned, x A n and x B n are two independent input vectors identically distributed, as well as x Au n and x Bu n . Notice the perfect symmetry of the formula which remains unchanged by switching the superscripts B and A. Incidentally, the superscript IA stands indifferently for Innovative Algorithm and Ivano Azzini, the first author of this article who initiated the work on these estimators [START_REF] Azzini | Sobol' main effect index: An innovative algorithm (IA) using dynamic adaptive variances[END_REF].

Interchanging (y Au , y Bu ) in Eq. ( 7) only changes the numerator and provides the estimator for ST IA v . The law of total variance implies that

S IA u + ST IA v = 1.
Therefore, the first-order sensitivity index S u is estimated as follows,

S IA u = 1 -ST IA v = 1 - N n=1 y B n -y Bu n 2 + y A n -y Au n 2 N n=1 (y A n -y B n ) 2 + (y Au n -y Bu n ) 2
which after some developments yields,

S IA u = 2 N n=1 y Au n -y B n y A n -y Bu n N n=1 (y A n -y B n ) 2 + (y Au n -y Bu n ) 2 . ( 8
)
Besides,

ST IA u -S IA u = N n=1 y B n -y Au n 2 + y A n -y Bu n 2 -2 y Au n -y B n y A n -y Bu n N n=1 (y A n -y B n ) 2 + (y Au n -y Bu n ) 2 which yields, ST IA u -S IA u = N n=1 y B n -y Au n + y A n -y Bu n 2 N n=1 (y A n -y B n ) 2 + (y Au n -y Bu n ) 2 ≥ 0 (9)
and proves that ST IA u ≥ S IA u for any N . Eq. ( 9) also shows that ST IA u = S IA u if and only if f (x) is additive with respect to u. In effect, an additive function with respect to u and v reads in this case,

y = f (u, v) = f 0 + f u (u) + f v (v)
and it is straightforward to prove that the numerator of Eq. ( 9) equals zero for any sample size N .

It turns out that the numerator of Eq. ( 8) is very similar to the one proposed by Owen in [START_REF] Owen | Better estimation of small Sobol' indices[END_REF]. Apart from the factor 2 due to the denominator, the difference is the use by the author of

y Cu n = f (u C n , v A n ) instead of y Au n
or y Bu n in Eq. [START_REF] Jansen | Monte Carlo estimation of uncertainty contributions from several independent multivariate sources[END_REF]. By doing so, the symmetry of the estimator is lost.

Interestingly, Lamboni [START_REF] Lamboni | Uncertainty quantification: A minimum variance unbiased (joint) estimator of the non-normalized Sobol' indices[END_REF] also derived unbiased estimators with minimum variance of the non-normalized Sobol' indices by leaning on the theory of U -statistics (see [START_REF] Hoeffding | A class of statistics with asymptotically distribution[END_REF][START_REF] Gamboa | Sensitivity analysis in general metric spaces[END_REF]). His construction led to estimators exactly equal to the numerators of Eqs. [START_REF] Janon | Asymptotic normality and efficiency of two sobol index estimators[END_REF][START_REF] Jansen | Monte Carlo estimation of uncertainty contributions from several independent multivariate sources[END_REF]. However, neither Lamboni in [START_REF] Lamboni | Uncertainty quantification: A minimum variance unbiased (joint) estimator of the non-normalized Sobol' indices[END_REF] nor Owen in [START_REF] Owen | Better estimation of small Sobol' indices[END_REF] paid attention to the estimation of the total variance. The proof that ST IA u ≥ S IA u does not depend on the choice of the total variance estimator (i.e. the denominator of Eq. ( 9)). Therefore, the estimators proposed by Lamboni in [START_REF] Lamboni | Multivariate sensitivity analysis: Minimum variance unbiased estimators of the first-order and total-effect covariance matrices[END_REF][START_REF] Lamboni | Uncertainty quantification: A minimum variance unbiased (joint) estimator of the non-normalized Sobol' indices[END_REF] also satisfy this property but they do not form a set of complementary formulas contrarily to the new estimators defined by Eqs.

(7-8) (also called the IA estimators in the sequel). Indeed, thanks to the special choice of the total variance estimator (i.e. the denominator of Eqs.

(7-8)), the IA estimators also satisfy the following equation,

S IA u + ST IA v = 1.
To our best knowledge, there are no other estimators of first-and total-order sensitivity indices that comply with both the equation S u + ST v = 1 and the inequation S u ≤ ST u .

Variances of the estimators

The performance of an estimator is characterized by its bias and its variance. For a given sample size N , the Sobol' indices can be computed several times by re-sampling the MC draws, thus, providing different estimates. Unbiased estimators provide replicates that on average yield the true values of the Sobol' indices. Estimators with small variances provide estimates that remain close to the true values of the Sobol' indices. In the sequel, the focus is on the variances of the estimators.

In the Appendices A and B, we establish the variances of the estimators discussed in the present paper under the asymptotic normality assumption [START_REF] Van Der Vaart | Asymptotic statistics[END_REF][START_REF] Janon | Asymptotic normality and efficiency of two sobol index estimators[END_REF]. They respectively read as follows,

σ 2 SS = V 2y A y Au -y B -S u y A -y B 2 4V [y] 2 (10) τ 2 SJ = V y Au -y B 2 -ST u y A -y B 2 4V [y] 2 (11)
and,

σ 2 IA = V 2 y A -y Bu y Au -y B -S u y A -y B 2 + y Au -y Bu 2 16V [y] 2 (12) τ 2 IA = V y A -y Bu 2 + y B -y Au 2 -ST u y A -y B 2 + y Au -y Bu 2 16V [y] 2 .(13)
In Eqs. (10-13), y A , y B , y Au , and y Bu are random variables. In practice, to compute the variances of the estimators, they are replaced by their samples y A , y B , y Au , and y Bu , (see Section 3).

It can be qualitatively speculated that the Sobol-Jansen estimator is more accurate than the one of Sobol-Saltelli. Indeed, we have (according to [START_REF] Sobol | Sensitivity estimates for nonlinear mathematical models[END_REF]),

y = f (u, v) = f 0 + f u (u) + f v (v) + f u,v (u, v).
This implies that,

y Au -y B = -f u (u B ) + f u (u A ) -f u,v (u B , v B ) + f u,v (u A , v B ) y A -y Bu = -f u (u B ) + f u (u A ) -f u,v (u B , v A ) + f u,v (u A , v A ).
Therefore, the variance of y Au -y B 2 is expected to be smaller than the one of 2y A y Au -y B because the former does not contain neither f 0 , nor f v contrarily to the latter with y A . What is worse, we guess that Eq. ( 4) may perform very poorly for high values of f 0 . For the same reason, the variance of y Au -y B y A -y Bu is expected to be smaller than the one of 2y A y Au -y B . Therefore, the IA estimator of the first-order Sobol' index should also perform better than Sobol-Saltelli, especially when f 0 is high

compared to V [y].
It is less obvious to infer whether τ 2 SJ is higher or lower than τ 2 IA . Therefore, this is investigated through numerical simulations in the next section.

Numerical examples

It is worth noticing that the current estimators Eqs. for the latter. This means that when we write that a sample of size N is used, this refers to the actual size of the samples for the new estimators while the sample size is 2N for the current estimators. In the following exercises, we exclusively use the latin hypercube sampler (lhs) because it allows for the replication of the Sobol' indices estimate and furthermore it performs better than random sampling. For the interested readers, more intensive numerical exercises are undertaken in [START_REF] Azzini | Sobol' main effect index: An innovative algorithm (IA) using dynamic adaptive variances[END_REF] with different sampling techniques and different estimators of first-order Sobol' index.

The Ishigami function

Let us consider the following three-dimensional function,

f (x 1 , x 2 , x 3 ) = f 0 + sin x 1 + 7 sin 2 x 2 + 0.1x 4 3 sin x 1 ( 14
)
where the input variables are independently and uniformly distributed over (-π, π) 3 . As compared to the original Ishigami function, we introduce a constant parameter f 0 which has no impact on the variance of the function. This simple function for which the exact Sobol' indices are known has the following features: x 1 and x 3 interact strongly while x 2 is additively influential, that is, S 2 = ST 2 ≃ 0.44. This allows to check whether, as previously guessed, we find S IA 2 = ST IA 2 . In this exercise, we numerically compare the performances of Eqs. (4-5) with Eqs. [START_REF] Janon | Asymptotic normality and efficiency of two sobol index estimators[END_REF][START_REF] Jansen | Monte Carlo estimation of uncertainty contributions from several independent multivariate sources[END_REF]. For this purpose, we set N = 64 (which means 128 for the current estimators) and we replicate one hundred estimates of the first-and total-order Sobol' indices with the estimators discussed in this paper.

Each replicate is obtained as follows:

1. Generate x A and x B two independent lhs samples

Case 1:

f 0 = 0
We first set f 0 = 0. The results are depicted in Fig. 1 which clearly shows that, as far as the first-order Sobol' indices are concerned, the new estimator Eq. ( 8) provides more robust estimates than Eq. ( 4); thus confirming our comments in § 2.4. Notably, the estimated first-order Sobol' indices of x 3 can be smaller than zero which is not consistent with the theory (Sobol' indices shall be within [0,1]). This is due to its interaction with x 1 . The new totalorder estimator, that is Eq. ( 7), has slightly lower variances for ST 1 and ST 2 than Eq. ( 5) and conversely for ST 3 . Fig. 2 ). Actually for some replicates, we even find S SS 2 > ST SJ 2 which is not consistent at all with the definition of first-and total-order Sobol' indices. We stress that, when x i has only an additive effect on the response, S IA i = ST IA i is independent of the sample size N . This information can be obtained at any sample size (even for N ∼ 1). This is also illustrated in Fig. 2. The red crosses were obtained with N = 1, 2, . . . , 10 without replicate. At these sample sizes, the pairs ( S IA 2 , ST IA 2 ) are also located along the diagonal S IA 2 = ST IA 2 .

Case 2: f 0 = 100

This case illustrates the sensitivity of the current first-order estimator to model responses with high expected value as compared with the total variance. We set f 0 = 100, which yields,

f (x 1 , x 2 , x 3 ) = 100 + sin x 1 + 7 sin 2 x 2 + 0.1x 4 3 sin x 1 ( 15
)
keeping in mind that the Ishigami function has a total variance approximately equal to V y = 13.84. One hundred lhs replicates of size N = 64 are employed.

The results are displayed in Fig. 3. They show that while the shift in the Ishigami function has no impact on the estimators of the total-order estimators and on the new first-order estimator (namely, Eq. ( 8)), it significantly deteriorates the performance of the current first-order estimator (Eq. ( 4))

when the variables highly interact with each other. Indeed, on the top of Fig. 3 we can notice that S SS 2 is not affected. This result is in line with our comments in Section 2.4.

One might propose to circumvent this issue by shifting the vectors of model responses by a factor µ ≃ E [y] before applying Eq. ( 4). This is indeed proposed in [START_REF] Sobol | Monte Carlo estimators for small sensitivity indices[END_REF]. However, by doing so, one would introduce another degree of freedom in the estimator (namely, the value of µ) that would vary from one replicate estimate to another (unless it is left invariant). Such a solution might impact the bias or the variance of the estimator so defined. As argued in [START_REF] Owen | Better estimation of small Sobol' indices[END_REF],

N n=1 y Au n -y B n y A n -y Bu n
(the numerator of Eq. ( 8)) can be seen as a random shifting whereas a constant shifting, namely, [START_REF] Sobol | Monte Carlo estimators for small sensitivity indices[END_REF]. The main finding of [START_REF] Owen | Better estimation of small Sobol' indices[END_REF] is that the former solution outperforms the latter for small first-order Sobol' indices.

N n=1 µ -y B n y A n -µ , is proposed in
Regarding the performance of the total-order estimators, it is not obvious to guess which one is the best. A glance at the plot on the bottom of Fig. 3 reveals that the new estimator has lower variance for ST 3 and higher or equal variances for the two others. One might conclude that the new total-order estimator is more accurate for high total-order Sobol' indices. We investigate this hypothesis further in the next numerical exercise.

The g-function

In this exercise, we study the performance of the two estimators of totalorder Sobol' index. Specifically, we investigate whether the variance of the new estimator is always smaller than the current one or if it depends on the value of ST i . For this purpose, we consider a ten-dimensional function whose total-order Sobol' indices of the input variables spread uniformly over (0, 1). Hence, we consider the Sobol' g-function defined as follows,

f (x) = 10 i=1 |4x i -2| + a i a i + 1
where x i ∼ U (0, 1) for all i = 1, . . . , 10 and the coefficients are chosen as follows: a = (-1.13, -1.24, -1.33, -1.42, -1.52, -1.64, -1.79, -2.00, -2.37, +1.52).

This choice approximately yields the following total-order Sobol' indices, (0.95, 0.85, . . . , 0.15, 0.05). Thus x 1 has the highest total-order effect and

x 10 the lowest. Using negative g-function coefficients a is somewhat unusual but in our case provides a set of total-order Sobol' indices that evenly spreads over (0, 1). In this way, we can investigate numerically whether the performance of the two estimators depends on the magnitude of the total-order Sobol' index. Indeed, the variance of some estimators may depend on the value of the targeted statistic while some may not.

The numerical setting is as follows: we compute one hundred lhs replicate estimates of the total-order sensitivity indices. Samples of size N = 2 is employed (2 21 for the current estimator Eq. ( 5)) in order to get accurate estimates with no overlapping of the ranges of variation. For each estimate, the asymptotic normal variances are evaluated by replacing in Eqs. [START_REF] Lamboni | Multivariate sensitivity analysis: Minimum variance unbiased estimators of the first-order and total-effect covariance matrices[END_REF][START_REF] Lamboni | Uncertainty quantification: A minimum variance unbiased (joint) estimator of the non-normalized Sobol' indices[END_REF][START_REF] Li | Relationship between sensitivity indices defined by variance-and covariance-based methods[END_REF] the exact Sobol' index (i.e. ST i ) and total variance (i.e., V [y]) by their estimated value. The lhs replicates provide also the empirical variances which can be confronted to the asymptotic normal variances.

The one hundred estimates are depicted in Fig. 4 with the exact totalorder Sobol' indices. Despite of the very large sample size, the ranges of variation of the high Sobol' indices estimates are rather large but do not overlap. This indicates that the studied function is a very difficult one for the Monte Carlo estimators. We note that the spread of the IA estimator for the total-order Sobol' index is slightly narrower than the SJ estimator in this example.

On the top of Fig. 5, we represent the estimated variance of the new estimator (namely, τ 2 IA ) versus the variance of the current estimator ( τ 2 SJ ). Because there are one hundred replicates of the sensitivity indices, for each sensitivity index ST i , i = 1, . . . , 10, we have one hundred estimates of the asymptotic normal variances. They are depicted in different coloured circles in the top plot. On the bottom of Fig. 5, we represent the empirical estimated variances obtained by computing directly the variance of the one hundred lhs replicates of each total-order Sobol' index. First, we can note that the y-and x-axes of the two plots (bottom and top) have the same ranges. This indicates that the asymptotic variances Eqs. are good proxies of the empirical variances for the function under study.

The continuous line in Fig. 5 represents τ 2 IA = τ 2 SJ . The scatter plots located below this line mean that τ 2 IA < τ 2 SJ . We observe that the scatter plots associated with the highest sensitivity indices (namely, from ST 1 to ST 4 ) are clearly below this line either for the asymptotic normal variances (top) or the empirical variances (bottom). This confirms that, likewise the Ishigami function, the new estimator Eq. ( 7) is more accurate than Eq. ( 5) at least for high sensitivity indices (say ST i > 0.55). Of course, this inference has been obtained numerically and cannot be generalised.

For the sake of completeness, the results for the first-order Sobol' indices are depicted in Fig. 6. We note that the first-order Sobol' indices are virtually zero which confirm that this case is a very difficult one for global sensitiviy analysis. Both estimators provide results rather centered on the true value.

Although one hundred replicates might not be sufficient to provide stable results, we can notice that, for the first inputs (x 1 , . . . , x 4 ) which have the relatively highest first-order effects, the ranges of variation of the replicates of the Sobol-Saltelli estimator are sligthly narrower than those produced by the IA estimator, for inputs (x 5 , x 6 , x 7 ) their results are rather similar, and for (x 8 , x 9 , x 10 ) (with the smallest sensitivity indices), the IA estimator provides narrower ranges of variation. mate change [START_REF] Charlson | Climate forcing by anthropogenic aerosols[END_REF][START_REF] Sokolik | Direct radiative forcing by anthropogenic airborne minerals aerosols[END_REF]. Radiative forcing models are developed to assess the impact of aerosols. In the present work, we study the direct forcing ∆F by sulfate aerosols in the analytical form provided by [START_REF] Charlson | Climate forcing by anthropogenic aerosols[END_REF],

∆F = - 1 2 S 0 (1 -A c ) T 2 (1 -R s ) 2 βψ e f ψe 3QY L A (16)
where, S 0 is the solar constant, T is the transmittance of the atmospheric layer above the aerosol, A c is the fractional cloud cover, R s is the mean albedo, β is the fraction of the radiation scattered upward by the aerosol, ψ e is the mass scattering efficiency, f ψe is the scaling factor that takes into account the dependence of ψ e to the relative humidity, Q is the global input flux of anthropogenic sulfur, Y is the fraction of SO 2 oxydized to SO 2- 4 , L is the sulfate lifetime in the atmosphere, and A is the area of the Earth [START_REF] Tatang | An efficient method for parametric uncertainty analysis of numerical geophisical model[END_REF].

The negative sign in Eq. ( 16) indicates that the forcing has a cooling effect.

The uncertainties associated with these parameters are taken from [START_REF] Tatang | An efficient method for parametric uncertainty analysis of numerical geophisical model[END_REF] and reported in Tab. 1. The log-normal distribution with geometric mean µ * and geometric standard deviation σ * is denoted LN (µ * , σ * ). If z is a standard normal variable, that is, z ∼ N (0, 1), by setting x = µ * (σ * ) z yields x ∼ LN (µ * , σ * ). According to this transformation, the lhs samples of the uncertain input parameters reported in Tab. 1 have been generated to carry out the sensitivity analysis of the direct forcing model. Note that S 0 and A are treated as parameters as they are known accurately. Therefore, nine uncertain input variables are considered in this study (i.e. D = 9). Samples of size N = 1, 000 have been chosen to perform the analysis as the model is given in an analytical form. This corresponds to a total of N t = 2N (9 + 1) = 20, 000 model runs. The aim of the analysis is to identify i) the irrelevant input variables and ii) possibly the inputs which do not interact with the other ones.

Results

The results are reported in the last column of Tab. 1. They correspond to the point estimate of the first-and total-order Sobol' indices at N = 1, 000.

They are associated with their 95% uncertainty range obtained from the IA estimators' variances as follows: ±1.96 σ IA for S IA and ±1.96 τ IA for

ST

IA . They show that (T, ψ e , Y, L) have first-and total-order Sobol' indices higher than 10%. Their estimated uncertainty ranges overlap which means that it cannot be inferred which parameter is the most important one. The backscattered fraction β has a total-order effect higher than 10% and a firstorder effect slightly lower. The other variables have smaller Sobol' indices but, with the exception of A c , they cannot be neglected ( ST IA i + 1.96 τ IA > 0.05). Actually, the total-order Sobol' index of each variable is virtually the double of the first-order effect. Hence, although the direct forcing model of sulfate aerosols has an apparent simple form, it entails strong interactions between the variables ( S IA i ≃ 0.72) and only one of them can be deemed irrelevant, namely, A c . Fig. 7 depicts the point estimates of the Sobol' indices versus the sample size (up to 1, 000). We can notice that the IA estimators require at least lhs samples of size N = 730 to provide stable results which is quite much for such a low dimensional model (D = 9). This is due, as previously discussed, to the presence of strong interactions between the variables and the relatively high effective dimension of the model (eight out of nine inputs have a totalorder effect greater than 5%). In the denomination of [START_REF] Kucherenko | The identification of model effective dimensions using global sensitivity analysis[END_REF], the direct forcing model of sulfate aerosols can be classified as a Type C model which is a very difficult case for variance-based global sensitivity analysis (see [START_REF] Kucherenko | The identification of model effective dimensions using global sensitivity analysis[END_REF]). [START_REF] Tatang | An efficient method for parametric uncertainty analysis of numerical geophisical model[END_REF]). The last two columns provide the estimated Sobol' indices at N = 1, 000. See body text for further comments.

Input

Conclusions

We have studied the properties of two MC estimators of the first-and total-order Sobol' indices. Their asymptotic variances have been derived

under the asymptotic normality assumption. The so-called IA estimators possess interesting features. One of these features is that the estimated firstorder Sobol' index is always smaller than or equal to the total-order Sobol' index while forming a set of complementary formulas (unlike the current estimators mostly in use by practitioners). By analysing their asymptotic normal variances and by conducting numerical exercises, we have shown that the new sampling strategy and its associated estimators perform better than the current sampling strategy originally introduced in [START_REF] Saltelli | Making best use of model evaluations to compute sensitivity indices[END_REF]. The improvement is especially significant for the first-order Sobol' index estimate. Hence, if one wishes to estimate both the first-and total-order Sobol' indices by Monte Carlo integral approximations, we recommend the use of the IA estimators and the associated sampling design.

  where, V [•] stands for the unconditional variance operator (resp. V [•|•] the conditional variance) and E [•] stands for the mathematical expectation (resp. E [•|•] the conditional expectation). The Sobol' indices range over [0, 1] and

  [START_REF] Cukier | Nonlinear sensitivity analysis of multiparameter model systems[END_REF][START_REF] Gamboa | Sensitivity analysis in general metric spaces[END_REF] require N (D + 2) model calls to estimate the overall set of first-and total-order Sobol' indices while Eqs. (7-8) require 2N (D + 1). To ensure a fair comparison, we set the sample size of the new estimators to half the one of the current estimators. In this way, the computational cost is 2N (D + 1) for the former and 2N (D + 2)

2 .

 2 Run the model with each sample and collect y A = f (x A ) and y B = f (x B ) 3. For all input factors, i = 1 to D a. Generate x A i (and also x B i for the IA estimators) by switching the i-th column of the x A sample matrix and the i-th column of the x B sample matrix b. Run the model and collect y A i = f (x A i ) (and y B i = f (x B i ) for the IA estimators) c. Compute S SS i and ST SJ i from Eqs. (4-5) with the sample set (y A , y B , y A i ) (or ST IA i and S IA i from Eqs. (7-8) with the sample set (y A , y B , y A i , y B i )) by setting u = {i} For a fair comparison, in step 1, the sample size is 2N for the current estimators and N for the new estimators.

Figure 1 :

 1 Figure 1: Boxplot of the 100 lhs replicates of first-and total-order Sobol' indices (resp. at the top and the bottom) with the current and new estimators for the classical Ishigami function. For fair comparison, the sample size is 128 for the current estimators and 64 for the new ones. The boxplot represents the following values: minimum, first-quartile, median, third-quartile and maximum.

Figure 2 :

 2 Figure 2: One hundred lhs replicates of the first-versus total-order Sobol' indices of x 2 obtained with the current and new estimators (blue circles and black crosses). The new estimators provide equal indices at any sample size N as x 2 does not interact with the other variables.

Figure 3 :

 3 Figure 3: Boxplot of the 100 lhs replicates of the first-and total-order Sobol' indices (resp. at the top and the bottom) with the current and new estimators for the modified Ishigami function (f 0 = 100). In this case, the current estimator for first-order Sobol' index performs poorly (top). For the represented values see Fig. 1.

Figure 4 :

 4 Figure 4: Boxplot of the 100 lhs replicates of the total-order Sobol' indices with the current and new estimators for the g-function. For fair comparison, the sample size is 2 20 for the current estimator and 2 21 for the new one.

4 .

 4 GSA of a radiative forcing model4.1. Problem settingAerosol particles influence the Earth's radiative balance directly by backscattering and absorption of solar radiation, thus, contributing to the global cli-

Figure 5 :

 5 Figure 5: Estimated variances of the total-order SI estimators. On the top, by using the asymptotic normal variance formulas, on the bottom, by evaluating the variances of the one hundred lhs replicates. The continuous lines represent τ 2 IA = τ 2 SJ .

Figure 6 :

 6 Figure 6: Boxplot of the 100 lhs replicates of the first-order Sobol' indices with the current and new estimators for the g-function. For fair comparison, the sample size is 2 20 for the current estimator and 2 21 for the new one. Small values indicate that interactions are preponderant in this case study.

Figure 7 :

 7 Figure 7: First-and total-order Sobol' indices estimate versus the sample size N . See discussion in the body text.

  depicts S 2 versus ST 2 for both sets of estimators (the current and new ones). We can see that the pairs ( S IA 2 , ST

					IA 2 ) spread along the
	line S IA 2 = ST	IA 2 contrarily to ( S SS 2 , ST	SJ 2 ). This is also in accordance with
	our findings in § 2.4 that S IA i	= ST	IA i	if x i does not interact with the other
	variables. This is not the case with ( S SS 2 , ST

SJ 2

Table 1 :

 1 Uncertainty assigned to the input parameters of the radiative forcing model (according to

		Baseline PDF	S IA ± 1.96 σ IA ST	IA	± 1.96 τ IA
	S 0 (W.m -2 )	1366	constant	NA	NA
	A (m 2 )	5.1 10 14 constant	NA	NA
	T	0.76	LN (0.76,1.2) 0.13 ± 0.04	0.23 ± 0.06
	1-A c	0.39	LN (0.39,1.1) 0.009 ± 0.003 0.02 ± 0.006
	1-R s	0.85	LN (0.85,1.1) 0.03 ± 0.02	0.06 ± 0.02
	β	0.30	LN (0.30,1.3) 0.07 ± 0.03	0.14 ± 0.04
	ψ e (m 2 .g -1 )	5.0	LN (5.0,1.4) 0.12 ± 0.04	0.23 ± 0.05
	f ψe	1.70	LN (1.70,1.2) 0.03 ± 0.02	0.06 ± 0.03
	Q (10 12 g.yr -1 ) 71	LN (71,1.15) 0.02 ± 0.01	0.04 ± 0.03
	Y	0.5	LN (0.5,1.5) 0.15 ± 0.05	0.27 ± 0.07
	L (days)	5.5	LN (5.5,1.5) 0.16 ± 0.06	0.29 ± 0.09

We denote by S SS u (N ) the estimator for a sample size N . In the sequel, we follow the steps of [START_REF] Janon | Asymptotic normality and efficiency of two sobol index estimators[END_REF] to establish that the asymptotic normality of this estimator is,

with σ 2 SS defined by Eq. [START_REF] Kucherenko | Estimation of global sensitivity indices for models with dependent variables[END_REF].

Proof. We set,

We also denote the associated random vector,

since their statistics do not depend on n.

We then have,

and from Eq. ( 4) we can write,

The so-called Delta method [START_REF] Van Der Vaart | Asymptotic statistics[END_REF] allows for evaluating the variance of the estimator as follows,

by accounting for the definition of (ᾱ, β) above.

Therefore, we find that the variance of this estimator is,

which can be rearranged as follows,

Replacing (α, β) by their expression provides the announced result.

Moreover, by noticing that in Eq. ( 18) α is the numerator of Eq. ( 4) and β the denominator, it is straightforward to demonstrate that the variance of estimator ( 5) is Eq. ( 11). This is merely established by setting α = y Au -y B 2 , β remaining unchanged.

Appendix B Asymptotic normality of S IA u and ST

IA u

In the same way, it can be established that the asymptotic normality of

with σ 2 IA given by Eq. ( 12).

Proof. From Eq. ( 8) we can write,

which yields,

We also denote the associated random vector,

since their statistics do not depend on n.

The so-called Delta method [START_REF] Van Der Vaart | Asymptotic statistics[END_REF] yields,

We find that,

by accounting for the definition of (ᾱ, β, γ) above.

Therefore, we find that the variance of the estimator is,

which can be rearranged as follows,

to finally give,

Furthermore, by replacing (α, β, γ) by their expression we find Eq. [START_REF] Lamboni | Uncertainty quantification: A minimum variance unbiased (joint) estimator of the non-normalized Sobol' indices[END_REF].

By changing (α, β, γ) accordingly we establish the variance of ST IA as,

which is Eq. ( 13).