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

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 liter-ature. 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 [25]. 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. The method to be used depends on the sensitivity indices (also called importance measures) that the analyst wants to compute. As recommended in [25] (see also [24]), 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 ( [28]).
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 [14,10,16,17,13,33]).
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, where, V [·] stands for the unconditional variance operator (resp. V [·|·] the conditional variance) and E [·] stands for the mathematical expectation (resp. indices, the more the group of inputs u is important for the model response. 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 contribute at all to the variance of y. More SA settings are discussed in [24].
There are typically two direct methods to estimate the first-and totalorder Sobol' indices. The first one uses Monte Carlo methods (e.g. [28,20]).
The second one casts the total variance onto orthogonal functions like the Fourier expansion (a.k.a. the Fourier amplitude sensitivity test, [4, 26, 15]) or the polynomial chaos expansion [32,2,27]. 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. [18] 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, 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, where x A n and x B n are two independent input vectors identically distributed, as well as x Au n and x Bu n . The u-values in vector x Au n (resp. x Bu n ) are identical to those in x A n (resp. x B n ) while the v-values are those of x B n (resp. x A n ).

Current estimators
The most popular sampling design to compute simultaneously the firstand total-order sensitivity indices as recommended by Saltelli et al. [23] 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, 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 [30] while Saltelli proposed an estimator similar to the numerator of Eq. (4) in [20].
The superscript SJ often refers to Sobol-Jansen although one can date back the numerator of Eq. (5) to Šaltenis and Dzemyda [36] and Jansen et al. [8] (see [22] 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 [20,21]. As defined, the denominator of Eqs.  at the denominator of Eqs. (4-5) is a positive scalar, we find that, 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 while the numerator is such that, Hence, the following symmetrical estimator for the total-order sensitivity index can be derived, 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 [1].
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 Therefore, the first-order sensitivity index S u is estimated as follows, Besides, and proves that ST 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 [19]. Apart from the factor 2 due to the denominator, the difference is the use by the author of y Cu Interestingly, Lamboni [12] also derived unbiased estimators with minimum variance of the non-normalized Sobol' indices by leaning on the theory of U -statistics (see [6,5]). His construction led to estimators exactly equal to the numerators of Eqs. (7)(8). However, neither Lamboni in [12] nor Owen in [19] 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 [11,12] also satisfy this property but they do not form a set of complementary formulas contrarily to the new estimators defined by Eqs.

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 9 [35,7]. They respectively read as follows, (11) and, . (13) In Eqs. (10)(11)(12)(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 [28]), This implies that, 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 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. (4)(5) require N (D + 2) model calls to estimate the overall set of first-and total-order Sobol' indices while Eqs.

The Ishigami function
Let us consider the following three-dimensional function, 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 In this exercise, we numerically compare the performances of Eqs. (4)(5) with Eqs. (7)(8). 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:  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 .

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, 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  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 [29]. 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 [19], 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, [29]. The main finding of [19] is that the former solution outperforms the latter for small first-order Sobol' indices.
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, 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.

16
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 19 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.

Problem setting
Aerosol particles influence the Earth's radiative balance directly by backscattering and absorption of solar radiation, thus, contributing to the global cli-  22 mate change [3,31]. 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 [3], 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 [34].
The negative sign in Eq. (16) indicates that the forcing has a cooling effect.
The uncertainties associated with these parameters are taken from [34] 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 deterministic input 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.   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 [9], 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 [9]

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 first- We denote by S SS u (N ) the estimator for a sample size N . In the sequel, we follow the steps of [7] to establish that the asymptotic normality of this estimator is, with σ 2 SS defined by Eq. (10).
Proof. We set, (α n , β n ) = 2y A n y Au n − y B n , y A n − y B n 2 .
We also denote the associated random vector, (α, β) = 2y A y Au − y B , y A − y B 2 since their statistics do not depend on n.
We then have, and from Eq. (4) we can write, The so-called Delta method [35] 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. In the same way, it can be established that the asymptotic normality of

Appendix B Asymptotic normality of S
with σ 2 IA given by Eq. (12). We also denote the associated random vector, (α, β, γ) = 2 y Bu − y A y B − y Au , y A − y B 2 , y Au − y Bu 2 since their statistics do not depend on n.
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. (12).