Robustness of the Sobol' indices to marginal distribution uncertainty

TL;DR

This paper introduces a new method to assess how sensitive Sobol' indices are to uncertainties in the marginal probability distributions of input variables in global sensitivity analysis, enhancing the robustness evaluation.

Contribution

It proposes a novel approach using optimal perturbations of marginal distributions to analyze Sobol' indices' robustness under distributional uncertainty.

Findings

Method effectively quantifies robustness of Sobol' indices

Illustrated with synthetic examples and a contaminant transport model

Provides insights into distributional assumptions impact on sensitivity analysis

Abstract

Global sensitivity analysis (GSA) quantifies the influence of uncertain variables in a mathematical model. The Sobol' indices, a commonly used tool in GSA, seek to do this by attributing to each variable its relative contribution to the variance of the model output. In order to compute Sobol' indices, the user must specify a probability distribution for the uncertain variables. This distribution is typically unknown and must be chosen using limited data and/or knowledge. The usefulness of the Sobol' indices depends on their robustness to this distributional uncertainty. This article presents a novel method which uses "optimal perturbations" of the marginal probability density functions to analyze the robustness of the Sobol' indices. The method is illustrated through synthetic examples and a model for contaminant transport.

Tables

Table 1. Table 1: Nominal values for parameters for the advection diffusion example. Each parameter is assumed to be independent and uniformly distributed with ± 30 % plus-or-minus percent 30 \pm 30\% uncertainty about these nominal values.

Parameter	$ϵ$	${\alpha }_1$	${\alpha }_2$	${\xi }_1$	${\xi }_2$	$\gamma$	$\beta$	${\nu }_1$	${\nu }_2$
Nominal Value	10	210	70	0.5	0.5	50	100	0.1	0.2

Equations

$$f(\mathbf{X})=f_0+\sum _{i=1}^p{f}_i(X_i)+\sum _{i=1}^p\sum _{j=i+1}^p{f}_{i,j}(X_i,X_j)+\cdots +f_{1,2,\dots ,p}(\mathbf{X})$$

$$f_0=\mathbb{E}[f(\mathbf{X})],$$

$$f_i(X_i)=\mathbb{E}[f(\mathbf{X})|X_i]-f_0,$$

$$f_{i,j}(X_i,X_j)=\mathbb{E}[f(\mathbf{X})|X_i,X_j]-f_i(X_i)-f_j(X_j)-f_0,$$

$$⋮$$

$$f_u({\mathbf{X}}_u)=\mathbb{E}[f(\mathbf{X})|{\mathbf{X}}_u]-\sum _{v\subset u}{f}_v({\mathbf{X}}_v).$$

$$\mathbb{E}[f_u({\mathbf{X}}_u)f_v({\mathbf{X}}_v)]=0\text{when}u\ne v.$$

$$\mathrm{Var}(f(\mathbf{X}))=\sum _{i=1}^p\mathrm{Var}(f_i(X_i))+\sum _{i=1}^p\sum _{j=i+1}^p\mathrm{Var}(f_{i,j}(X_i,X_j))+\cdots +\mathrm{Var}(f_{1,2,\dots ,p}(\mathbf{X})).$$

$$S_u=\frac{\mathrm{Var}(f_u(\mathbf{X}))}{\mathrm{Var}(f(\mathbf{X}))},$$

$$T_k=\sum _{k\in u}{S}_u,k=1,2,\dots ,p.$$

$$S_k=\frac{\mathrm{Var}(\mathbb{E}[f(\mathbf{X})|X_k])}{\mathrm{Var}(f(\mathbf{X}))}\text{and}{T}_k=\frac{\mathbb{E}[\mathrm{Var}(f(\mathbf{X})|{\mathbf{X}}_{\sim k})]}{\mathrm{Var}(f(\mathbf{X}))}.$$

$$\mathrm{Var}(\mathbb{E}[f(\mathbf{X})|X_k])=\mathbb{E}[f(\mathbf{X})(f({\mathbf{X}}_k^{\prime },{\mathbf{X}}_{\sim k})-f({\mathbf{X}}^{\prime }))],$$

$$\mathbb{E}[\mathrm{Var}(f(\mathbf{X})|{\mathbf{X}}_{\sim k})]=\frac{1}{2}\mathbb{E}[(f(\mathbf{X})-f(X_k^{\prime },{\mathbf{X}}_{\sim k}){)}^2],$$

$$\displaystyle\Big{|}\Big{|}\psi_{i}\Big{|}\Big{|}_{V_{i}}=\Big{|}\Big{|}\frac{\psi_{i}}{\phi_{i}}\Big{|}\Big{|}_{L^{\infty}(\Omega_{i})},$$

$$\displaystyle\Big{|}\Big{|}(\psi_{1},\psi_{2},\dots,\psi_{p})\Big{|}\Big{|}_{V}=\max\limits_{1\leq i\leq p}\Big{|}\Big{|}\psi_{i}\Big{|}\Big{|}_{V_{i}}.$$

$$\frac{{\varphi }_i+{\psi }_i}{1+{\int }_{{\Omega }_i}{\psi }_i(x_i)d{x}_i}$$

$$S_k(\eta )=\frac{F_k(\eta )}{H_k(\eta )}\text{and}{T}_k(\eta )=\frac{G_k(\eta )}{H_k(\eta )},$$

$$F_k(\eta )=\frac{1}{\prod _{i=1}^p{\int }_{{\Omega }_i}{\eta }_i(y)dy}{\int }_{\Omega \times \Omega }f(\mathbf{x})\left(f(x_k^{\prime },{\mathbf{x}}_{\sim k})-f({\mathbf{x}}^{\prime })\right)\prod _{i=1}^p{\eta }_i(x_i){\eta }_i(x_i^{\prime })d\mathbf{x}d{\mathbf{x}}^{\prime },$$

$$G_k(\eta )=\frac{1}{2}\frac{1}{{\int }_{{\Omega }_k}{\eta }_k(y)dy}{\int }_{\Omega \times {\Omega }_k}{\left(f(\mathbf{x})-f({\mathbf{x}}^{\prime })\right)}^2\prod _{i=1}^p{\eta }_i(x_i){\eta }_k(x_k^{\prime })d\mathbf{x}d{\mathbf{x}}_k^{\prime },$$

$$H_k(\eta )={\int }_{\Omega }f(\mathbf{x}{)}^2\prod _{i=1}^p{\eta }_i(x_i)d\mathbf{x}-\frac{1}{\prod _{i=1}^p{\int }_{{\Omega }_i}{\eta }_i(y)dy}{\left({\int }_{\Omega }f(\mathbf{x})\prod _{i=1}^p{\eta }_i(x_i)d\mathbf{x}\right)}^2.$$

$$\frac{\prod _{i=1}^p{\eta }_i}{\prod _{i=1}^p{\int }_{{\Omega }_i}{\eta }_i(y)dy}.$$

$$\mathcal{D}{S}_k(\varphi )\psi =\frac{\mathcal{D}{F}_k(\varphi )\psi }{H_k(\varphi )}-S_k(\varphi )\frac{\mathcal{D}{H}_k(\varphi )\psi }{H_k(\varphi )},\mathcal{D}{T}_k(\varphi )\psi =\frac{\mathcal{D}{G}_k(\varphi )\psi }{H_k(\varphi )}-T_k(\varphi )\frac{\mathcal{D}{H}_k(\varphi )\psi }{H_k(\varphi )},$$

$$\mathcal{D}{F}_k(\varphi )\psi ={\int }_{\Omega \times \Omega }f(\mathbf{x})(f(x_k^{\prime },{\mathbf{x}}_{\sim k})-f({\mathbf{x}}^{\prime }))\left(\sum _{i=1}^p\frac{{\psi }_i(x_i)}{{\varphi }_i(x_i)}+\frac{{\psi }_i(x_i^{\prime })}{{\varphi }_i(x_i^{\prime })}\right)\prod _{i=1}^p{\varphi }_i(x_i){\varphi }_i(x_i^{\prime })d\mathbf{x}d{\mathbf{x}}^{\prime }$$

$$-\left({\int }_{\Omega }\left(\sum _{i=1}^p\frac{{\psi }_i(x_i)}{{\varphi }_i(x_i)}\right)\prod _{i=1}^p{\varphi }_i(x_i)d\mathbf{x}\right){\int }_{\Omega \times \Omega }f(\mathbf{x})(f(x_k^{\prime },{\mathbf{x}}_{\sim k})-f({\mathbf{x}}^{\prime }))\prod _{i=1}^p{\varphi }_i(x_i){\varphi }_i(x_i^{\prime })d\mathbf{x}d{\mathbf{x}}^{\prime },$$

$$\mathcal{D}{G}_k(\varphi )\psi =$$

$$-$$

$$\mathcal{D}{H}_k(\varphi )\psi =$$

$$+\left({\int }_{\Omega }\left(\sum _{i=1}^p\frac{{\psi }_i(x_i)}{{\varphi }_i(x_i)}\right)\prod _{i=1}^p{\varphi }_i(x_i)d\mathbf{x}\right){\left({\int }_{\Omega }f(\mathbf{x})\prod _{i=1}^p{\varphi }_i(x_i)d\mathbf{x}\right)}^2$$

$$-2\left({\int }_{\Omega }f(\mathbf{x})\prod _{i=1}^p{\varphi }_i(x_i)d\mathbf{x}\right)\left({\int }_{\Omega }f(\mathbf{x})\left(\sum _{i=1}^p\frac{{\psi }_i(x_i)}{{\varphi }_i(x_i)}\right)\prod _{i=1}^p{\varphi }_i(x_i)d\mathbf{x}\right).$$

$$({\varphi }_1,{\varphi }_2,\dots ,{\varphi }_p)\mapsto {\int }_{\Omega }\prod _{i=1}^p{\varphi }_i(x_i)d\mathbf{x},$$

$${\int }_{\Omega }\sum _{j=1}^p\frac{{\psi }_j(x_j)}{{\varphi }_j(x_j)}\prod _{i=1}^p{\varphi }_i(x_i)d\mathbf{x}.$$

$${\int }_{\Omega }f(\mathbf{x})\prod _{i=1}^p{\varphi }_i(x_i)d\mathbf{x}\text{and}{\int }_{\Omega }f(\mathbf{x}{)}^2\prod _{i=1}^p{\varphi }_i(x_i)d\mathbf{x}.$$

$${\psi }_i^j(\mathbf{X})={\chi }_{R_i^j}(X_i)e_i,j=1,2,\dots ,M_i,i=1,2,\dots ,p$$

$$V_M=\text{span}\{{\psi }_i^j|j=1,2,\dots ,M_i,i=1,2,\dots ,p\}\subset V,$$

Full text

{AMS} 65C60, 62E17

1 Introduction

Uncertainties are ubiquitous in mathematical models of complex processes. Understanding how each source of input uncertainty contributes to the model output uncertainty is critical to the development and use of such models. Global sensitivity analysis (GSA) addresses this challenge by considering the inputs of the mathematical model as random variables (or implicitly assuming they are uniformly distributed) and computing a measure of their relative importance in determining the model output. A common tool for this purpose is the Sobol’ indices [19, 20, 17]. In their classical framework, one assumes that the model inputs ${\bf X}=(X_{1},X_{2},\dots,X_{p})$ are independent random variables and $f({\bf X})$ is some function mapping ${\bf X}$ to a scalar valued quantity of interest (QoI), a function of the model output. The Sobol’ indices apportion how much of the variance of $f({\bf X})$ is contributed by each input.

Computing the Sobol’ indices requires the user to specify a probability distribution for ${\bf X}$, draw samples from this distribution, and evaluate $f$ at these samples. However, the distribution of ${\bf X}$ is typically estimated from limited data and/or knowledge and hence is uncertain itself. This article studies the robustness of the Sobol’ indices to perturbations in the marginal distributions of ${\bf X}$, i.e. the distributions of $X_{i}$, $i=1,2,\dots,p$. This study is motivated by problems where the only knowledge of ${\bf X}$ is a uncertain nominal estimate, for instance a least squares estimate with noisy data, and possibly loose upper and lower bounds. This is common in many applications. For such problems, the user often assumes that each $X_{i}$ follows a symmetric distribution centered at the nominal estimate with some variance, usually chosen with limited knowledge. An example is centering a uniform distribution at the nominal estimate. The approach proposed in this article may be applied by (i) giving each $X_{i}$ a distribution with support on a wide interval with greater probability near the nominal estimate, (ii) computing the Sobol’ indices using these distributions, (iii) and assessing the Sobol’ indices robustness with respect to changes in the distributions.

In the authors previous work [10], the robustness of Sobol’ indices was studied by taking perturbations of the joint probability density function (PDF) of ${\bf X}$. This approach was useful for studying robustness with respect to statistical dependencies in ${\bf X}$. This article focuses on cases where the inputs are independent and the uncertainty is exclusively in their marginal distributions. We modify the framework introduced in [10] to develop a novel method to study the robustness of Sobol’ indices with respect to changes in the marginal distributions. It is not a generalization, specialization, nor improvement of [10], but rather a novel method which leverages ideas from [10].

The robustness of Sobol’ indices with respect to changes in the marginal distributions of ${\bf X}$ has been considered in the life cycle analysis literature [14] and ecology literature [16]. “Robust sensitivity indicators” are defined in [6] to account for uncertainty in the marginal distributions of ${\bf X}$. A generalization of the analysis of variance decomposition is considered in [2] to facilitate studying the dependence of Sobol’ indices on the distribution of ${\bf X}$. Other work such as [7, 8, 4] explore the robustness of Sobol’ indices to distributional uncertainty and [15, 11, 3, 1] explore more general questions of robustness for other computed quantities. In all of these cases, the user is required to either specify the parametric forms of distributions a priori, compute additional evaluations of $f$ beyond those needed to estimate Sobol’ indices, or both. The approach proposed in this article requires neither.

Section 2 reviews the Sobol’ indices and sets the notation for the article. We present the mathematical foundations of our proposed marginal PDF perturbation method in Section 3 and an algorithmic summary is given in Section 4. User guidelines are considered in Section 5. Numerical results are presented in Section 6. Concluding remarks are made in Section 7.

2 Review of Sobol’ indices

Let $X_{i}$ be a random variable supported on $\Omega_{i}\subset{\mathbb{R}}$, $i=1,2,\dots,p$, $\Omega=\Omega_{1}\times\Omega_{2}\times\Omega_{p}$, and $f:\Omega\to{\mathbb{R}}$. Assume that the inputs ${\bf X}=(X_{1},X_{2},\dots,X_{p})$ are independent and $f({\bf X})$ is square integrable. For $u=\{i_{1},i_{2},\dots,i_{k}\}\subset\{1,2,\dots,p\}$, let ${\bf X}_{u}=(X_{i_{1}},X_{i_{2}},\dots,X_{i_{k}})$ denote the corresponding subset of input variables. Then $f({\bf X})$ admits the following ANOVA (analysis of variance) decomposition,

[TABLE]

where

[TABLE]

The independence of the inputs $(X_{1},X_{2},\dots,X_{p})$ implies that

[TABLE]

Computing the variance of both sides of (1) and applying (2) yields

[TABLE]

Hence the variance of $f({\bf X})$ may be decomposed into contributions from each subset of input variables. The Sobol’ index for the variables ${\bf X}_{u}$ is defined as

[TABLE]

and may be interpreted as the relative contribution of ${\bf X}_{u}$ to the variance of $f({\bf X})$. In practice, one frequently computes the first order Sobol’ indices $\{S_{k}\}_{k=1}^{p}$, and the total Sobol’ indices $\{T_{k}\}_{k=1}^{p}$ defined by

[TABLE]

It easily observed that $0\leq S_{k}\leq T_{k}\leq 1$. The first order Sobol’ index $S_{k}$ may be interpreted as the relative contribution of variable $X_{k}$ to the variance of $f({\bf X})$ by itself and the total Sobol’ index $T_{k}$ may be interpreted as the relative contribution of $X_{k}$ to the variance of $f({\bf X})$ by itself and through interactions with other variables. The total Sobol’ index may also be defined for a subset of variables ${\bf X}_{u}$ in an analogous fashion. This article will focus on the indices $\{S_{k},T_{k}\}_{k=1}^{p}$; though the results are easily extended to general indices $S_{u}$ or $T_{u}$, $u\subset\{1,2,\dots,p\}$.

There are a plurality of ways to estimate the first order and total Sobol’ indices [17]. We focus on the Monte Carlo integration scheme presented in [17]. The indices $S_{k}$ and $T_{k}$ may be written probabilistically as

[TABLE]

The properties of conditional expectation and Fubini’s theorem yield

[TABLE]

and

[TABLE]

where ${\bf X}^{\prime}$ is an independent copy of ${\bf X}$.

The indices $\{S_{k},T_{k}\}_{k=1}^{p}$ may be estimated by Monte Carlo integration using $(p+2)N$ evaluations of $f$, where $N$ is the number of Monte Carlo samples. Specifically, the samples are generated by constructing 2 matrices of size $N\times p$, call them $A$ and $B$, where the rows correspond to independent samples of ${\bf X}$. Then $p$ additional matrices of size $N\times p$, denote them as $C_{k}$, $k=1,2,\dots,p$, are generated by setting the $i^{th}$ column of $C_{k}$ equal to the $i^{th}$ column of $A$, $i\neq k$, and the $k^{th}$ column of $C_{k}$ equal to the $k^{th}$ column of $B$. The set of indices $\{S_{k},T_{k}\}_{k=1}^{p}$ may be estimated using evaluations of $f$ at the rows of $A$, $B$, and $C_{k}$, $k=1,2,\dots,p$. This scheme is presented in Lines 1-6 of Algorithm 1.

The reader is directed to [19, 20, 12, 9, 17, 21, 18, 13, 22], and references therein, for additional discussion about the Sobol’ indices.

3 Marginal PDF Perturbations

We study the robustness of the Sobol’ indices to distributional uncertainty by computing a derivative of the Sobol’ index with respect to the PDFs of $X_{1},X_{2},\dots,X_{p}$. To this end we make the following assumptions:

$\Omega_{i}$ is a compact interval, $i=1,2,\dots,p$, 2. 2.

$X_{i}$ admits a PDF $\phi_{i}$, $i=1,2,\dots,p$, 3. 3.

$\phi_{i}$ is continuous on $\Omega_{i}$, $i=1,2,\dots,p$, 4. 4.

$\phi_{i}(x_{i})>0$ $\forall x_{i}\in\Omega_{i}$, $i=1,2,\dots,p$, 5. 5.

$f$ is bounded on $\Omega$.

Without loss of generality, under the assumptions above, assume that $\Omega=[0,1]^{p}$. Distributions with unbounded support may be truncated to satisfy these assumptions, hence our method is applicable on a large class of problems.

In order to perturb the PDFs $\phi_{1},\phi_{2},\dots,\phi_{p}$ while preserving the non-negativity and unit mass properties of PDFs, we define the Banach spaces $V_{i}$, $i=1,2,\dots,p$, as all bounded functions on $\Omega_{i}$ with norm

[TABLE]

where the $L^{\infty}(\Omega_{i})$ norm is the traditional supremum norm. The product Banach space $V=V_{1}\times V_{2}\times\cdots\times V_{p}$ is equipped with norm

[TABLE]

Then for $\psi=(\psi_{1},\psi_{2},\dots,\psi_{p})\in V$ with $||\psi-\phi||_{V}\leq 1$, $\phi=(\phi_{1},\phi_{2},\dots,\phi_{p})$, the functions

[TABLE]

are PDFs supported on $\Omega_{i}$, $i=1,2,\dots,p$, i.e. are non-negative and integrate to 1 over $\Omega_{i}$.

For $\eta=(\eta_{1},\eta_{2},\dots,\eta_{p})\in V$, define the first order and total Sobol’ indices as the operators $S_{k},T_{k}:V\to{\mathbb{R}}$, $k=1,2,\dots,p$,

[TABLE]

where $F_{k},G_{k},H_{k}:V\to{\mathbb{R}}$ are the operators

[TABLE]

[TABLE]

[TABLE]

For $\eta\in V$ with $||\eta-\phi||_{V}\leq 1$, $S_{k}(\eta)$ and $T_{k}(\eta)$, $k=1,2,\dots,p$, are the first order and total Sobol’ indices computed when ${\bf X}$ has the joint PDF

[TABLE]

Theorem 3.1 gives the Fréchet derivative of the Sobol’ index with respect to the PDFs.

Theorem 3.1**.**

*The operators $S_{k}$ and $T_{k}$, $k=1,2,\dots,p$, are Fréchet differentiable at the nominal PDF $\phi=(\phi_{1},\phi_{2},\dots,\phi_{p})$ and the Fréchet derivatives $\mathscr{D}S_{k}(\phi):V\to{\mathbb{R}}$ and $\mathscr{D}T_{k}(\phi):V\to{\mathbb{R}}$, $k=1,2,\dots,p$, are the bounded linear operators given by *

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

Proof 3.2**.**

One may easily observe that $H_{k}(\eta)>0$ in a neighborhood of $\phi$ (assuming $f({\bf x})$ is not constant). It is sufficient to compute the Fréchet derivatives of $F_{k},G_{k}$, and $H_{k}$ as the Fréchet derivatives of $S_{k}$ and $T_{k}$ will follow from the quotient rule for differentiation. The Fréchet derivative of the operator

[TABLE]

mapping from $V$ to ${\mathbb{R}}$, acting on $\psi=(\psi_{1},\psi_{2},\dots,\psi_{p})\in V$, is easily shown to be

[TABLE]

Similar results hold for the Fréchet derivatives of

[TABLE]

*The Fréchet derivative of $H_{k}$ follows from the sum/difference, product, and chain rules for differentiation. Similar arguments may be used to compute Fréchet derivatives of $F_{k}$ and $G_{k}$. *

Theorem 3.1 can be generalized to Sobol’ indices $S_{u}$ or $T_{u}$ for any $u\subset\{1,2,\dots,p\}$; we omit this generalization for simplicity. If the Sobol’ indices are estimated with the Monte Carlo integration scheme presented in Section 2, then the action of the Fréchet derivatives on any $\psi\in V$ may be estimated using the existing evaluations of $f$. Hence the additional computation needed to estimate the Fréchet derivative is negligible.

We would like determine an optimal perturbation of the marginal PDFs which yields the greatest change in the Sobol’ indices. For notational simplicity, let ${\mathcal{I}}$ denote a generic Sobol’ index, which may be $S_{k}$ or $T_{k}$ for any $k=1,2,\dots,p$. The locally optimal perturbation direction is given by the $\psi\in V$, $||\psi||\leq 1$, which maximizes the absolute value of the Fréchet derivative of ${\mathcal{I}}$. To estimate it, we discretize $V$ by partitioning each interval $\Omega_{i}=[0,1]$ into $M_{i}$ disjoint subintervals, $\{R_{i}^{j}\}_{j=1}^{M_{i}}$, i.e. $\Omega_{i}=\cup_{j=1}^{M_{i}}R_{i}^{j}$, $i=1,2,\dots,p$, and define basis functions $\psi_{i}^{j}:\Omega\to\{0,1\}^{p}$ as

[TABLE]

where $\chi$ is the indicator function of a set and $e_{i}$ is the $i^{th}$ canonical unit vector in ${\mathbb{R}}^{p}$. Then

[TABLE]

is a subspace of $V$ on which we will perform our robustness analysis. Discretizing $V$ with piecewise constant functions may be viewed as approximating all bounded functions by discretizing them on a coarse grid which is constrained by the existing samples and evaluations of $f$.

Then the locally optimal perturbation is solution of the optimization problem

[TABLE]

Exploiting properties of the basis we chose for $V_{M}$, (4) may be solved in closed form yielding the optimal perturbation direction

[TABLE]

where

[TABLE]

Because of the norm on $V$, this formulation gives equal weight to each marginal distribution. In actuality, we would like to focus our effort on taking perturbations on the “most important” marginals. To this end, we define the weights

[TABLE]

and define the perturbed PDF as

[TABLE]

where $\delta\in[-\overline{\delta},\overline{\delta}]$, $\overline{\delta}=\left(\max_{i=1,2,\dots,p}\delta_{i}\right)^{-1}$; the determination of $\delta$ will be discussed in Section 4. The perturbed PDF depends upon the partitions $\{R_{i}^{j}\}_{j=1}^{M_{i}}$, $i=1,2,\dots,p$, and perturbation size $\delta$, which depends on user specified parameters, see Section 4.

We will refer to the distribution of ${\bf X}$ as the nominal distribution and the distribution with PDF (5) as the perturbed distribution. Likewise, we refer to the Sobol’ indices computed using the nominal distribution as nominal Sobol’ indices. We compute perturbed Sobol’ indices, where the inputs follow the perturbed distribution, by reweighing the evaluations of $f$ used to estimate the nominal Sobol’ indices. Hence the computational effort to compute perturbed Sobol’ indices is negligible.

4 Algorithm

Evaluating the perturbed PDF requires estimation of the Sobol’ index Fréchet derivative to compute $\psi^{\star}$, and the determination of $\delta$. As previously discussed, the Sobol’ index Fréchet derivative is easily estimated using existing evaluations of $f$. The scalar $\delta$ may be determined by finding the largest and smallest possible values in $[-\overline{\delta},\overline{\delta}]$ for which the perturbed index estimator’s sample standard deviation is sufficiently small. Specifically, a threshold should be defined by the user and a value of $\delta$ is accepted if the ratio of sample standard deviations for the perturbed indices and nominal indices is below the threshold.

Algorithm 1 provides an overview of the Sobol’ index computation with marginal distribution robustness post processing. The algorithm requires four inputs:

1. $\bullet$

$n$, the number of Monte Carlo samples 2. $\bullet$

$\{R_{i}^{j}\}_{j=1}^{M_{i}}$, a partition of $\Omega_{i}$, $i=1,2,\dots,p$ 3. $\bullet$

$\tau$, the threshold for admissible increases in the Sobol’ index estimator’s sample standard deviation 4. $\bullet$

$r$, an integer used to compute admissible value of $\delta$

The determination of $n$ is problem dependent and must be done by the user to compute Sobol’ indices.

The partitions $\{R_{i}^{j}\}_{j=1}^{M_{i}}$, $i=1,2,\dots,p$, may be a chosen in many ways. A general partitioning strategy is to specific a number of subintervals $M$ and partition each $\Omega_{i}$ into $M$ subintervals according to the quantiles of the marginal distribution of $X_{i}$. Alternatively, the user may use problem specific information to customize the partition, for instance, taking a finer partition in regions of greater uncertainty or adapting $M_{i}$ based on the uncertainty in each marginal.

The threshold $\tau$ bounds the increase in estimation error introduced by reweighing samples to compute perturbed Sobol’ indices. If the original Sobol’ index estimation has small estimator standard deviation, then $\tau=1.5$ is considered a reasonable threshold to permit nontrivial perturbations while preserving reasonable accuracy of the estimator.

To determine admissible values of $\delta$, we consider a function $\Delta:[-\overline{\delta},\overline{\delta}]\to{\mathbb{R}}$ which inputs a perturbation size $\delta$ and returns the maximum (over all $p$ indices) ratio of sample standard deviations of the perturbed and nominal Sobol’ index estimators. We evaluate $\Delta$ by subsampling the data and computing a collection of estimators and their standard deviations. Finding admissible values of $\delta$ is equivalent to finding solutions of $\Delta(\delta)\leq\tau$. Intuition (and direct calculations) indicates that $\Delta$ is approximately a quadratic function centered at $\delta=0$. Admissible values of $\delta$ are determined by evaluating $\Delta$ at $r$ equally spaced points in $[-\overline{\delta},\overline{\delta}]$ and accepting those for which $\Delta\leq\tau$. Since it is discretizing an interval, $r$ does not need to be very large, $r=60$ is suggested as a reasonable choice to balance the computational cost (the for loop in Line 14 of Algorithm 1) and the necessary discretization to determine admissible values of $\delta$. The evaluations of $\Delta$ at $r$ points may be plotted to confirm that $r=60$ is sufficient.

Our choices for the inputs described above are based on numerical experimentation, intuition, and experience from the authors previous work in [10].

Algorithm 1 returns the nominal Sobol’ indices and a collection of perturbed Sobol’ indices. The user may query these perturbed Sobol’ indices and visualize measures of robustness alongside the indices themselves. In addition, samples may be drawn from the perturbed distributions. The ratios of the nominal PDF and perturbed PDFs are computed (though not explicitly written) in Algorithm 1. Using candidate samples from the nominal distribution and the ratio of PDFs provides an easy way to generate samples from the perturbed distributions with rejection sampling. Visualizing the perturbed marginal distributions provides additional insights into the regions of parameter space which cause a lack of robustness. The examples in Section 6 illustrate efficient ways to visualize the output of Algorithm 1 and make inferences with it.

5 User guidelines

The method presented in Section 3 only considers perturbations of the one dimensional marginals PDFs, we refer to this as the marginal perturbation approach. In the authors previous work [10], perturbation were taken on the joint PDF, we refer to this as the joint perturbation approach. Other forms of perturbations may be considered as well. For instance, assuming a block dependency structure motivates a form of PDF perturbations on the marginal PDFs for each block. This approach would be a hybrid between the marginal perturbation approach of this article and joint perturbation approach of [10]. In this section, we compare and contrast important features of the marginal perturbation and joint perturbation approaches to provide some principles for a user to consider when deciding which method is appropriate for a given application. Since the two approaches seek to solve different problems, the user must be cognizant of their differences when determining which to use in an application.

If the inputs are assumed to be independent for lack of better knowledge, but dependencies are suspected, then the joint perturbation approach is more appropriate. If the inputs are believed to be independent, the joint perturbation approach may give pessimistic results by introducing dependencies which are not realistic, hence the marginal perturbation approach is superior.

In the joint perturbation approach, the user must generate a partition of $\Omega$ which ensures a sufficient number of samples in each subset. If $p$ (the number of uncertain inputs) is large then the joint perturbation approach yields a partition which is coarse in each input domain $\Omega_{i}$. This problem is amplified if the nominal distribution of ${\bf X}$ is far from being uniform. In contrast, the partitioning is trivial in the marginal perturbation approach because it only requires partitioning each $\Omega_{i}$, $i=1,2,\dots,p$, separately. Hence it allows a richer class of perturbations in the marginal distributions. In addition, the user may insert a priori knowledge about each marginal distribution separately, for instance, in the partitions input to Algorithm 1.

In the marginal perturbation approach, each marginal is weighted by $\delta_{i}$, a measure of the Fréchet derivative’s magnitude. This weighting strategy manifests itself differently in the joint perturbation approach, where the input variables are prioritized by the partitioning algorithm. The optimality of the weighting scheme is unclear in both cases; however, the simplicity of the weighting scheme for the marginal perturbation approach makes it more appealing.

The user must prioritize whether there is greater uncertainty in the marginal distributions or the dependency structure. The number of uncertain inputs and characteristics of their nominal distribution influence the performance of the two approaches and must also be considered. Subsection 6.1 gives a numerical illustration comparing the two approaches.

6 Numerical results

Three examples are given below. The first example compares the marginal perturbation approach of this article with the joint perturbation approach of [10]. The second is a synthetic example serving to illustrate our marginal perturbation method. The third is an application of robustness analysis to an advection diffusion contaminant transport problem. For simplicity and conciseness, we only report results for the total Sobol’ indices.

6.1 Comparison of marginal perturbation and joint perturbation approaches

Let

[TABLE]

where $a_{i}=11-i$, $i=1,2,\dots,10$. We consider two cases, first when $X_{i}$ has a uniform distribution on the interval $[0,1]$, $i=1,2,\dots,10$, and second, when $X_{i}$ has a truncated normal distribution on the interval $[0,1]$, $i=1,2,\dots,10$. In the truncated normal case, the mean of $X_{i}$ is $0.5$ and its standard deviation is $0.73-0.04i$, $i=1,2,\dots,10$.

The marginal perturbation and joint perturbation approaches are compared by applying them to analyze the robustness of the total Sobol’ indices of (6), computed using $5,000$ Monte Carlo samples, for each of the two cases. Figure 1 displays the nominal, maximum, and minimum Total Sobol’ indices, where the maximum and minimum are determined by searching all perturbed total Sobol’ indices with $\Delta\leq\tau$. The top of the green bar and bottom of the magenta bar denote the maximum and minimum total Sobol’ index, respectively, as determined by the marginal perturbation approach; the blue region between them denotes the nominal total Sobol’ index. The top of the orange bar and bottom of the yellow bar denote the maximum and minimum total Sobol’ index, respectively, as determined by the joint perturbation approach; likewise the blue region between them denotes the nominal total Sobol’ index. The left panel corresponds to the uniform marginals case and the right panel corresponds to the truncated normal marginals case.

In both cases, the marginal perturbation approach found larger perturbation of the total Sobol’ indices than the joint perturbation approach, i.e. the marginal perturbation approach discovered a greater lack of robustness than the joint perturbation approach. This is because it was able to take a finer partition of $\Omega$ and hence find “better” perturbations.

The total Sobol’ indices in truncated normal marginals case are more robust than in the uniform marginals case. This occurs because the largest perturbations of $T_{i}$ corresponds to putting greater probability near the boundaries, i.e. for $X_{i}$ near 0 and 1. The truncated normal marginals case puts less probability near the boundaries so perturbations of this nominal PDF will not change the total Sobol’ indices as much as in the uniform marginals case. Also notice that the joint perturbation approach yields very small perturbation for the truncated marginals case. This is because there are low probability tails in several marginal distributions which hampers the partitioning algorithm. The marginal perturbation approach does not suffer this limitation.

Using the marginal perturbation approach, the maximum and minimum perturbed total Sobol’ indices are approximately symmetric about the nominal total Sobol’ indices. In the joint perturbation approach, statistical dependencies are introduced which result in a decrease of the magnitude of the total Sobol’ indices, see [9, 10]. As a result, the joint perturbation approach yields a skewness where the largest perturbations correspond to decreasing the total Sobol’ indices.

6.2 Synthetic example

To emulate a problem where there is uncertainty in the support of ${\bf X}$, let $X_{i}$ be uniformly distributed on the interval $[A_{i},B_{i}]$, where $A_{i}$ and $B_{i}$ are themselves random variables which are uniformly distributed on $[0,.1]$ and $[.9,1]$, respectively, $i=1,2,3$. Consider $f:[0,1]^{3}\to{\mathbb{R}}$ defined by

[TABLE]

It is apparent from the structure of (7) that $f$ will vary more in the $X_{2}$ direction when $X_{1}$ is close to 0 and less when $X_{1}$ is close to 1. We study the sensitivity of $f$ to $X_{2}$, and its robustness with respect to changes in the marginal distributions. Figure 2 displays the total Sobol’ indices of (7) in blue, along with the total Sobol’ indices under perturbations of the marginal distributions which maximize $T_{2}$ in red and minimize $T_{2}$ in cyan.

Figure 3 displays samples from the nominal distribution in blue (left), the perturbed distribution which maximizes $T_{2}$ in red (center), and the perturbed distribution which minimizes $T_{2}$ in cyan (right). In each column, the top row is the marginal distribution for $X_{1}$, the middle row is the marginal distribution for $X_{2}$, and the bottom row is the marginal distribution for $X_{3}$. Notice that, as expected, the distribution of $X_{1}$ is perturbed to give greater probability near 0 (1) to maximize (minimize) $T_{2}$. Further, to maximize (minimize) $T_{2}$, the distribution of $X_{2}$ is concentrated toward the boundary (interior) which has the effect of maximizing (minimizing) the variance of $X_{2}$. Conversely, the distribution of $X_{3}$ is concentrated toward the interior (boundary) to minimize (maximize) the variance of $X_{3}$, since decreasing (increasing) the variance of $X_{3}$ will increase (decrease) $T_{2}$.

In practice, the support of $X_{i}$ is frequently unknown and in many cases chosen ad hoc. Using a mixture distribution to model $X_{i}$ and coupling it with marginal perturbation robustness analysis, is a simple way to account for this uncertainty. Further, visualizing the perturbed distributions indicates whether the support will significantly influence the Sobol’ indices.

6.3 Advection diffusion example

In this example, we analyze the sensitivity of a contaminant concentration to uncertainties in its source and dynamics. Specifically, let $u(y_{1},y_{2})$ be the steady state contaminant concentration at a given point in the domain $[0,1]^{2}$ and our QoI be the average concentration in the subregion $[0.5,0.7]\times[0.5,0.7]$. To evaluate the QoI $f({\bf X})$, we solve

[TABLE]

where the velocity field and source are given by

[TABLE]

and $\mathbf{n}$ denotes the outer normal vector to the boundary. The partial differential equation is solved using Chebfun [5]. The QoI $f({\bf X})$ is then computed by integrating $u(y_{1},y_{2})$ over $[0.5,0.7]\times[0.5,0.7]$. In this example,

[TABLE]

is a vector of 9 uncertain inputs which are assumed to be independent uniform random variables with $\pm 30\%$ uncertainty about their nominal values, given in Table 1.

Figure 4 illustrates several features of this advection diffusion example. The left panel of Figure 4 shows contours of the average source (computed from 10,000 samples), with the average velocity field overlaid with red arrows. The right panel of Figure 4 displays contours of the average solution $u(y_{1},y_{2})$. In both panels, cyan lines indicate the Robin boundaries and magenta lines indicate the Dirichlet boundaries and a black box encloses the region over which we integrate to evaluate $f$.

The total Sobol’ indices of $f$ are computed via Monte Carlo integration with $n=10,000$ samples, followed by marginal distribution robustness analysis post processing. In a similar fashion to Figure 1, Figure 5 displays the nominal, maximum, and minimum total Sobol’ indices observed over all perturbations with $\Delta\leq\tau$. The top of the green bar and bottom of the magenta bar denote the maximum and minimum values, respectively; the blue region between them denotes the value of the nominal total Sobol’ index. The green and magenta bars in Figure 5 do not correspond to a particular set of Sobol’ indices for a fixed perturbation, but rather the extreme values observed over all perturbations. We observe that $\alpha_{2}$, $\nu_{1}$, and $\nu_{2}$ exhibit minimal influence on our QoI and do not change significantly when the distribution of ${\bf X}$ is perturbed. The source magnitude $\beta$ is found to be the most important parameter (under the nominal distribution), but its total Sobol’ index varies significantly when the marginal distributions are perturbed.

To further explore the influence of $\beta$, we analyze the particular perturbations which cause the greatest change in its total Sobol’ index. Figure 6 displays the nominal total Sobol’ indices in blue, the perturbed total Sobol’ indices which maximize the total Sobol’ index for $\beta$ in red, and the perturbed total Sobol’ indices which minimize the total Sobol’ index for $\beta$ in cyan. This is complemented by Figure 7 which shows the marginal distributions for each parameter using the nominal distribution in blue, the perturbed distribution which maximize the total Sobol’ index for $\beta$ in red, and the perturbed distribution which minimize the total Sobol’ index for $\beta$ in cyan. We observe that in the maximum case (red), the total Sobol’ index for $\beta$ is much larger than all other total Sobol’ indices. The perturbed distribution in this case exhibit several features:

1. $\bullet$

greater probability near the boundaries of $\beta$’s support, 2. $\bullet$

greater probability near the interior of $\gamma$’s support, 3. $\bullet$

greater probability for smaller values of $\alpha_{1}$ and $\xi_{2}$, 4. $\bullet$

greater probability for larger values of $\alpha_{2}$ and $\xi_{1}$.

The minimum case (cyan) gives a different conclusion where $\beta$ is the fifth most important parameter. The features of the perturbed distribution in this case are opposite of the maximum case (red). In particular, the source location and advection velocity become more important by giving greater probability to the source being in the top left region of the domain (smaller $y_{1}$, larger $y_{2}$), and having having greater probability in the tails of the horizontal advection speed $\alpha_{1}$.

7 Conclusion

A novel method has been presented which considers perturbations of the one dimensional marginal PDFs to measure the robustness of Sobol’ indices to distributional uncertainty. In its construction, this method assumes that the inputs are independent and the uncertainty is strictly in their one dimensional marginal distributions. This is a common occurrence in practice when, for instance, application specific information indicates that the inputs do not posses statistical dependences, but their statistical characterization is poorly known because of sparse and/or noisy data. In contrast, the authors previous work [10] provides a method to study robustness with respect to changes in the joint PDF, thus accounting for possible statistical dependencies. The source of distributional uncertainty is problem dependent and must be determined by the user. Since both methods are post processing steps after computing Sobol’ indices, the user may apply both methods and compare their results. Future work may include exportation of other PDF perturbation forms such as block dependency structures, and extensions of the PDF perturbation approach to other global sensitivity analysis methods.

Software Availability

Matlab codes are available at

https://github.com/jlhart352/Sobol_Index_Robustness

which implement the joint perturbation approach of [10] and the marginal perturbation approach presented in this article.