# Robustness of the Sobol' indices to distributional uncertainty

Global sensitivity analysis (GSA) is used to quantify the influence of uncertain variables in a mathematical model. Prior to performing GSA, the user must specific a probability distribution to model the uncertainty, and possibly statistical dependencies, of the variables. Determining this distribution is challenging in practice as the user has limited and imprecise knowledge of the uncertain variables. This article analyzes the robustness of the Sobol' indices, a commonly used tool in GSA, to changes in the distribution of the uncertain variables. A method for assessing such robustness is developed which requires minimal user specification and no additional evaluations of the model. Theoretical and computational aspects of the method are considered and illustrated through examples.

## Authors

• 7 publications
• 5 publications
12/17/2018

### Robustness of the Sobol' indices to marginal distribution uncertainty

Global sensitivity analysis (GSA) quantifies the influence of uncertain ...
08/07/2020

### An information geometry approach for robustness analysis in uncertainty quantification of computer codes

Robustness analysis is an emerging field in the domain of uncertainty qu...
01/15/2018

### An Efficient Method for Uncertainty Propagation in Robust Software Performance Estimation

Software engineers often have to estimate the performance of a software ...
09/25/2019

### Temperature expressions and ergodicity of the Nosé-Hoover deterministic schemes

Thermostats are dynamic equations used to model thermodynamic variables ...
02/19/2019

### New statistical methodology for second level global sensitivity analysis

Global sensitivity analysis (GSA) of numerical simulators aims at studyi...
02/06/2018

### Quantifying dependencies for sensitivity analysis with multivariate input sample data

We present a novel method for quantifying dependencies in multivariate d...
07/27/2021

### A Tale Of Two Long Tails

As machine learning models are increasingly employed to assist human dec...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

Global sensitivity analysis (GSA) aims to quantify the relative importance of the input variables to a function [25]. Such quantification is a crucial step in the development of predictive models. The Sobol’ indices [26, 27] are a commonly used tool for this task, though there are many alternative methods as well, see [25, 13, 8] and references therein for overviews of GSA. Performing GSA in practice consists of:

1. defining a probability distribution for ,

2. evaluating at samples from this distribution,

3. and computing statistics using these evaluations.

Much of the GSA research has focused on steps (ii) and (iii) assuming that step (i) is done by the user. In practice, defining a probability distribution for is a challenging task which carries uncertainties itself. This raises the question, “how robust is my GSA to changes in the distribution of ?” In this article we address this question for the Sobol’ indices.

The robustness of the Sobol’ indices with respect to changes in the marginal distributions of , , is considered in the Life Cycle Analysis literature [16]. Their approach requires the user to specify various admissible marginal distributions for the variables and compute the Sobol’ indices for each possibility. This requires many evaluations of thus limiting its use more broadly. The robustness of the Sobol’ indices with respect to changes in the correlation structure is highlighted in [9] where the authors seek to quantify the risk of ignoring correlations between the variables. Imprecise probabilities are used in [10] to quantify uncertainty in the sensitivity indices. This approach requires the user to parameterize admissible distributions and collect additional samples, i.e. additional evaluations of . In the ecology literature [22]

, the robustness of Sobol’ indices to changes in the means and standard deviations defining normally distributed inputs is examined. The challenges of deep uncertainty, i.e. uncertainty in the distributions of

, , are identified in [7] where the authors compute Sobol’ indices for different distributions and define “robust sensitivity indicators” as a function of the Sobol’ indices from different distributions. Changes in the marginal distributions were shown to change the ordering of the Sobol’ indices in [6]. Similar questions about the robustness of computed quantities with respect to distributional uncertainty may be found in [20, 12, 5, 1].

In [3]

, the functional analysis of variance decomposition (the building block of Sobol’ indices) is analyzed when

does not have a unique distribution but rather multiple possible distributions. The authors provide a framework for analyzing the robustness of the Sobol’ indices which depends upon the user specifying a prior on the space of possible distributions of . In line with [2], the robustness of the Sobol’ indices may be determined by considering a set of possible distributions, sampling from their mixture distribution, and computing the Sobol’ indices with respect to each distribution using a weighting scheme. This approach does not require additional evaluations of , but the user must specify the set of possible distributions, which is challenging in practice.

All of the aforementioned approaches require user specification of possible distributions, additional evaluations of

, or both. In this article, we present a method to measure the robustness of the Sobol’ indices to distributional uncertainties without requiring either. In particular, we consider perturbations of the probability density function (PDF) of

and compute the extreme scenarios when the Sobol’ indices differ most from those computed with the user specified PDF. A judicious formulation allows us to determine these extreme scenarios as the solution of an optimization problem which is solved in closed form. The Sobol’ indices with a perturbed PDF are then computed using weighted averages. Our proposed method is a post processing step which requires minimal user specification and no evaluations of beyond those already taken to compute the Sobol’ indices.

Section 2 provides a review of the Sobol’ indices and establishes notation for the article. Section 3 develops the theoretical and computational framework for our method of assessing robustness. In some cases, the Sobol’ indices may not be robust; however, the importance of the input variables relative to one another may be. This motives us to define and study “normalized Sobol’ indices” in Section 4. In Section 5, we provide an algorithmic description of our method and give a guide for visualizing and interpreting the results. A variety of examples are given in Section 6 to demonstrate the proposed approach and highlight its properties. Section 7 concludes the article by elaborating on limitations and possible extensions of this work.

## 2 Review of Sobol’ Indices

Let , , be a function or model and let be the input variables of that model. The Sobol’ indices [26, 27] measure the importance of a variable (or group of variables) by apportioning to the variable (or group of variables) its relative contribution to the variance of .

Let be a subset of and be its complement. We refer to the group of variables corresponding to as . Following [15], assume that is square integrable and consider the law of total variance decomposition,

 Var(f(X))=Var(E[f(X)|Xu])+E[Var(f(X)|Xu)]. (1)

Using (1), the Sobol’ index and total Sobol’ index for are defined as

 Su=Var(E[f(X)|Xu])Var(f(X))andTu=1−S∼u,

respectively. The Sobol’ index may be interpreted as the proportion of contributed by alone; hence and larger values indicate that is influential. The total Sobol’ index may be interpreted as the proportion of remaining if is known, hence and larger values indicate that is influential. Under the assumption that are independent, and their difference may be interpreted as a measure of the interaction between and . They possess other useful statistical properties and are a preferred method for global sensitivity analysis in many applications. However, much of the statistical theory does not generalize when possess dependencies. In [11], is shown to have a useful approximation theoretic interpretation (with both independent or dependent inputs). Specifically, corresponds to the squared relative error when is optimally approximated by a function which does not depend on . In other words, the influence of is measured by the error when is approximated by a function which does not depend on . Using this interpretation, the total Sobol’ index provides a useful measure of the importance of , with independent or dependent variables. Statistical dependencies in may effect the magnitude of and care must be taken in how one measures the relative importance of the variables. In Section 4, we introduce the normalized Sobol’ index as a means to measure the relative importance of variables when dependencies exist.

We direct the reader to [26, 27, 13, 11, 23, 28, 24, 15, 30, 18]

for a fuller discussion of the Sobol’ indices, total Sobol’ indices, their interpretations, and their estimation. The reader may also consider

[14, 21, 29] for additional discussion of Sobol’ indices with dependent variables, and an alternative approach, the Shapley value.

## 3 Robustness of the Sobol’ Index to PDF Perturbations

Assume that admits a PDF . For simplicity, and because of its approximation theoretic interpretation with dependent variables, we focus on the robustness of the total Sobol’ index to changes in ; the Sobol’ index may be analyzed in a similar fashion. For the remainder of the article, we will use “Sobol’ index” to refer to the total Sobol’ index.

There are multiple ways to express and estimate ; a useful expression from [15] is

 Tu=12∫Ω×Ωu(f(x)−f(x′))2ϕ(x)ϕx|x∼u(x′|x∼u)dxdx′u∫Ωf(x)2ϕ(x)dx−(∫Ωf(x)ϕ(x)dx)2 (2)

where ), ), is the conditional density for , and is the Cartesian product of each , . Note that ) is not a permutation of the entries of but rather a partitioning of them. Then may be estimated by drawing samples from , drawing a second set of samples from , and estimating (2) via Monte Carlo integration of the numerator and denominator separately.

The basic idea of the proposed approach is to view as an operator which inputs the PDF and returns the Sobol’ index. We compute the Fréchet derivative of this operator at and use it to analyze the robustness of . To this end, we make the following assumptions throughout the article:

1. is a Cartesian product of compact intervals,

2. ,

3. is continuous on ,

4. is bounded on .

Some of the results below may be shown with weaker assumptions, these overarching assumptions are made now for conciseness and simplicity. Without loss of generality, under the assumptions above, assume that . We revisit our assumption on the compactness of in Section 7.

We seek to perturb the PDF, so it is essential that the perturbations preserve properties of PDF’s, specifically, that every PDF is non negative and its integral over equals one.

Since is continuous and is compact, is bounded above and below by positive real numbers. Define the Banach space as the set of all bounded functions on equipped with the norm

 ||ψ||V=∣∣∣∣∣∣ψϕ∣∣∣∣∣∣L∞(Ω),

where is the supremum norm on , the set of bounded function on . This norm ensures that for every with , the non negativity property of PDF’s.

To ensure that the integral over equals one, we introduce a normalization operator which takes and returns . Composing this normalization operator with (2) yields the Sobol’ index as an operator on . Define by

 F(η)=12∫Ω×Ωu(f(x)−f(x′))2η(x)η(x′)1∫Ωuη(x)dxudxdx′u, (3)
 G(η)=∫Ωf(x)2η(x)dx−1∫Ωη(x)dx(∫Ωf(x)η(x)dx)2, (4)
 Tu(η)=F(η)G(η). (5)

It is easily observed that multiplying the numerator and denominator of (5) by yields that (5) and (2) coincide with replaced by . In this framework, every such that is nonnegative and corresponds to the Sobol’ index computed with respect to the PDF .

Having defined the Sobol’ index as an operator which inputs bounded PDF’s, Theorem 3.1 below gives the Fréchet derivative of the Sobol’ index at .

###### Theorem 3.1

The operator is Fréchet differentiable at with Fréchet derivative given by the bounded linear operator

 DTu(ϕ)ψ=DF(ϕ)ψG(ϕ)−Tu(ϕ)DG(ϕ)ψG(ϕ), (6)

where

 DF(ϕ)ψ= 12∫Ω×Ωu(f(x)−f(x′))2ψ(x′)ϕ(x′)ϕ(x)ϕx|x∼u(x′|x∼u)dxdx′u + 12∫Ω×Ωu(f(x)−f(x′))2ψ(x)ϕ(x)ϕ(x)ϕx|x∼u(x′|x∼u)dxdx′u − 12∫Ω×Ωu(f(x)−f(x′))2∫Ωuψ(x)dxu∫Ωuϕ(x)dxuϕ(x)ϕx|x∼u(x′|x∼u)dxdx′u

and

 DG(ϕ)ψ= ∫Ωf(x)2ψ(x)ϕ(x)ϕ(x)dx −2∫Ωf(x)ϕ(x)dx∫Ωf(x)ψ(x)ϕ(x)ϕ(x)dx +(∫Ωψ(x)ϕ(x)ϕ(x)dx)(∫Ωf(x)ϕ(x)dx)2.

A proof for Theorem 3.1 is given in the appendix. Without the assumption that is compact (or at least bounded), an infinitesimal perturbation may yield that is not integrable. Hence it is a theoretical necessity to assume that is bounded. In Section 7, we revisit this discussion and highlight that this assumption is less restrictive in practice.

If the Sobol’ index is computed using Monte Carlo estimators of (2), then may be estimated using these samples and evaluations of ; the only additional work is evaluating and at the sample points. Hence may be estimated at any with negligible computational cost. This is why, as previously mentioned, our method is a post processing step which requires no additional evaluations of beyond those taken to compute the Sobol’ indices.

We seek to find an “optimal” perturbation of in the sense that it causes the greatest change in the Sobol’ index. The locally optimal perturbation is the , , which maximizes . To estimate this , we define a finite dimensional subspace and compute the operator norm of the restriction of to . When choosing , there is a trade off to consider between the approximating properties of functions from , our ability to use existing samples to estimate the action of on functions from , and the ease of computing the operator norm of restricted to . In what follows, we choose to be a subspace generated by the span of a set of locally supported piecewise constant functions.

Let , , be a partition of into open hyperrectangles, i.e. and for ; denotes the closure of . Define

 ψi(x)={1x∈Ri0x∉Ri

to be the indicator function of , , and , a dimensional subspace of . The partition may be efficiently constructed using Regression Trees [4]; we will elaborate on this in Section 5.

Constructing in this way has computational and approximation theoretic advantages. Its computational advantage, demonstrated below, is that it enables a closed-form solution to an otherwise challenging optimization problem. Its approximation theoretic advantage is that is provides a mechanism to constrain the subspace with the existing evaluations of , see Section 5 for more details. Piecewise constant functions defined on a partition of the domain are useful in the sense that they permit perturbations of any functional form, constrained by the coarseness of the partition. The coarseness of the partition depends on the number of existing evaluations of , hence the form of functions from are constrained by the existing data, not the users imposition of functional forms. The proposed subspace is optimal in the (informal) sense that is provides the most flexibility in functional forms given the constraint of existing data.

The operator norm of restricted to is given by

 ||DTu(ϕ)||L(VM,R) =maxψ∈VM||ψ||V≤1|DTu(ϕ)ψ| =maxa∈RM||∑Mi=1aiψi||V≤1∣∣ ∣∣DTu(ϕ)(M∑i=1aiψi)∣∣ ∣∣ =maxa∈RM||∑Mi=1aiψi||V≤1∣∣ ∣∣M∑i=1aiDTu(ϕ)ψi∣∣ ∣∣

Since the basis functions have disjoint support, it follows that

 ∣∣ ∣∣∣∣ ∣∣M∑i=1aiψi∣∣ ∣∣∣∣ ∣∣V=∣∣ ∣∣∣∣ ∣∣M∑i=1ai1ϕψi∣∣ ∣∣∣∣ ∣∣L∞(Ω)=maxi=1,2,…,M|ai|∣∣∣∣∣∣1ϕ∣∣∣∣∣∣L∞(Ri),

which implies

 ∣∣ ∣∣∣∣ ∣∣M∑i=1aiψi∣∣ ∣∣∣∣ ∣∣V≤1

is equivalent to , where is the infimum of on , .

Let be defined by , . Then we have

 ||DTu(ϕ)||L(VM,R)=maxa∈RM|ai|≤bii=1,2,…,M|dTa|.

This problem may be solved in closed form to get

 ai=sign(di)bi

and

 ||DTu(ϕ)||L(VM,R)=||d||1.

In what follows, we refer to , , which maximizes the Fréchet derivative, as the optimal perturbation. Finding the optimal perturbation and the corresponding operator norm simplifies to evaluating for , which may be estimated with negligible additional computation. However, estimating is typically more challenging than estimating the Sobol’ index. Rather than inferring robustness with , we propose to:

1. estimate , ,

2. use weighted averaging with the existing evaluations of and to estimate the Sobol’ indices with respect to the optimally perturbed PDF, which we define as

 ϕ+δM∑i=1aiψi1+δM∑i=1aivol(Ri), (7)

where is a parameter to scale the size of the perturbation and is the volume of the set ; the determination of will be discussed in Section 5. We will refer to the Sobol’ indices computed with as the nominal Sobol’ indices and the Sobol’ indices computed with the optimally perturbed PDF (7) as the perturbed Sobol’ indices.

In practice, it is suggested to estimate the terms , , in (7) with a Monte Carlo estimator from the existing data. They may be computed analytically since is known; however, if they are computed exactly then the weights used to estimate perturbed Sobol’ indices may not sum to one because of Monte Carlo error in the estimate. This can bias the resulting analysis. Estimating , , from the existing data diminishes this potential bias.

Our weighted averaging approach is an improvement from traditional derivative based robustness analysis in several ways:

1. Estimating is easy. Since is known, may be computed numerically by querying the existing evaluations of (or possibly analytically, for instance if then for every ), we consider this negligible. Assuming that we have enough samples for the Sobol’ index estimation to converge, determining the sign of with these samples is relatively easy. Additionally, when we do not determine the sign of correctly it is frequently because , in which case this error is benign in the scope of our analysis.

2. The user must determine ; however, various values of may be tested at negligible computation cost. The sample standard deviation of the weighted average may be compared with the sample standard deviation in the original estimator to determine admissible values of . Additional details are given in Section 5.

3. Computing the perturbed Sobol’ indices estimates a realized worst case. This is superior to worst case bounds, error bars, or confidence intervals, which in many cases are overly pessimistic. Further, computing error bars for each Sobol’ index individually may yield misleading results. For instance, error bars for two variables may yield large intervals for each Sobol’ index, but their magnitude relative to one another is nearly constant for any PDF perturbation. In this case the user would incorrectly conclude that the relative importance of the variables to one another is uncertain.

As previously highlighted, one way to test for robustness is to use weighted averages to estimate the Sobol’ indices with different PDF’s. The challenge with this approach is that the user must specify the perturbed PDF’s. Our method may be viewed as an improvement on this idea by automating the choice of perturbed PDF’s. The Fréchet derivative operator norm yields a locally optimal perturbation, which will likely reveal greater changes in the Sobol’ indices when compared with a user manually selecting a small set of perturbed PDF’s. However, our method does not have the danger of finding unrealistic worst cases since it only seeks perturbations in a neighborhood of the existing PDF and is constrained to use the existing samples.

## 4 Robustness of the Normalized Sobol’ Index to PDF Perturbations

As highlighted in [11], the Sobol’ indices are frequently smaller when possesses stronger dependencies. Since the magnitude of the Sobol’ indices may change when the PDF is perturbed, we seek to analyze the robustness of the relative importance of the variables rather than the magnitude of the Sobol’ indices. This is useful in practice when, for instance, a modeler (for lack of better knowledge) assumes the variables are independent, computes the Sobol’ indices to rank the importance of the variables, but is uncertain of their ranking because dependencies that were potentially ignored. For clarity and notational simplicity, the remainder of the article will focus on the Sobol’ indices when is a singleton, i.e. the set of Sobol’ indices .

To measure the relative importance of the variables as the PDF varies, define the normalized Sobol’ index as

 Tk(ϕ)=Tk(ϕ)p∑i=1Ti(ϕ) (8)

for . The example below illustrates the behavior of the Sobol’ indices and normalized Sobol’ indices.

### Example

Let

 f(X)=1.5X1+1.25X2+X3 (9)

and follow a multivariate normal distribution with mean and covariance matrix given by

 μ=⎡⎢⎣000⎤⎥⎦,Σ=⎡⎢⎣1ρρρ1ρρρ1⎤⎥⎦,0≤ρ≤1.

Figure 1 shows the Sobol’ indices and normalized Sobol’ indices , , as a functions of . As increases (the correlations strengthen), the Sobol’ indices decrease while the normalized Sobol’ indices are constant. The trend of decreasing as correlations increase is general; [11] explains it with an approximation theoretic perspective of the Sobol’ indices. The normalized Sobol’ indices are constant indicating that, though the Sobol’ indices decrease, the relative importance of the variables does not change. In general the normalized Sobol’ indices will change as the distribution of the input variables changes.

Applying Theorem 3.1 and the quotient rule to (8) yields that is Fréchet differentiable with Fréchet derivative

 DTk(ϕ)ψ=(p∑i=1Ti(ϕ))DTk(ϕ)ψ−Tk(ϕ)(p∑i=1DTi(ϕ)ψ)(p∑i=1Ti(ϕ))2.

Since is a linear combination of the operators from Section 3, we may easily estimate using the same results previously presented. In fact, in Section 3 a subspace is defined and we compute for . Using this computation we may easily compute , , at no additional cost. The same procedure from Section 3 may be adopted to compute perturbed PDF’s and perturbed Sobol’ indices via weighted averaging. Since the cost is negligible, it is suggested to compute the optimal perturbation using and , and estimate the perturbed Sobol’ indices for each perturbation.

Definition 1 below aids to identify perturbations which change the Sobol’ indices but not the relative importance of the variables.

###### Definition 1

Let and denote the perturbed Sobol’ indices and perturbed normalized Sobol’ indices, respectively, for some perturbation of . The absolute change in the Sobol’ indices is

 p∑k=1|Tk−~Tk|

and the relative change in the Sobol’ indices is

 p∑k=1|Tk−~Tk|.

It is suggested to consider both the perturbed Sobol’ indices which yield the largest absolute and relative changes. This will be further described in Section 5 and demonstrated in Section 6. We emphasize that the normalized Sobol’ indices are a tool to find perturbations changing the relative importance of the variables, but the user will typically use the traditional Sobol’ indices to make inferences.

## 5 Algorithmic Description

Algorithm 1 below summarizes our proposed method. In this section, we discuss the user inputs of Algorithm 1 in detail, highlight important algorithmic features of our method, and consider the visualization and interpretation of the results.

It was previously suggested to generate the partition , , with a Regression Tree [4]. This is a judicious choice because the minimum number of samples in the sets is easily specified. An integer may be input and the Regression Tree will recursively partition ensuring that each set of the partition contains at least samples. This simplicity make Regression Trees attractive in our context. Taking small values of typically results in being a larger subspace, but will create error when estimating (since there will be fewer samples to estimate the integrals). The determination of is discussed below. The relationship between and depends on the algorithm used to generate the partition; a Regression Tree will uniquely determine as a function , typically a decreasing function of .

As highlighted in Section 3, piecewise constant functions from permit general functional forms for the PDF perturbations. Because the partition size is constrained by the existing evaluation of , it will typically be coarse by approximation theoretic standards. In order to find the largest possible changes in the Sobol’ indices, constrained by the coarseness of the partition, it is advantageous to partition the domain finely in regions where varies more. This is precisely what a Regression tree, trained to predict , seeks to do. The Regression Tree algorithm recursively partitions the domain, through a greedy algorithm, which minimizes the difference between and a piecewise constant approximation at each iteration. It begins with the entire domain and refines the partition by considering splits along the directions of the coordinate axes. When is large it may not split in every coordinate direction. This is acceptable, and in many cases beneficial, as it adapts the partition according to the variability of .

In some cases, as illustrated in Subsection 6.3, partitioning the domain according to the variability of may result in sets , , where most of the ’s are small. This is problematic because it limits the size of admissible perturbations. To mitigate this, a Regression Tree may be trained to generate a coarser partition which can be refined by the user to ensure that only a few ’s are small. We discuss this further in Subsection 6.3.

The norm of the perturbed PDF in (7) depends on . It was suggested to try various values of (equally spaced points in ) and accept those which meet a convergence tolerance. If Monte Carlo integration is used to estimate the Sobol’ indices , then the sample standard deviation may be used as a metric for convergence. Let and , , denote the sample standard deviation for the nominal and perturbed Sobol’ indices, respectively. For the results presented in this article, the sample standard deviation is estimated by computing the standard derivation of 50 estimates generated by randomly subsampling half of the function evaluations. Assuming that , , are sufficiently small to ensure convergence of the nominal Sobol’ indices, it is required that be less than a threshold. Define

 t=maxj=1,2,…,p(~σj/~Tj)(σj/Tj)

and specify a threshold . The perturbed Sobol’ indices are accepted if .

The inputs of Algorithm 1 are:

1. , the number of Monte Carlo samples,

2. , the minimum number of samples in each set of the partition,

3. , an integer denoting how many values of to consider,

4. and , the acceptance threshold for the perturbed Sobol’ indices.

The results in Section 6 use , , and ; the number of Monte Carlo samples required depends on the problem. Numerical evidence, and intuition, indicate that is approximately a quadratic function of centered at . To determine , we may solve the scalar nonlinear equation by evaluating at equally spaced points in . It is not necessary to take large values for ; the choice introduces negligible computation and provides sufficient resolution for our purposes. The choice is considered a reasonable threshold to permit non trivial perturbations without introducing significant numerical errors. Our choice of is the least intuitive of the inputs. To justify this choice, a numerical experiment was performed varying , . The results, omitted from this article for conciseness, indicate that our method is robust to changes in . If necessary, the user may easily verify the particular choice of inputs used in their application by varying them. The computational cost of this numerical experiment is small.

Lines 2-5 of Algorithm 1 is the Sobol’ index estimation and Lines 6-17 is our robustness analysis. In many applications, Line 4 dominates the computational cost and hence the cost of robustness analysis is negligible. Lines 6 and 8 may be done analytically in many applications. The computation in Lines 9-18 is primarily taking sample averages of data on memory so its cost is small. In particular, the nested for loops may appear burdensome, but the operations inside of them are sufficiently simple that they may be executed quickly.

Algorithm 1 returns a collection of sets perturbed Sobol’ indices. We suggest extracting the perturbed Sobol’ indices with the largest absolute and relative changes to visualize alongside the nominal Sobol’ indices, denote them as where the superscripts and identify the Sobol’ indices with largest absolute and relative changes, respectively. This may be done by querying the collection of perturbed Sobol’ indices and creating a bar plot of , see Figure 3 for an illustration of this. There are several possible scenarios the user may observe:

1. If , , then the user may confidently make inferences with the Sobol’ indices.

2. If , , but , , then the user may confidently make inferences about the relative importance of the variables but not the magnitude of the Sobol’ indices.

3. If there are variables such that then they may be considered unimportant.

4. If but then the user should excise caution treating as unimportant.

5. If but then the user may not be certain of the importance of and relative to one another.

If a particular Sobol’ index is of interest, the collection of perturbed Sobol’ indices may be queried to asses its robustness. The user may easily visualize all of the perturbed indices in a histogram.

## 6 Numerical Results

In this section, three examples are presented to highlight different properties of our proposed method. The first example analyzes how our robustness analysis changes as more samples are collected. The second example expands on Section 4 by highlighting a case when the largest absolute change in the Sobol’ indices yields a small relative change. The final example is an application of our method to the Lorenz system [17]. We consider two cases in this example to demonstrate the effect of the partitioning on our robustness analysis.

### 6.1 g-function Example to Demonstrate Convergence in Samples

Let

 f(X)=10∏k=1|4Xk−2|+ak1+ak, (10)

where each

is independent and uniformly distributed on

, and for . This is the “g-function” [27] commonly used in the GSA literature.

We compute the nominal Sobol’ indices and perturbed Sobol’ indices of (10). The number of Monte Carlo samples is varied to analyze the convergence behavior of our robustness estimation, specifically, we use 1,000, 5,000, 10,000, and 50,000 Monte Carlo samples. For each fixed sample size, 32 repetitions of the calculation is performed to understand sampling variability. Figure 2 below displays box plots for the estimation of the largest Sobol’ index, . The center panel is our estimation of

; the median estimation is nearly constant and the quantiles shrink as the number of samples increases, this reflects convergence of the estimation. The perturbation size

is varied between -1 and 1 and it is determined that is the maximum admissible perturbation size for the threshold . The left and right panels show the convergence of with perturbations and , respectively. The shrinking quantiles are very similar to those in the center panel demonstrating that the estimation error in is comparable to the estimation error in . The left and right panels have slight decreasing and increasing trends, respectively. This is because the subspace is larger when more samples are taken, thus the perturbations yield larger changes in the Sobol’ indices. For this example, the trend is relatively small reflecting the fact that taking a larger subspace does not yield significant changes in the Sobol’ index.

### 6.2 Linear Example to Demonstrate the Normalized Sobol’ Indices

This example illustrates the difference in the largest absolute and relative perturbations of the Sobol’ indices. Let be defined by (9) and each be independent and uniformly distributed on , . The Sobol’ indices are estimated with 5,000 Monte Carlo samples. Figure 3 displays the nominal Sobol’ indices of in blue, the perturbed Sobol’ indices with the largest absolute differences change in cyan, and the perturbed Sobol’ indices with the largest relative change in yellow.

The largest absolute change of the Sobol’ indices corresponds to the case when they are all shifted down but their relative importance does not change. The largest relative change identifies a case where decreases while and increase. The relative importance of the variables change with this perturbation, demonstrating the benefit of considering the largest absolute and relative perturbations.

### 6.3 Lorenz System

This example applies our method to the well known Lorenz system [17], a model for atmospheric convection. Sobol’ indices were considered for this system in [19]

. The Lorenz system is described by the system of ordinary differential equations

 dy1dt=σ(y2−y1) dy2dt=y1(ρ−y3)−y2 dy3dt=y1y2−βy3

with initial conditions , . Letting denote the uncertain parameters, we compute the Sobol’ indices of the function

 f(X)=y3(1)y2(1), (11)

the ratio of the states and at time . This choice of corresponds to a ratio of temperature variations after a duration of 1 time unit.

The distribution of is chosen to reflect uncertainty about nominal values of the parameters. Two different cases, in the sub-subsections below, are considered to highlight different features of our method. For each case, Monte Carlo samples are taken for the Sobol’ index estimation.

#### 6.3.1 Lorenz System Case 1

In this first case we assume the parameters are independent with the uniform distributions given in Table 1 below. Figure 4 displays the nominal Sobol’ indices in blue, the perturbed Sobol’ indices with the largest absolute change in cyan, and the perturbed Sobol’ indices with the largest relative change in yellow. Several inferences may be drawn from this result,

1. and are the most influential parameters, although their Sobol’ indices and relative importance is uncertain,

2. the Sobol’ indices for , , and and their importance relative to one another is robust,

3. has little influence and its small Sobol’ index is robust, it may be considered a non-influential parameter.

#### 6.3.2 Lorenz System Case 2

In this second case we assume the parameters are independent and that all parameters have the same marginal distribution given in Table 1 with the exception of . Instead of being uniformly distributed on as in Case 1, we take

to have a Beta distribution on

with shape parameters 111For shape parameters

, a Beta random variable

on has PDF .
. This corresponds to giving greater probability to .

A partition is generated by training a Regression Tree to predict . The left panel of Figure 5 displays the nominal Sobol’ indices in blue, the perturbed Sobol’ indices with the largest absolute change in cyan, and the perturbed Sobol’ indices with the largest relative change in yellow. The results indicate that the Sobol’ indices are robust, a different conclusion than was reached in Case 1. This occurs because the Regression Tree never partitioned on so each set contained the entire support of . Because the marginal PDF for takes small values on part of its support, namely near , the infimum of on each is small. The partition generated by the Regression Tree yielded very small perturbations and as a result did not produce significant changes in the Sobol’ indices.

To alleviate this problem, a partition is generated by a Regression Tree trained to predict using all of the variables except . A minimum of samples are requested in each hyperrectangle rather than , as requested previously. This yields a coarser discretization of the other 5 variables. The resulting partition is refined by splitting each set into 4 subsets defined by partitioning at the quantiles of . This yields a partition with approximately samples per subset and a sufficient discretization of to enable larger perturbations. Figure 5 displays the nominal Sobol’ indices in blue, the perturbed Sobol’ indices with the largest absolute change in cyan, and the perturbed Sobol’ indices with the largest relative change in yellow. Larger changes in the Sobol’ indices are observed, as is expected. However the changes are smaller than what was observed in Case 1. This is because the partition used in Case 1 was generated by a Regression Tree which better approximated , and hence allowed for larger perturbations of the Sobol’ indices. The general conclusion from this example is that the partition should be generated so that the Regression Tree approximates as well as possible. If small values of prohibit taking large perturbations, then the partition may be generated with fewer hyperrectangles, followed by a refining of this coarse partition to sufficiently discretize the necessary regions. This may result in a failure to discover the largest possible perturbations, as demonstrated by comparing Case 1 and Case 2.

## 7 Conclusion

This article presents a novel framework in which robustness of the Sobol’ indices with respect to the input variables distribution may be assessed. The proposed method permits such analysis to be done at negligible computational cost. For a modeler using Sobol’ indices, this robustness analysis can be obtained as a by-product of computing Sobol’ indices and may be easily visualized along with the indices themselves. Understanding the robustness of the Sobol’ indices to distributional uncertainty prevents the user from making incorrect inferences which have significant consequences. For instance, reducing dimensions by fixing variables with small Sobol’ indices—which are not robust—may result in model variations which are not explained in the lower dimensional space.

The method suffers four primary limitations, namely,

1. the nominal PDF must be compactly supported,

2. perturbations may not change the support of the nominal PDF,

3. the perturbations are only locally optimal,

4. generating a partition is difficult if the distribution of is far from being uniform.

The first limitation prohibits a direct application of our method to many commonly used PDF’s. This occurs because the Fréchet derivative is not well defined if we allow perturbations in the tail of the distribution. However, one may defined a compact subset of the domain where the PDF assumes most of its mass and allow perturbations on this subset while keeping the tail fixed. The compact subset may be chosen large enough that this truncation error is practically irrelevant, for instance, having the probability of the tails less than machine epsilon. This limitation is primarily theoretical and is not a significant practical concern. The greater limitation will be regions of small probability within the compact subset, see the fourth limitation. A theory analogous to what is presented in the article may be developed when the domain is truncated. In practice, the modeler will take a compact subset which contains all of the existing samples and, on the discrete level, the robustness analysis will be identical to what is presented in this article.

The second limitation occurs because we have formulated the method to work with existing samples. If the support of the PDF increases, then we would need additional evaluations of in these unexplored regions.

The third limitation arises because the perturbation direction is determined by maximizing a derivative, which is local. If the PDF to Sobol’ index mapping is highly nonlinear, this may not be an adequate. However, finding a globally optimal perturbation requires far more computational effort. A locally optimal perturbation is useful and appropriate, particularly for its computational advantages.

The fourth limitation, as demonstrated in Subsection 6.3, may arise when marginal distributions differ significantly from being uniform. If a small number of marginals do so, this limitation may be mitigated by the approach described in Subsection 6.3. There is ongoing work to develop a variant of our proposed method which removes this limitation by partitioning and taking perturbations on the marginal distributions separately. The method proposed in this article will be most effective when the distribution of is approximately uniform, which is a common occurrence in applications where statistical information is poorly known but bounds may be provided from the physics of the problem.

## Appendix

Proof of Theorem 3.1.

###### Proof

One may easily observe that in a neighborhood of (assuming is non constant). It is sufficient to compute the Fréchet derivatives of and , the Fréchet derivative of follows from the quotient rule. The Fréchet derivatives of

 ∫Ωf(x)ϕ(x)dx,∫Ωf(x)2ϕ(x)dx, and ∫Ωϕ(x)dx,

when considered as operators from to , acting on , are easily shown to be

 ∫Ωf(x)ψ(x)dx,∫Ωf(x)2ψ(x)dx, and ∫Ωψ(x)dx,

respectively, using the definition of the Fréchet derivative. The Fréchet derivative of

follows from the sum/difference, product, and chain rule of differentiation.

The Fréchet derivative of may be computed by first defining an operator
,

 H(η)=η(x)η(x′)1∫Ωuη(x)dxu,

where . The Fréchet derivatives of

 η(x),η(x′), and% ∫Ωuη(x)dxu,

when considered as operators from to , acting on , are easily shown to be

 ψ(x),ψ(x′), and% ∫Ωuψ(x)dxu,

respectively, using the definition of the Fréchet derivative. The Fréchet derivative of follows from the product and quotient rules of differentiation. The Fréchet derivative of may be easily computed using the Fréchet derivative of , the boundedness of , and the chain rule of differentiation.

## References

• [1] R. J. Beckman and M. D. McKay, Monte Carlo estimation under different distributions using the same simulation, Technometrics, 29 (1987), pp. 153–160.
• [2] R. Bolado-Lavin and A. C. Badea, Review of sensitivity analysis methods and experience for geological disposal of radioactive waste and spent nuclear fuel, tech. report, European Commission, 2009.
• [3] E. Borgonovo, M. D. Morris, , and E. Plischke, Functional ANOVA with multiple distributions: Implications for the sensitivity analysis of computer experiments, SIAM/ASA J. Uncertain. Quantif., 6 (2018), pp. 397–427.
• [4] L. Breiman, J. Friedman, R. Olshen, and C. Stone, Classification and regression trees, Wadsworth Advanced Books and Software, 1984.
• [5] S. E. Chick, Input distribution selection for simulation experiments: Accounting for input uncertainty, Operations Research, 49 (2001), pp. 744–758.
• [6] A. Cousin, A. Janon, V. Maume-Deschamps, and I. Niang, On the consistency of Sobol’ indices with respect to stochastic ordering of model parameters
• [7] L. Gao, B. A. Bryan, M. Nolan, J. D. Connor, X. Song, and G. Zhao, Robust global sensitivity analysis under deep uncertainty via scenario analysis, Environmental Modelling & Software Software, 76 (2016), pp. 154–166.
• [8] R. Ghanem, D. Higdon, and H. Owhadi, eds., Handbook of Uncertainty Quantification, Springer, 2017.
• [9] E. Groena and R. Heijungs, Ignoring correlation in uncertainty and sensitivity analysis in life cycle assessment: what is the risk?, Environmental Impact Assessment Review, 62 (2017), pp. 98–109.
• [10] J. Hall, Uncertainty-based sensitivity indices for imprecise probability distributions, Reliability Engineering and System Safety, 91 (2006), pp. 1443–1451.
• [11] J. Hart and P. Gremaud, An approximation theoretic perspective of Sobol’ indices with dependent variables, International Journal for Uncertainty Quantification, 8 (2018), pp. 483–493. https://arxiv.org/abs/1801.01359.
• [12] Z. Hu, J. Cao, and L. J. Hong, Robust simulation of global warming policies using the dice model, Management Science, 58 (2012), pp. 2190–2206.
• [13] B. Iooss and P. Lemaître, A review on global sensitivity analysis methods, in Uncertainty management in simulation-optimization of complex systems, G. Dellino and C. Meloni, eds., Springer, 2015, ch. 5, pp. 543–501.
• [14] B. Iooss and C. Prieur, Shapley effects for sensitivity analysis with dependent inputs: comparisons with sobol’ indices, numerical estimation and applications, https://hal.inria.fr/hal-01556303/file/RESS17-ioossPrieur.pdf.
• [15] S. Kucherenko, S. Tarantola, and P. Annoni, Estimation of global sensitivity indices for models with dependent variables, Computer Physics Communications, 183 (2012), pp. 937–946.
• [16] M. Lacirignola, P. Blanc, R. Girard, P. Pérez-López, and I. Blanc, LCA of emerging technologies: addressing high uncertainty on inputs’ variability when performing global sensitivity analysis, Science of the Total Environment, 578 (2017), pp. 268–280.
• [17] E. N. Lorenz, Deterministic non-periodic flow, Journal of the Atmospheric Sciences, 20 (1963), pp. 130–141.
• [18] T. Mara and S. Tarantola, Variance-based sensitivity analysis of computer models with dependent inputs, Reliability Eng. Sys. Safety, 107 (2012), pp. 115–121.
• [19] C. Marzban, Variance-based sensitivity analysis: An illustration on the Lorenz’63 model, American Meteorological Society, (2013), pp. 4069–4079.
• [20] A. Millner, S. Dietz, and G. Heal, Scientific ambiguity and climate policy, Environ Resource Econ, 55 (2013), pp. 21–46.
• [21] A. Owen and C. Prieur, On shapley value for measuring importance of dependent inputs, SIAM/ASA J. Uncertain. Quantif., 5 (2017), pp. 986–1002.
• [22] L. Paleari and R. Confalonieri, Sensitivity analysis of a sensitivity analysis: We are likely overlooking the impact of distributional assumptions, Ecological Modelling, 340 (2016), pp. 57–63.
• [23] C. Prieur and S. Tarantola, Variance-based sensitivity analysis: Theory and estimation algorithms, in Handbook of Uncertainty Quantification, R. Ghanem, D. Higdon, and H. Owhadi, eds., Springer, 2017.
• [24] A. Saltelli, P. Annoni, I. Azzini, F. Campolongo, M. Ratto, and S. Tarantola, Variance based sensitivity analysis of model output. design and estimator for the total sensitivity index, Computer Physics Communications, 181 (2010), pp. 259–270.
• [25] A. Saltelli, M. Ratto, T. Andres, F. Campolongo, J. Cariboni, D. Gatelli, M. Saisana, and S. Tarantola, Global sensitivity analysis: the primer, Wiley, 2008.
• [26] I. Sobol’, Sensitivity estimates for non linear mathematical models, Math. Mod. Comp. Exp., 1 (1993), pp. 407–414.
• [27] I. Sobol’, Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates, Mathematics and Computers in Simulation, 55 (2001), pp. 271–280.
• [28] I. Sobol’, Theorems and examples on high dimensional model representation, Reliability Eng. Sys. Safety, 79 (2003), pp. 187–193.
• [29] E. Song, B. L. Nelson, and J. Staum, Shapley effects for global sensitivity analysis: Theory and computation, SIAM/ASA J. Uncertain. Quantif., 4 (2016), pp. 1060–1083.
• [30] S. Tarantola and T. A. Mara, Variance-based sensitivity indices of computer models with dependent inputs: The Fourier Amplitude Sensitivity Test, International Journal for Uncertainty Quantication, 7 (2017), pp. 511–523.