Sobol's method
The method of Sobol' [21] is a variance-based technique. The model is y = f(x1, x2, ..., xk), where x1, x2, ..., xk are independent input factors and y is the model output. The joint pdf of the input factors is:
(1)
Mean and variance of y are:
(2)
(3)
If one of the input
factors j is
fixed to a generic value
j, the
resulting (conditional) variance of y is:
(4)
For the purpose of
sensitivity analysis one is interested in eliminating the dependence
upon the value j by
integrating V(y|
j=
j) over
the
probability density function of
j,
obtaining:
(5)
The dependence on j is
dropped from the left-hand side of eq. 5, due to the integration.
Subtracting Eq. 5 from Eq. 3 one obtains:
(6)
Thanks to a statistical
identity, the left hand side of Eq. 6 is also equal to V(E(y|j)),
and represents a measure of
the sensitivity of y
with respect
to factor j. If
one divides it by the unconditional variance, one obtains the so-called
first order sensitivity index . The Si’s
are nicely scaled in [0,1]. Eq. 6 is computationally impractical. In a
Monte Carlo framework, it implies a double loop: the inner one to
compute E2(y|
j=
j) ,
and the outer to compute
the integral over d
j.
For this reason the integral in 6 has been rewritten by Ishigami and
Homma [8] as:
(7)
The expedient of using the additional integration variable primed, allows us to realize that the integral in Eq. 7 is the expectation value of the function F of a set of (2k-1) factors:
(8)
The integral (7) can hence
be
computed using a single Monte Carlo loop. The following procedure was
proposed by Sobol’ 1993 and reviewed in Saltelli (2002).
Two input sample matrices M1 and M2 are generated:
(9)
where n is the sample
size used for the Monte Carlo estimate. In order to estimate the
sensitivity measure for a generic factor j,
i.e.
(10)
we need an estimate for both E(y) and Uj. The former can be either obtained from values of y computed on the sample in M1 or M2. Uj can be obtained from values of y computed on the following matrix Nj:
i.e. by:
(12)
If one thinks of matrix M1 as the
“sample“ matrix, and of M2 as the
“re-sample” matrix, then Ûj is obtained from
products of values of f
computed from the sample matrix times values of computed from
Nj, i.e. a matrix
where all factors except j are
re-sampled. In this way the computational cost associated with a full
set of first order indices Si is n(k+1). One set of n evaluations of f is needed to
compute , and sets of
evaluations of are needed for the second term in
the product (12).
The problem setting of Sobol’ was that of identifying a
subset of the k
factors that could account for most of the variance of y. Imagine that the
factors have been partitioned into a trial set u=(xi 1, xi 2, ..., xi m) and the
remaining set v=(xl 1, xl 2, ..., xl k-m). Then according
to Sobol’, the overall effect of the subset u
on the variance of the output can be estimated from:
(13)
(14)
(15)
In Equation (15), V(E(y|u) is the first order effect of the set u, while V(E(y|u,v) is the interaction term between the sets u and v. If V(y) ≅ V(E(y|v), then u is non-influent, and all factors in u can be fixed in a subsequent analysis of the model. Formula (13) shows the same expedient of the additional integration variables already described. The Monte Carlo estimate of Uv is:
(16)
i.e. to estimate the total
effect of set u
one must now re-sample the variables in the set u.
One can easily see that (12) is a particular case of (16). Error
estimates for Ûj’s are
discussed in the original reference of Sobol’.
Equation (15) is a particular case of a general variance decomposition
scheme proposed by Sobol’, whereby the total unconditional
variance can be decomposed as:
(17)
where
and so on. The development in (17) contains k terms of the first order Vi, k(k-1)/2 terms of the second order Vij and so on, till the last term of order k, for a total of 2k-1 terms. The Vij terms are the second order (or two-way) terms, analogous to the second order effects described in experimental design textbooks. The Vij terms capture that part of the effect of xi and xj that is not described by the first order terms. Sobol’s version of formula (17) is based on a decomposition of the function f itself into terms of increasing dimensionality, i.e.:
(18)
where each term is function only of the factors in its index, i.e. fi = fi(xi), fij = fij(xi,xj) and so on. Decompositions (17, 18) are unique provided that the input factors are independent and that the individual terms fi1i2...is in (18) are square integrable and have zero mean over the domain of existence.
Homma and Saltelli (1996) introduced the new estimate U-j :
As before:
(19)
where V(E(y|x-j)) is the
total contribution to the variance of y due to non-j. This implies that
the difference V(x) -
V(E(y|x-j)) is equal to the
sum of all terms in the variance decomposition (15) that include
j. For
a case with k=3
:
(20)
where e.g. S1= V(E(y|x1)) / V(y), and analogous expressions can be written for S2T,S3T. The SjT’s are called Total effect terms. The total effects are useful to identify input factors that are non-influential. If the total effect for a given input factor is negligible, then that factor can be fixed to any value within its range of uncertainty (and the dimensionality of the space of input can be reduced accordingly).
Observe here that when the
model is purely additive, ∑i=1,k Si =1, while for a given
factor j an important difference between SjT
and Sj flags an important role of interactions for that factor
on y.
A complete sensitivity analysis can be obtained by estimating all terms
in (17), but these are as many as 2k-1. This problem is called the curse of
dimensionality. The computational cost of estimating all effects in
(17) is in fact as high as n2k, where again n is the sample size used to estimate one individual effect ( n(2k-1) would be needed to compute all effects, and n more to compute Ê(y) ).
For this reason we customarily prefer to compute the set of
all Si plus the set of all SiT, that gives a
fairly good description of the model sensitivities. With the method of
Sobol’ this entails a computational cost
of n(2k-1) model evaluations, i.e. nk for
the first order terms, nk for the total effect terms,
plus n for Ê(y).
(21)
rather than
(22)
Equation (21) is a legitimate estimate
of the squared sample mean given the independence of the two sample
vectors. It is clear from (10) and (12) that the estimate of Sj goes more naturally to zero for a non-influential factor j when (21) is used, as can be seen from:
(23)
On the other hand the computation of the total effect indices is better achieved using (22). V(y) is computed from M1 for all indices. In conclusion, the standard computational strategy so far implemented to compute the full set of total and first order indices entails a total of n(2k+2) model evaluations, two samples being used to estimate Ê2.
PARAMETERS
The user must select which of the analyzes has to be performed. Estimation of first order and total effect indices or estimation up to a given order (custom). The total number of indexes calculated is available to the user when he has selected the Custom calculation. Number of model executions: The user must provide the sample size n
IMPORTANT: The method of Sobol’ as implemented in SimLab can only be used with noncorrelated sets of input factors.