SOBOL' INDICES FOR PROBLEMS
DEFINED IN NON-RECTANGULAR DOMAINS
S. Kucherenko, O.V. Klymenko, N. Shah
Imperial College London, London, SW7 2AZ, UK
e‐mail: s.kucherenko@imperial.ac.uk
Abstract. A novel theoretical and numerical framework for the estimation of Sobol’ sensitivity indices for
models in which inputs are confined to a non-rectangular domain (e.g., in presence of inequality
constraints) is developed. Two numerical methods, namely the quadrature integration method which may
be very efficient for problems of low and medium dimensionality and the MC/QMC estimators based on
the acceptance-rejection sampling method are proposed for the numerical estimation of Sobol’ sensitivity
indices. Several model test functions with constraints are considered for which analytical solutions for
Sobol’ sensitivity indices were found. These solutions were used as benchmarks for verifying numerical
estimates. The method is shown to be general and efficient.
Keywords: Global sensitivity analysis, Sobol’ sensitivity indices, Dependent inputs, Constrained global
sensitivity analysis

1.

Introduction

High complexity of models in physics, chemistry, environmental studies, reliability and structural
analysis and other fields results in the increased uncertainty in model parameters and model structures.
Uncertainty in inputs is reflected in uncertainty of model outputs and predictions. Uncertainty and
sensitivity analysis has been recognized as an essential part of model applications. Global sensitivity
analysis (GSA) offers a comprehensive approach to model analysis by quantifying how the uncertainty in
model output is apportioned to the uncertainty in model inputs [1,2]. Unlike local sensitivity analysis, GSA
estimates the effect of varying a given input (or set of inputs) while all other inputs are varied as well, thus
providing a measure of interactions among variables. GSA is used to identify key parameters whose
uncertainty most affects the output. This information then can be used to rank variables, fix unessential
variables and thus decrease problem dimensionality. The variance-based method based on Sobol'
sensitivity indices has become very popular among practitioners due to its efficiency and ease of
interpretation [3,4]. Most of the developed techniques for GSA are designed under the hypothesis that
inputs are independent. However, in many cases there are dependencies among inputs, which may have
significant impact on the importance results.
There have been a number of attempts to extend GSA to models with dependent inputs. We present
a brief survey of only the most recent and easy-to-follow developments. Xu and Gertner [5] suggested to
split the contribution of an individual input to the uncertainty of the model output into two components:
the correlated contribution and the uncorrelated one. A regression-based method was used for estimating
the correlated and uncorrelated contributions of the inputs. The approach developed originally only for
linear models was later extended to nonlinear models [6]. Li et al [7] proposed a generalization of the
ANOVA-HDMR decomposition by including covariances to the model variance. They distinguished between
structural and correlative contributions of a given input. This method presents some critical issues such as
1

non-uniqueness of the functional decomposition and hence difficulties in interpreting the results. Da Veiga
et al. [8] used the estimator of the main effect index firstly proposed in [9] and further developed in [10].
They suggested a new procedure based on a heteroscedastic regression model and local polynomial
metamodel.
A novel approach for the estimation of Sobol’ sensitivity indices for models with dependent variables
using generalizations of direct Sobol’ formulas was developed in Kucherenko et al. [11]. Both the main
effect and total sensitivity indices were derived as generalizations of Sobol’ sensitivity indices. Formulas and
Monte Carlo (MC) numerical estimators similar to Sobol’ formulas were proposed. A Gaussian copula-based
approach was used for sampling from arbitrary multivariate probability distributions. The generalization
does not involve the use of surrogate models, data-fitting procedures or orthogonalization of the input
factor space.
Mara and Tarantola [12] introduced a set of sensitivity indices which relies on orthogonalisation of
correlated inputs. The computation of sensitivity indices was performed using a parametric method,
specifically the polynomial chaos expansion. Mara et al. [13] extended the development of sensitivity
indices suggested in [12] by proposing two sampling strategies for non-parametric, numerical estimation
using the Rosenblatt transformation (RT) [14]. RT, and hence values of sensitivity indices, is not unique: for
a model with n inputs there are n! possibilities corresponding to all possible permutations of the elements
of the input vector x

x1 ,..., xn . The authors considered only the RT obtained after circularly reordering

the set x1 ,..., xn , resulting in n RT transformations. The authors also established the link with the indices
proposed by Kucherenko et al. [11]. W.Hao et al. [15] suggested a detailed interpretation of indices
proposed in [11]. Taking a quadratic polynomial model with normal inputs as an illustration, they derived
explicit expressions for sensitivity indices and considered contributions to the values of sensitivity indices
from all components.
In this work we propose GSA formulations for an even wider class of problems with dependent
variables, namely for models involving inequality constraints (which naturally leads to the term
‘constrained GSA’ or cGSA). Such constraints impose structural dependences between model variables in
addition to potential correlations between them. This implies that the parameter space may no longer be
considered as an n-dimensional hypercube or hyperrectangle as considered within GSA so far, but may
assume any shapes (including disconnected domains) depending on the number and nature of constraints.
This class of problems encompasses a wide range of situations encountered in the natural sciences,
engineering, design, economics and finances where model variables are subject to certain limitations
imposed e.g. by conservation laws, geometry, costs, quality constraints etc. An important particular case
within this paradigm corresponds to imposing a minimum (maximum) threshold for the model output, i.e.,

f x1 ,..., xn

f min ,

in

g ( x1 ,..., xn )

f x1 ,..., xn

which

case

the

constraint

function

can

be

defined

as

f min and the corresponding constraint as g ( x1 ,..., xn ) 0 .

The development of efficient computational approaches for cGSA is challenging because of
potentially arbitrary shape of the feasible domain of the model variables’ variation, thus requiring the
development of special MC or quasi MC (QMC) sampling techniques and approaches for computing
sensitivity indices. This becomes especially difficult for models of high dimensionality. Another challenge
consists in analysing and interpreting model variance and sensitivity indices in such domains, since in this
case sensitivity indices will carry structural information imposed by model constraints. The interpretation of
sensitivity indices in such circumstances is crucial for the efficient design of experiments as well as for
potential model reduction by fixing or eliminating nonessential variables.

2

In this paper we have developed and tested several approaches for the numerical estimation of main
effect and total sensitivity indices in the cGSA setting. It is organized as follows: The next Section presents
formulas for the main effect and total sensitivity indices for models with dependent variables and
acceptance–rejection method which can be used to avoid sampling from conditional distributions. Section
3 considers numerical aspects of the method, namely it presents general MC estimators, acceptancerejection estimation of Sobol’ sensitivity indices using grid quadrature formulas and MC estimators. Section
4 presents the application of the proposed methodology to three test cases. Numerical and analytical
results are compared and the convergence rates of the MC and QMC methods are discussed. Finally,
conclusions are presented in the last Section.

2.
2.1.

Variance-based sensitivity indices for models with dependent inputs
Main and total effect formula for inputs defined in non-rectangular domain
defined in R n . The input vector x

Consider a model function f x1 ,..., xn

n

x1 ,..., xn consists of

real-valued random variables with a continuous joint probability distribution function (pdf) p x1 ,..., xn . It
is assumed that f x1 ,..., xn

y

xi1 ,..., xis , 1 s

has a finite variance. Consider an arbitrary subset of the variables

n and a complementary subset z

xis 1 ,..., xin , so that x ( y, z) .

The ratio

Sy

Dy [ Ez ( f ( y, z ))]

(1)

D
is known as the main effect index of the subset y . The quantity
S Ty

Ez [ Dy ( f ( y , z ))]

(2)

D
T

is known as the total effect index of the subset y . Collectively S y , S y are known as Sobol’ indices. In this
paper our notations follow closely those of Kucherenko et al. [11].

Ez ( f ( y, z )) in Eq. (1) denotes a conditional expectation with respect to z with y being fixed and

Dy ( f ( y , z )) in (2) is a conditional variance with respect to y with z being fixed. We use notations z and

z (correspondingly y and y ) to distinguish a random vector z (correspondingly y ) generated from a
joint pdf p( y, z ) and a random vector z (correspondingly y ) generated from a conditional pdf
p( y, z | y ) (correspondingly p( y , z | z ) ).
A formula for the main effect index can be explicitly written as
2

Sy
Here f0

1
D

p( y )dy
Rs

f ( y, z ) p( y, z | y)dz
Rn

f 02 .

(3)

s

E( f ( x)) , p ( y ) is the marginal pdf of subset of inputs y . This is the so-called double

loop formula (DL) [1]. It can be transformed into a different formula which in some cases can be more
efficient [11]:

3

1
D

Sy

f ( y , zˆ) p( y , zˆ | y )dzˆ

f ( y , z ) p( y , z )dy dz
R

n

R

n s

f ( y, z ) p( y, z )dydz
R

.

(4)

n

The following notation is used: z , z are two different random vectors generated from the joint pdf
p( y, z ) , a random vector ẑ is generated from a conditional pdf p( y , zˆ | y ) .
T

A formula for computing the total effect index S y (2) can be written as

S Ty

1
[ f ( y, z )
2 D Rn s

f ( y , z )]2 p( y, z ) p( y , z | z )dydydz .

(5)

Consider a more general situation when the variation of model inputs is subject to an inequality
constraint (or a finite number thereof) of the form g ( x1 ,..., xn )

0 giving rise in general to an arbitrary

(not necessarily connected) domain n Rn (Fig. 1). Despite the apparent complication of the problem of
evaluating sensitivity indices due to the introduction of constraints, the above formulas for main effect and
n

total sensitivity indices remain to be applicable in any domain
Eqs. (3)-(5) is the replacement of pdf’s with those defined in

n

Rn . The only modification required in

, which we denote in the following with a

' ( e.g. p ( y, z) ).

superscript '

For simplicity, in the following we assume that the model function f x1 ,..., xn is defined in the unit
hypercube H n (in the general case this can be achieved by applying a corresponding coordinate
transformation to the original model variables) while the permissible domain
n

n

is enclosed by H n :

Hn .

n

Fig. 1. Domain

2.2.

Hn .

Transformation from conditional to joint and marginal distributions

We assume that there is a known procedure for the generation of vectors ( y, z ) from a joint pdf

p ( y, z) defined in a non-rectangular domain
( y, z)

n

n

, which means in particular that p ( y, z )

0,

. Computation of S y and S Ty requires also sampling from conditional pdf’s. Consider

generation of a vector ( y , z ) for the evaluation of f ( y, z ) from a conditional probability distribution
function p ( y, z | y) . Using the Bayes formula p ( y, z | y) p ( y)
transformed as
4

p ( y, z ) , p ( y, z | y) can be

p ( y, z )
.
p ( y)

p ( y, z | y )

(6)

Using similar transformations for p ( y , zˆ | y ) and p ( y , z | z ) formulas (3)-(5) for S y and S Ty in
n

domain

can be presented as follows
2

Sy

Sy
S Ty

1
D

1
D

p ( y )dy
s

n s

f ( y, z )
p ( y, z )dz
p ( y)

f ( y , z ) p ( y , z )dy dz
n

1
2D

n s

n

(7)

f ( y , z)
p ( y , z )dz
p (y )

s

Here the lower-dimensional integrals are understood in the sense that
and

n s

F ( y, z )dz

( y,z)

n

,

f ( y, z ) p ( y, z )dydz

L1 (

F ( y, z )dz for any F ( y, z)

n

(8)

n

p ( y , z)
dydy dz .
p ( z)

f ( y , z )]2 p ( y, z )

[ f ( y, z )

f 02 ,

(9)

s

F ( y, z )dy

( y,z )

n

F ( y, z )dy

).

The usefulness of these formulas is based on the ability to sample from the joint pdf p ( y, z ) in

n

. Although there is a technique based on the sequential sampling from an inverse cumulative distribution
presented in Appendix A, it has a limited applicability. One example of an application of this technique is
given in Appendix B. In the next section we present a more general approach based on the acceptance –
rejection method.
2.3.

Acceptance – rejection method

Consider the joint pdf p( y, z ) in the absence of constraints in H n . We assume that constraining the
variables to an area

Hn \

n

n

H n implies that their joint pdf becomes p ( y, z ) , which takes zero values in

and is proportional to p( y, z ) within

n

functional form of the joint pdf p( y, z ) inside
n

p ( y, z )dydz

Hn

. Hence, imposing constraints does not affect the
n

apart from appropriate rescaling such that

p ( y, z )dydz 1 . Then the ‘constrained’ joint pdf can be defined through the

‘unconstrained’ one via the relationship:

p ( y, z )

p ( y, z ) I ( y, z )
n

p ( y, z )dydz

p ( y, z ) I ( y, z )
,
I

(10)

where

I

n

p ( y, z )dydz

Hn

p( y, z ) I ( y, z )dydz

is a scaling factor and I ( y, z ) is an indicator function of the subset

I ( y, z )

1, ( y, z )

n

0, ( y, z )

n

(11)

n

.

(12)

The latter can be represented through the Heaviside step function U ( ) applied to the constraint(s):

I ( y, z ) U ( g ( y, z )) .
A constrained marginal pdf is then defined as
5

(13)

p ( y)

n

1
I

p ( y, z )dz

Hn

s

p( y, z ) I ( y, z )dz ,

(14)

while a conditional pdf takes the form

p ( y, z )
p ( y)

p ( y, z | y )

p ( y, z ) I ( y, z )
H

n s

p( y, z ) I ( y, z )dz

.

(15)

T
y

The marginal pdf p ( z ) required for the evaluation of S can be expressed similarly to the formula (14).
n

Consider now a model function f ( y, z) defined in

. We note that it can be artificially extended to

H n by defining a continuation of f ( y, z) such that f ( y, z) 0 for ( y, z)

Hn \

n

to benefit from

simpler boundaries of the integration domain. Then the expected value and the variance of f ( y, z) are
given by

1
I

f0

n

f ( y, z ) p ( y, z )dydz

D

n

f 2 ( y, z ) p ( y, z )dydz

Hn

f 02

f ( y, z ) p( y, z ) I ( y, z )dydz ,

1
I

Hn

(16)

f 2 ( y, z ) p( y, z ) I ( y, z )dydz

f 02 ,

(17)

correspondingly.
Using the expressions above for pdf’s as well as the expected value and total variance of the
function, the main effect and total indices can be computed through their integral formulations presented
in (7)-(9). The important difference, however, is that even if the relevant pdf’s are known only within the
enveloping hypercube (i.e., for the unconstrained case), they can be used to directly compute their
constrained counterparts on the basis of an acceptance-rejection approach invoking the indicator function
of the feasible subdomain.
The DL formula (7) can be rewritten in the following simplified form
2

Sy

1
D

H
H

n s

f ( y, z ) p ( y, z )dz
dy

p ( y)

s

f 02

(18)

while equations (8)-(9) remain unchanged. Note that integration in (18) is performed over lowerdimensional projections H s and H n
of potentially complex shape.

s

of the unit hypercube H n as opposed to those of the domain

n

T
An alternative expression for the total effect index S y can be obtained using the identity

DTy

D Dz E y f ( y , z ) ,

(19)

where
2

Dz E y f ( y , z )

n s

p ( z )dz

s

f ( y , z ) p ( y, z | z )dy

f 02 .

(20)

Using the Bayes formula and transformations similar to presented above we obtain a DL-like formula
for total indices:
2

S Ty

1

1
D

Hs
Hn

s

f ( y, z ) p ( y, z )dy
p ( z)

6

dz

f 02 .

(21)

2.4.

Interpretation of indices

We start with Interpretation of the total index. In the case of independent inputs

DTy

Dy

Dyz ,

(22)
T

where Dyz is an interaction term between subsets y and z . From definition (19) it follows that D y is the
part of D which remains after deduction of Dz E y f ( y , z )

. In the case of dependent inputs

contains contribution to the total variance corresponding to the effect of subset z on

Dz E y f ( y , z )

its own and its dependence with subset y due to correlation or dependence via inequality constrains,
C

T

which we’ll call for brevity dependence contribution D yz . Hence, D y contains dependence contribution

D yzC with the negative sign that is D Ty does not contain correlated/dependence contribution. On the other
hand, the main effect index includes contribution corresponding to the effect of subset y on its own plus
C

dependence contribution D yz . This new meaning of Dy lead the authors of [6] to call S y the total
T

correlated contribution and S y - the total uncorrelated contribution.

3.

Numerical methods

In this section we present general MC estimators, acceptance-rejection estimation of Sobol’
sensitivity indices using grid quadrature formulas and MC estimators based on the acceptance-rejection
method.
3.1.

MC estimators

We consider MC estimators for the evaluation of integral expressions in (7)-(9) assuming that all pdf’s
are defined in

n

and there is a sampling procedure able to generate random vectors within this domain.

Formula (7) can be used to derive the DL MC estimator for S y for a single-variable index ( y

xi ). In this

case N points x(l ) , l 1, 2,..., N are generated from the joint probability distribution p ( y, z ) . For each
random variable y

xi , the sample set x(l ) , l 1, 2,..., N is sorted in the ascending order on the interval

[0,1] with respect to this variable and then subdivided into N y equally populated partitions (bins) each
containing N z

EZ f ( y, z ) | y Aj

N / N y points ( N y
1
Nz

N ). Within each bin we calculate the local mean value

Nz

f ( y jk , z jk ) , where j is the index of a bin containing N z points, j 1,..., N y
k 1

and y Aj is a mean value of y within this bin. We note, that y Aj is not actually computed and it is used only
for notation purposes. Finally, the variance of all conditional averages is computed to yield the following
double loop reordering (DLR) formula:

7

Sy

1 1
D Ny

Ny

j 1

1
Nz

2

Nz

f ( y jk , z jk )
k 1
Nz

1
Nz

2

Nz

f

1 1
D N

2
0

p ( y jk )

Ny

j 1

k 1

f ( y jk , z jk )
k 1
Nz

f 02 .

(23)

p ( y jk )
k 1

The subdivision into bins is done in the same way for all inputs using the same set of sample points. A
critical issue is the link between N and N y . It was suggested in [16] to use as a “rule of thumb”

Ny

N.

We note that it is practically impossible to generalise the DLR formula for more than one index,
hence there is no similarly efficient DLR formula which allows to compute total Sobol’ sensitivity indices.
The expected value and total variance are computed using the MC estimators

1
N

f0

1
N

D

N

f ( yl , zl ) .

(24)

l 1

N

f 2 ( yl , zl )

f 02 .

(25)

l 1

The MC estimator for the main effect index based on formula (8) has a form

Sy

1
DN

N

f ( yl , zl )
p ( yl )

f yl , zl
l 1

f yl , zl

,

(26)

and an estimator for the total effect index based on formula (9) can be written as:

S Ty

1
2 DN

N
l 1

1
f yl , zl
p ( zl )

f yl , zl

2

.

The usefulness of these estimators relies on the ability to sample from the joint pdf p ( y, z ) in

(27)
n

.

In the next two sub-sections we present a practical approach for such sampling based on the acceptancerejection method.

3.2.

Acceptance-rejection estimation of sensitivity indices using grid quadrature formulas

Integrals in (16)-(18) and (21) can be efficiently evaluated using grid quadrature methods for
problems of low and medium dimensionality. Application of grid quadrature methods is rather
straightforward when integration domains are hyperrectangular as in the case of the acceptance-rejection
approach presented in section 2.3. However, given that the feasible domain n can be arbitrary in shape,
efficient higher-order quadrature methods (such as those based on Gaussian or Clenshaw–Curtis
quadrature [17]) lose their advantages as the integrands are not differentiable at the boundary of the
permissible domain n (the indicator function has a jump discontinuity at the boundary of n ). Thus the
use of lower-order integration formulas such as the second-order multidimensional trapezoidal rule used in
this paper is preferable because they are less sensitive to non-smooth or discontinuous integrands.
Numerical integration can gain additional efficiency by performing preliminary bracketing of the
domain n (or finding its ‘minimum bounding box’) within H n to minimize the number of rejected points
during sampling. This is especially important in higher dimensions, when the number of grid points in each
dimension can’t be large due to the “curse of dimensionality”. The total number of grid points is N k n ,
where k is the number of points in each dimension. k includes both internal and boundary points and it
8

can’t be smaller than 4. We don’t consider the case of the different values of k for different directions due
the lack of prior information about importance of individual inputs in the general case.
The bracketing is performed by first using a uniform grid in H n to evaluate new lower ( ximin
upper (

ximax

where x

i

min
i

1 ) bounds for each variable. x

max
i

and x

ximin

inf xi :

ximax

sup xi :

are defined as
Hn

1

Hn

is the vector of all model variables but xi : x

0 ) and

i

I ( x)dx

I ( x)dx

1

0 ,

i

i

0 ,

(28)
(29)

( x1 ,..., xi 1 , xi 1 ,..., xn ) . Once the new bounds
n

n
have been determined the found hyperrectangle H min

ximin , ximax

which can be seen as a tight

i 1

enclosure of

n

is covered with a new grid, which is used to sample values of the model function f ( x ) at

points where the indicator function I ( x) is non-zero.
This approach was successfully implemented for problems of dimensionality up to 10. The results are
presented in Section 4. The use of this approach for higher dimensional models is computationally
prohibitive. This is also the main reason why the grid quadrature approach is not applicable to the modified
formulas (8)-(9) since the effective number of dimensions for integration in this case is (3n s ) and

(n s ) , respectively. Extension to higher dimensions may be possible with application of sparse grid
methods.

3.3.

MC estimators for the acceptance-rejection method
T

In this subsection we consider MC estimators of S y and S y valid when p ( y, z ) is known but the
constrained pdf p ( y, z ) is not known explicitly. We also assume that there is no explicit technique for
sampling in

n

of arbitrary shape. Firstly we derive an MC estimator for S y defined by (18). Its

approximation leads to the DLR approach discussed in subsection 3.1. Based on a sequence of points

( yl , zl ) , l 1,..., N sampled from the joint pdf p( y, z ) , the mean value f 0 and total variance D of
f ( y, z) are estimated as

f0
D

1
I N
1
I N

N

f ( yl , zl ) I ( yl , zl ) ,

(30)

l 1
N

2

f ( yl , zl )

f 0 I ( yl , zl ) ,

(31)

l 1

where

I

1
N

N

I ( yl , zl )

(32)

l 1

is the scaling factor.
Further we assume that y

xi . The integral F ( y )

Hn

s

f ( y, z ) p ( y, z )dz in the numerator of

the integrand of (18) can be approximated by the following estimator

9

1
I Nz

F ( y Aj )

Nz

f ( y jk , z jk ) I ( y jk , z jk ) ,

(33)

k 1

A
i.e., by averaging within each of N y bins, where j is the index of a bin and y j represents an average

value of y in the j -th bin. The marginal distribution function p ( y) in the denominator of (18) is
approximated within each bin by

1
I Nz

p ( y Aj )

Nz

I ( y jk , z jk ) .

(34)

k 1

Combining the above we obtain a DLR MC estimator for the main effect index defined by formula
(18) in the following form:

Sy

1 1
D Ny

Ny

j 1

F 2 ( y Aj )

f 02 .

A
j

p (y )

(35)

An MC estimator for S y based on the modified formula (8) has the form:

Sy

1
2
I DN

N

f yl , zl I ( yl , zl )
l 1

f ( yl , zl ) I ( yl , zl )
p ( yl )

f yl , zl I ( yl , zl )

.

(36)

Similarly an estimator for the total effect index (9) can be written as:

S Ty

1
2
2 I DN

N

f yl , zl I ( yl , zl )

f yl , zl I ( yl , zl )

l 1

2

1
,
p ( zl )

(37)

where the estimation of the marginal probability distribution p ( z ) requires additional sampling as
described below.
The marginal distribution p ( y) approximated by (38) is estimated by subdividing the whole set of
sample points ( yl , zl )

N
l 1

into N y bins according to the values of yl so that the probability of points

ending up in any of the bins is the same and equal to

1 N y (we recall that y

py

xi is scalar, y [0,1]

). To keep the same resolution in the estimation of p ( z ) an (n s ) -dimensional “ z -bin” can be defined
by choosing probabilities
earlier N z

Ny

pzi 1 N z for each i 1,..., n s . Owing to the “rule of thumb” introduced

N thus leading to similar probabilities defining a bin in all dimensions, whether y or

zi , i 1,..., n s . However, such a definition results in a dramatically decreasing number of points of the
original sample falling within such a bin when the model dimensionality n (and hence that of the bin,
n s ) increases. Indeed, the probability that a random point sampled from the joint pdf falls within a z bin defined in such a way is given by

pbin, z

n s
i 1

N z (n

pzi

s)

of the original sample expected to be found in such a bin is N bin, z
than one only if n s

N

n s
2

so that the number of points

N pbin, z

N

n s
1
2

, which is greater

2 or, in the most typical case of s 1 considered here, if n 3 . Therefore, to

estimate p ( zl ) with acceptable accuracy a relatively small additional number of samples of the form
( yq , zl ), q

1,..., N z is required which yields the estimator

p ( zl )

1
I Nz

Nz

I ( yq , zl ) .
q 1

10

(38)

We note that only the indicator function is evaluated at each of these points and although this
additional sampling at each zl increases the computational cost of the estimator (37), this increase may
not be detrimental if the computational cost of the evaluation of model constraints is significantly smaller
than that of the model function itself (which is often the case in practical problems). A value N z

64 was

used in the test cases reported below, which was found to be high enough to provide good convergence
rates.
T

There are two different ways of computing the set {Si , Si }, i

1,..., n . The first one is based on the

modified formulas (36), (37). The required number of function evaluations in this case is NCPU

N (n 2) .

The second one is based on using the DLR formula (35) for computation Si (which requires NCPU
function evaluations) and formula (37) for computation SiT (which requires NCPU

N

N (n 2) function

evaluations). However, when (35) is used in conjunction with the modified formula (37) for total effect
sensitivity indices, the former can benefit from using all the available NCPU

N (n 2) samples required

to compute SiT to enhance the accuracy of Si estimation. All derived MC estimators can be computed
using MC or QMC sampling.

4.

Test cases

For all test cases considered in this section we found analytical values of sensitivity indices, so that
they can be used as benchmarks for verification of numerical estimates.
4.1.

Product function in triangle domain

1.

Consider the function

f ( x1 , x 2 )

with both variables uniformly distributed in an upper triangle

p 1 ( x1 , x2 )
1

is defined by constraint g ( x1 , x2 )

2, ( x1 , x2 )
0, ( x1 , x2 )

(39)

x1 x 2
1
1

shown in Fig. 2 with joint pdf

.

(40)

1

x1 x2 1 0 .

Fig. 2. Definition of area
11

1.

Using the Bayes theorem the marginal pdf of x1 and the pdf of x2 conditional on x1 can be
computed explicitly:
1

1

p1 ( x1 )

p

1

1

( x1 , x2 ) dx2

1 x1

2 x1 ,

2 dx2

(41)

1 x1

p 1 ( x1 , x2 )
p1 1 ( x1 )

p2 1 ( x2 | x1 )

1
.
x1

(42)

Expectation and total variance for this function can also be computed explicitly:
1 1

f0

E{ f ( x1 , x2 )}

5
,
12

x1 x2 2 dx2 dx1
0 1 x1

1

1

5
12

2

D

( x1 x2 ) 2 dx2 dx1
0 1 x1

2

(43)

3
.
80

(44)

Using the definition (7) of the main effect index in a non-rectangular area and the relationship
T

between the main effect and total indices S x2

1 S x1 the sensitivity indices for the product function (39)

are evaluated as follows:

S x1

1
D

1

p1 1 ( x1 )dx1

1 S x1

f 02

f ( x1 , x2 ) p2 1 ( x2 | x1 )dx2

0

1 x1
1

80
3

S xT2

2

1

2

1

2 x1dx1
0

1

1
x1 x2 dx2
x1
x1

5
12

2

7
,
27

20
.
27

(45)

(46)

Owing to the symmetry of the function and the area

1

S x2

S x1 and S xT1

S xT2 .

Note that in the absence of the constraint, i.e. when the variables are uniformly distributed in the
unit square, the sensitivity indices are S x1
total variance are f0

4.2.

1 4 and D

S x2

3 7 and S xT1

S xT2

4 7 , while the mean value and

7 144 .

g-function in triangular domains

1

and

2

.

Consider the so-called g-function which is often used in GSA for illustration purposes [1]:
2

f

4 xi

1 ai

i 1

with parameters a1

0 and a 2

2

ai

(47)

1 (see Fig. 3) and x1 , x2 uniformly distributed in the unit square. Below

we compute its sensitivity indices when the feasible domain is restricted by a parametric linear (Fig. 4) or
nonlinear (Fig. 9) constraint.

12

Fig. 3. Plot of the g-function (47)

4.2.1. A linear constraint
In this subsection we consider a more general parametric linear constraint (illustrated in Fig. 4) than
that used in the test case presented in the previous subsection. The shape and size of the feasible domain
2

are defined by the variable angle

between the top side of the unit square and the line defined by

the linear constraint

0.

(48)

6 and

4 . While the latter value leads

g ( x1 , x2 ) 1 tan( ) x1 x2
Firstly we obtain analytical results for two values of
to a symmetrical domain (complementary to

1

:

in Fig. 2), the former does not result in any specific

simplification of the problem. The reference values of the mean, total variance and main effects are given
in Table 1. Exact solution for

6 was obtained using symbolic integration in Maple®, however the

resulting expressions are too cumbersome to report them here, so the values are given in the decimal form
to sufficient precision for error analysis (see below). As noted above, the values of total sensitivity indices
T

are readily obtained from the relationships S x2

1 S x1 , S xT1

1 S x2 which are valid in the 2D case.

Figs. 5, 6 show the variations of the expected value, total variance and the main effect and total
sensitivity indices with the angle

changing from 0 to

2 , which corresponds to the whole range from

the completely unconstrained 2D problem to the degenerate 1D case x1

0 , x2

0,1 . Analytical exact

values given in Table 1 are denoted by overlaid symbols. The results perfectly respect limiting cases given in
Table 1. The numerical results were obtained with the use of the grid quadrature method.
Table 1. Exact values of mean, total variance and main effects for selected values of

f0

D

S x1

S x2

0.7703487112

0.3214987100

6

0.9714128589

0.4483218079

4

1

4/9

93 9
ln 2
40 2
13

9 9
ln 2
20 8

Fig. 4. Domain
As the value of

2

as defined by the linear constraint (48).

increases the levels of importance of the variables swap with x1 being initially

significantly more important. However, as

x1

2 the domain

2

degenerates into the segment

0 , 0 x2 1 so that the model function clearly ceases to depend on x1 . This is reflected by the value
T

T

of the total effect index S x1 dropping to zero and a simultaneous increase of S x2 towards unity.

Fig. 5. Values of (a) f 0 and (b) D for the 2D g-function (47) versus angle

14

.

T

T

Fig. 6. Values of (a) S x1 and S x1 , and (b) S x2 and S x2 for the 2D g-function (47) versus angle

.

T
Fig. 7. Convergence of the MC estimators of Si (35) and Si (37) for the 2D g-function (47) defined in

for

2

6 obtained using (a) MC (b) QMC methods.

We also tested the performance of the MC estimators. The root-mean-square error (RMSE) obtained
using the MC estimator (35) for the main effect indices and the DLR estimator (37) for total effects indices
for a fixed value of

6 is presented in Fig. 7. To reduce the scatter in the error estimation the values

of RMSE were averaged over L = 50 independent runs:

i

1
L

L

I

*
i ,l

I0

2

I0

l 1

1
2

.

Here I i* is the numerical value of a particular estimator, I 0 is the corresponding analytical value. The RMSE
is approximated by a trend line cN k . Values of k are given on the plots. The convergence rates of Si and

SiT for the MC method (Fig. 7a) are close to theoretically expected value of 1/2. On the other hand, the
convergence rates for the QMC method (Fig. 7b) are significantly higher (not worse than 0.85 for the Si
T
and 0.7 for Si ). It should be noted that both approaches suffer from decreased convergence rates when

the area of

2

diminishes (data not shown) and the domain becomes elongated along the x 2 axis.

However, even in this case the QMC method still outperforms MC.

15

Fig. 8. Convergence of the DLR (35) and SK (36) MC estimators for Si for the 2D g-function (47) defined in
2

6 using (a) MC (b) QMC methods.

for

Finally, we compare the performance of the two presented estimators for the main effect indices:
namely, the DLR (35) and the estimator (36) for the modified formula (8) by S. Kucherenko et al [11], the
latter is denoted as SK. The results presented in Fig. 8 show significant advantages of the DLR approach
when considering the convergence rate versus the total number of sampling points N CPU .
4.2.2. A parabolic constraint
In this subsection we consider a more complex example involving a nonlinear constraint of the form:

g ( x1 , x2 )

x2

x1 1 x1

0,

(49)

which describes the part of the unit square above a parabola (Fig. 9). Depending on the value of
permissible domain is either connected (

4 ) or disconnected (

that the domain is non-convex for any value of

the

4 ) as illustrated in Fig. 9. We note

.

Fig. 9. A series of parabolic constraints defined by (49) for different values of

and the contour lines of the

g-function (47). See Fig. 3 for a 3D plot of the g-function
Variations of the function mean, total variance and main and total effect sensitivity indices reveal
highly nonlinear behaviour shown in Figs. 10 and 11. The dashed lines in these Figs indicate the critical
16

value

4 for the connectedness of the feasible domain. It is also worth noting that the limits

and

0

represent the two extreme cases corresponding either to the unconstrained situation (unit

square) or to the degenerate case when the feasible domain is the union of two disjoint segments

x1

0, 1, x2

0,1 , respectively. In the latter limit the numerical values of f 0 , D , main and total effect

sensitivity indices tend toward the same limit as those under a linear constraint when

2 (compare

with Figs. 5, 6).

Fig. 10. Values of (a) f 0 and (b) D for the 2D g-function (47) versus parameter

from parabolic

constraint (49).

T

T

Fig. 11. Values of (a) S x1 and S x1 , and (b) S x2 and S x2 for the 2D g-function (47) versus parameter

from

parabolic constraint (49).
Consider the particular case of

4 which is the smallest value of

when the feasible domain is

essentially disconnected. The exact values for the mean, total variance and main effect indices are:

f0

3 2,

D

521
1260

4
2,
35
17

3 384 2 575
S x1

S x2

,

144 2 521

3 264 2 293
.
2 144 2 521

(50)

T
Fig. 12. Convergence of the MC estimates of Si (35) and Si (37) for the 2D g-function (47) under the

4 obtained using (a) MC (b) QMC methods.

parabolic constraint (49) with

Fig. 12 shows the variation of RMSE versus the total number of sampling points N CPU for MC and
QMC methods. The latter is clearly superior to MC.
4.3.

K-function

K-function is defined as
i

n

K

( 1)
i 1

where variables x j , j

i

xj ,

(51)

j 1

1,..., n are independent uniformly distributed random variables in [0, 1]. K-

function is also used in GSA for illustration purposes (see f.e.[ 18]).
We consider four different cases for domain definitions. The first one is an unconstrained problem (

x

H n ). In the other three cases the unit hypercube is divided by a hyperplane into two parts one of which

is the permissible region for the problem variables x j , j
four-dimensional space n

1,..., n . All of these cases are considered in the

4 . The constraints are as follows:

I1 : x1 x2 1,

(52)

I 2 : x3 x4 1,

(53)

I3 : x1 x3 1.

(54)

These constraints can be represented using the following indicator functions: I1

U (1 x1 x2 ) ,

I 2 U (1 x3 x4 ) , I3 U (1 x1 x3 ) . See Fig. 13 for a schematic plot illustrating I1 constraint in the 3D
space.
18

For numerical estimates we make use of grid quadrature formulas (multidimensional trapezoidal
rule) presented in subsection 3.2. In order to assess the accuracy of numerical computations in the
unconstrained case the exact solution for the total effect indices reported in [19] was used while the
analytical solution for the main effect sensitivity indices was derived in this paper:

Si

1
2
3
2

2i 2

3
1
( 1) n
5
2

1
2
n 1

n i 2

1 1
10 3

1
2
n 3

2n

1
3
2

2n

.

(55)

For the situations involving the constraints (52)-(54) exact values of Si and SiT were obtained with
the aid of Maple® software.

Fig. 13. Schematic representation of permissible regions for the K-function (shaded area) in the 3D case.
The exact values of Si and SiT for all four cases are presented in Fig. 14. For the unconstrained case
the most influential input is x1 followed by x 2 with x3 and x 4 having equal and much less significant
contributions to the function variance.
Compared to the unconstrained case, the introduction of the first constraint defined by (52) (the
indicator function I 1 ) leads to a significant increase of the main effect indices S1 and S 2 accompanied by
a simultaneous decrease of S 3 and S 4 . It reflects the fact that when the constraint (52) is imposed an even
more substantial part of the variance is contributed by the first two terms of the K-function compared to
the unconstrained case.
The second constraint (53) (the indicator function I 2 ) introduces additional interaction between the
less important variables x3 and x 4 . Since x1 and x 2 are not affected by I 2 their sensitivity indices (both
main and total effects) remain at the same level as in the unconstrained case. However, variable x3
features more prominently than x 4 compared with the unconstrained case.
The last constraint (54) (the indicator function I 3 ) links an influential variable x1 and a significantly
less important x3 . It results in the decreased values of sensitivity indices of x1 , while the values of
sensitivity indices of x 2 are increased. On the other hand, the main effect index of x3 becomes
significantly more important than in the unconstrained case drastically outrunning x 4 . However, the
corresponding total effects S 3T and S4T turn out to be the same.
19

Fig. 14. (a) Main effect and (b) total sensitivity indices of the K-function in H 4 for the unconstrained and
constraints (52)-(54) cases.
Fig. 15 shows the performance of the grid quadrature integration method. Absolute errors of the
estimates of Si (Fig. 15,a) and SiT (Fig. 15,b) decrease in all cases at a rate ~ O ( N 1 / 2 ) . This is indeed
as expected since the error of the trapezoidal rule (regardless of the dimensionality) is of the form

~ O ( h 2 ) , where h is the grid spacing assumed to be equal in all dimensions. Since the number of
function evaluations in an

n -dimensional space N ~ (1 1 / h ) n (with equality instead of proportionality in

the unconstrained case) we obtain that in general

~ O( N

2/n

).

Fig. 15. Relative error in (a) Si and in (b) SiT evaluated using grid quadrature versus the number of
function evaluations NCPU , for the K-function under constraint (54) (defined by the indicator I 3 ).

20

Fig. 16. Convergence rates of the (a, b) DLR estimator (35) of Si and (c, d) modified Sobol’ estimator (37) of

SiT for the 4D K-function under constraint (52) (defined by the indicator I 1 ) obtained using (a, c) MC and
(b, d) QMC.
Fig. 16 illustrates the performance of the MC estimators for DLR approach for the main effect and the
modified Sobol’ formula for the total effect indices. Similarly to the previous examples, the use of QMC
sampling results in faster convergence.
Conclusions
In this work, we have proposed a novel concept of constrained GSA which adds the ability to analyse
model output variance in arbitrarily shaped n-dimensional domains. This amounts to greatly expanding the
scope of GSA by allowing model variables to be subject to inequality constraints, which is common in a
range of situations of practical importance.
The proposed formulas build upon Sobol’ sensitivity indices and their recent development for models
with dependent variables [11]. The advantage of the presented formulations is that no prior knowledge of
conditional or marginal distributions is assumed. All the required dependences are derived from the joint
pdf in the presence of constraints. It is shown that the knowledge of the joint pdf corresponding to the
unconstrained formulation is sufficient to build numerical estimators for sensitivity indices.

21

Three types of numerical estimators of the sensitivity indices have been proposed. Grid quadrature
may compete with MC estimators for low and medium dimensional models. However, its convergence rate
rapidly degrades with increasing model dimensionality. Despite its simple concept, the DLR approach
demonstrates good convergence when applied to the evaluation of main effects sensitivity indices
outperforming the MC estimator of the modified Sobol’ formulas. On the other hand, DLR is not a viable
approach for the evaluation of total effect indices for which the modified Sobol’ formula gives good results.
Further work is needed to develop a clear interpretation of the cGSA results, in particular by
decomposing the variance contribution into correlated and structural (uncorrelated) contributions.
Acknowledgments
The authors would like to thank Dr Stefano Tarantola and Dr. Cleo Kontoravdi for their contributions
and support of this work. One of the authors (O.V.K.) also gratefully acknowledges the financial support by
the EPSRC grant EP/K038648/1.

Appendix A. Sampling from joint and conditional distributions in a non-rectangular domain using
cumulative distributions
To sample from joint distribution p( x1 ,..., xn ) we use the following Theorem from [20]:
Theorem 1. Let
defined on

n

1

,...,

n

be independent random numbers. The set of random values { 1 ,...,

n

}

obtained from

F1 ( 1 )
F2 (

2

1

,

| 1)

2

........
Fn ( n | 1 ,...,
has the pdf p ( 1 ,...,

n

corresponding to { 1 ,...,

) . Here F1 ( 1 ) , F2 (
n

,

(56)

n 1

)

n

| 1 ) … Fn (

2

n

| 1 ,...,

n 1

) are cumulative distributions

}.

The same approach can be used to sample from conditional distributions.
Appendix B. Constant distribution in a triangle
Consider an upper triangle

1

1

shown in Fig. 2 in which random variables x1 , x2 have a joint pdf

p 1 ( x1, x2 ) 2 , ( x1 , x 2 )

1.

(57)

From the Bayes theorem:

p 1 ( x1 , x2 )

p 1 ( x1 ) p 1 ( x2 | x1 ) .

(58)

The marginal pdf for the first variable and the conditional pdf for the second variable can be
computed explicitly:
1

p1 1 ( x1 )

1

p 1 ( x1 , x2 ) dx2
1 x1

p2 1 ( x2 | x1 )

2 dx2

2 x1 ,

(59)

1 x1

p 1 ( x1 , x2 )
p1 1 ( x1 )
22

1
.
x1

(60)

From this we can find cumulative distributions:
x1

2 x1dx1 x12 ,

F1 ( x1 )
0

x2

F2 ( x2 | x1 )
1 x1

1
dx
x1

x2 (1 x1 )
x1

x1 x2 1
.
x1

(61)

Following (57) we use the following procedure to sample ( x1 , x2 ) :

x12

1

,

x1 x2 1
x1

2

(62)

.

After transformation

x1
x2

1

x1 (

,
2

(63)

1) 1.

'
'
To sample ( x1 , x2 ) and ( x1 , x2 ) two independent sets of ( 1 ,

2

) and ( 1' ,

'
2

) are needed. To

'
'
sample ( x1 , x2 | x1 ) which is needed to use formula (8) the following procedure is used:

x1'
x2

'
1

x1' (

,
2

(64)

1) 1.

This explicit sampling is limited to the case of simple geometries of n H n allowing computation
of cumulative distributions. We note, that the same test case was used in [13] to illustrate the Rosenblatt
transformation.

23

References
[1] A. Saltelli, S. Tarantola, F. Campolongo, M. Ratto, Sensitivity analysis in practice, Wiley, London,
2004.
[2] I. Sobol', S. Kucherenko, Global sensitivity indices for nonlinear mathematical models. Review,
Wilmott Mag 1 (2005) 56-61.
[3] I.Sobol’, Sensitivity estimates for nonlinear mathematical models, Matem. Mod. 2(1) (1990) 112118 [in Russian; translated in I.M. Sobol’, Sensitivity estimates for nonlinear mathematical models, Math.
Mod. Comp. Exp. 26 (1993) 407-414].
[4] I.Sobol’, Global sensitivity indices for nonlinear mathematical models and their Monte Carlo
estimates, Math. Comput. Simulat. 55(1-3) (2001) 271-280.
[5] C. Xu, G.Z. Gertner, Uncertainty and sensitivity analysis for models with correlated parameters,
Reliab. Eng. Syst. Safe. 93 (2008) 1563-1573.
[6] C. Xu, Decoupling correlated and uncorrelated uncertainty contributions for nonlinear models,
Appl. Math. Model. 37 (2013) 9950-9969.
[7] G. Li, H. Rabitz, P. Yelvington, O. Oluwole, F. Bacon, C. Kolb, J. Schoendorf, Global Sensitivity
Analysis for systems with independent and/or correlated inputs, J. Phys. Chem. A 114 (2010) 6022-6032.
[8] S. Da Veiga, F. Wahl, F. Gamboa, Local polynomial estimation for sensitivity analysis of models
with correlated inputs, Technometrics 51 (2009) 452-463.
[9] M. McKay, Variance-Based Methods for Assessing Uncertainty Importance in NUREG-1150
Analyses, in: LA-UR-96-2695, Los Alamos National Laboratory, NM, 1996, pp. 1–27.
[10] M. Ratto, S. Tarantola, A. Saltelli, Estimation of Importance Indicators for Correlated Inputs, in:
Proceedings of ESREL2001, Torino, Italy, 2001, pp. 157–164.
[11] S. Kucherenko, S. Tarantola, P. Annoni, Estimation of global sensitivity indices for models with
dependent variables, Comput. Phys. Commun. 183 (2012) 937-946.
[12] T.A. Mara, S. Tarantola, Variance-based sensitivity indices for models with dependent inputs,
Reliab. Eng. Syst. Safe. 107 (2012) 115-121.
[13] T.A. Mara, S. Tarantola, P. Annoni, Non-parametric methods for global sensitivity analysis of
model output with dependent inputs, Environ. Modell. Softw. 72 (2015) 173-183.
[14] M. Rosenblatt, Remarks on a multivariate transformation, Ann. Math. Stat. 23(3) (1952) 470472.
[15] W. Hao, Z. Lu, L. Li, A new interpretation and validation of variance based importance measures
for models with correlated inputs, Comput. Phys. Commun. 184 (2013) 1401-1413.
[16] E Plischke, An adaptive correlation ratio method using the cumulative sum of the reordered
output, Reliab. Eng. Syst. Safe.107 (2012) 149-156.
[17] L.N. Trefethen, Is Gauss Quadrature Better than Clenshaw–Curtis? SIAM Rev. 50 (2008) 67-87.
[18] S. Kucherenko, B. Feil, N. Shah, W. Mauntz, The identification of model effective dimensions
using global sensitivity analysis, Reliab. Eng. Syst. Safe. 96 (2011) 440-449.
[19] A. Saltelli, P. Annoni, I. Azzini, F. Campolongo, M. Ratto, S. Tarantola. Variance based sensitivity
analysis of model output. Design and estimator for the total sensitivity index, Comput. Phys. Commun. 181
(2010) 259-270.
[20] I.M. Sobol’, Numerical methods of Monte Carlo, Nauka, Moscow, 1973.

24

