- Email: [email protected]
A comparison of two sampling methods for global sensitivity analysis
Computer Physics Communications 183 (2012) 1061–1072
Contents lists available at SciVerse ScienceDirect
Computer Physics Communications www.elsevier.com/locate/cpc
A comparison of two sampling methods for global sensitivity analysis Stefano Tarantola, William Becker ∗ , Dirk Zeitz Joint Research Centre of the European Commission, Institute for the Protection and Security of the Citizen, Via E. Fermi 2749, 21027 Ispra, Italy
a r t i c l e
i n f o
a b s t r a c t
Article history: Received 12 May 2011 Received in revised form 5 December 2011 Accepted 21 December 2011 Available online 24 December 2011 Keywords: Latin supercube Sobol’ sequence Monte Carlo Quasi-random sequence Effective dimension
We compare the convergence properties of two different quasi-random sampling designs – Sobol’s quasiMonte Carlo, and Latin supercube sampling in variance-based global sensitivity analysis. We use the non-monotonic V-function of Sobol’ as base case-study, and compare the performance of both sampling strategies at increasing sample size and dimensionality against analytical values. The results indicate that in almost all cases investigated here, the Sobol’ design performs better. This, coupled with the fact that effective Latin supercube sampling requires a priori knowledge of the interaction properties of the function, leads us to recommend Sobol’ sampling in most practical cases. © 2012 Elsevier B.V. All rights reserved.
1. Introduction The approach used when performing a sensitivity analysis is very dependent on the problem at hand. In the simplest cases, analytical solutions may be available. In very large models, emulator approaches are often the only option – recent popular approaches include Gaussian processes [1], polynomial chaos expansions [2], and smoothing splines [3], with comparisons and summaries in [4, 5]. However, where computational cost permits, Monte Carlo (MC) and Quasi-Monte Carlo (QMC) approaches are preferred since they make the fewest assumptions about the model, and can estimate sensitivity indices of any order (unlike HDMR approaches, for instance). When implementing such (Q)MC analyses, practitioners are required to make a number of decisions. The first concerns the sampling strategy, which could be either simple random sampling, or quasi-random (low-discrepancy) sampling, including stratification approaches (which involve division of the input space into strata and random sampling within each one [6]). The former can be implemented by various available random number generators, while the latter can be realised by Sobol’, Halton, Faure or Niederreiter sequences, to name but a few – see [7] for a summary. Once the sampling strategy is decided, the second choice concerns the sampling design, which aims to prepare a specific sample to feed the model and generate outputs. Finally, the resulting sample data is used to estimate the various sensitivity
*
Corresponding author. E-mail address: [email protected] (W. Becker).
0010-4655/$ – see front matter doi:10.1016/j.cpc.2011.12.015
©
2012 Elsevier B.V. All rights reserved.
indices, which requires a choice of estimator. Of the many sensitivity estimators available, we apply here the most efficient Monte Carlo estimators [8] to calculate first and total order sensitivity indices. In this paper, we test two alternative sampling designs and investigate their performance in terms of convergence rate for global sensitivity analysis. The two designs under consideration are the randomised Sobol’ Quasi-Monte Carlo sampling routine [9,10] and Latin Supercube Sampling (LSS) [11]. LSS was developed because the efficacy of QMC has been observed to deteriorate for highdimensional problems [12,13]. The LSS design has proven to be more efficient than standard Monte Carlo [11] and the Sobol’ QMC design has proven to be more efficient than Latin hypercube sampling [14]. However, to date, no comparison between LSS and Sobol’ QMC designs has been made in the context of sensitivity analysis. The test model used in this performance assessment is the V-function of Sobol’ (also sometimes known as the G-function), which is a base case study for which analytical values of first and total order sensitivity indices are available. Different settings of the V-function are tested, varying the number of input variables and their relative importance. The paper is organised as follows. In Section 2 we give an overview of variance-based sensitivity analysis, introducing the sensitivity estimators that will be used in the rest of the paper. Section 3 describes the sampling designs used, while Section 4 contains a description of the test models and the performance measurements used. We present our results in Section 5 and we conclude in Section 6.
1062
S. Tarantola et al. / Computer Physics Communications 183 (2012) 1061–1072
2. Variance-based sensitivity analysis 2.1. Background Global sensitivity analysis is the study of how the uncertainty in the output of a computational model can be decomposed according to uncertainty in individual inputs. Recently, variance-based methods have become very popular among practitioners of global sensitivity analysis [15]. Variance-based methods have a long history in sensitivity analysis, starting with a Fourier implementation in the seventies [16], having a milestone in the work of Sobol’ [17], followed by the introduction of total sensitivity indices by Homma and Saltelli [18]. Variance-based methods have a number of interesting features such as (see [19]):
• Model independence: the sensitivity analysis can be performed for any kind of model;
• The capacity to capture the influence of the full range of variation of each input variable;
• The quantification of overall interaction effects between input variables via total sensitivity indices; • Sensitivity analysis can be performed on groups of input variables. In this section we offer an overview of variance decomposition methods (full details can be found in [17]). Given a squareintegrable function Y = f ( X 1 , X 2 , . . . , X k ) defined over Ω , the kdimensional unit hypercube:
Ω = ( X | 0 x i 1; i = 1, . . . , k )
Vi +
i
i
j >i
V i j = V X i X j E X ∼i j ( Y | X i , X j ) − V X i E X ∼i ( Y | X i )
− V X j E X∼ j ( V | X j )
(2)
and so on for higher-order terms. If we divide both sides of (1) by V (Y ), we obtain:
Si +
i
i
N 1
N
(i )
f (B) j f A B
j =1
j
− f ( A) j
which is a slight modification of the estimator proposed in [20], more suited to Sobol’ sequences. The same matrices may be used to estimate S T i . Here we use the estimator provided by Jansen [21, 22] which expresses the numerator in Eq. (3):
E X ∼i V X i ( Y | X ∼ i ) ≈
N 1
2N
j =1
(i )
f ( A) j − f A B
2 j
.
S i j + · · · + S 12...k = 1,
while the outputs f ( A B ) j are excluded from the computation as they are not independent.
A number of quasi-random sampling procedures are available in the literature. In the following, we adopt the randomised Sobol’ quasi-random sequences for both sample designs, since the LSS design requires a quasi-random number generator. We make here a distinction between “sampling strategy”, which refers to the practice of generating uniformly-distributed sample points across Ω , and the “sample design”, which we regard as the assembly of the matrices for use in the estimation of sensitivity indices. 3.1. The sampling strategy The sampling strategy refers to the selection of an appropriate random number generator. In this paper, Sobol’ quasi-random
j >i
E X ∼i ( V X i (Y | X ∼i )) V (Y )
=1−
(5)
Thus, all that is needed to compute both sets of S i and S T i for the k inputs using formulas (4) and (5) is the triplet of matrices A, B (i ) and A B . The convergence rate and bias of these estimators is of course dependent on the ability to generate a uniform sample – this is addressed in the following section. Note that the total output variance V (Y ) is estimated using only the outputs f ( A ) j and f ( B ) j (since they are independent),
where the S i , S i j , etc., are the widely-known sensitivity indices for first and higher orders. Since this variance decomposition has 2k − 1 terms, it is typical to compute only two sets of k indices: the k so-called ‘first-order’ effects ( S i ) and the k so-called ‘total effects’ S T i , which are defined as:
STi =
(4)
(1)
where
3. Sampling
V i j + · · · + V 12...k
V i = V X i E X ∼i ( Y | X i ) ,
V X i E X ∼i ( Y | X i ) ≈
(i )
the decomposition of the model output variance V (Y ), which holds for independent input variables, is:
V (Y ) =
Let us generate a sample of N points within Ω , the input space of dimension k (we ignore for the moment the issue of the sample design). The coordinates of the N points can be summarised by a matrix A with N rows and k columns. Let us now generate a further set of N points within Ω , independently from the previous set, and summarise the coordinates in a second N × k matrix, B. (i ) We now define the following matrices { A B : i = 1, . . . , k}, where all columns are from A except the ith column (i.e. the ith input variable) which is from B. Saltelli et al. suggest the following estimator for V i [8]:
V X ∼i ( E X i (Y | X ∼i )) V (Y )
.
(3)
Details on total effects can be found in [18]. In the following investigation, we will use the estimation of both the first-order and total sensitivity indices to test the two sampling designs. 2.2. On the implementation of sensitivity estimators First-order and total sensitivity estimates are based on generated samples, subsequent model evaluations, and the use of specific estimators for Eqs. (2) and (3). Here we describe the best available practice to compute S i and S T i .
Fig. 1. The Sobol’ V -function for various coefficient values.
S. Tarantola et al. / Computer Physics Communications 183 (2012) 1061–1072
1063
Fig. 2. Convergence (top) and i-wise errors at N = 8192 (bottom) for A1-1, k = 10, approximate convergence rates in parentheses.
sequences have been chosen. Quasi-random sequences are also known as “low-discrepancy” sequences, having the property of near-uniformity over Ω (strictly speaking a low-discrepancy sequence is defined by a star discrepancy below a certain measure – see [6] for further information). Unlike random numbers, quasirandom points are deterministic sequences that know about the position of previously sampled points and are constructed to avoid the presence of gaps and clusters (discrepancies) as much as possible. They generally outperform crude Monte Carlo sampling in the estimation of multidimensional integrals [23]. The Sobol’ quasi-random sequences satisfy, in particular, the following three main requirements [10]:
• Best uniformity of distribution as N → ∞. • Good distribution for fairly small initial sets. • A very fast computational algorithm. A full description of the Sobol’ sequence will not be presented here since it would require a lengthy discourse. However, the underpinnings of the theory can be found in [9], and an implementation of the algorithm and further practical details can also be found in [24]. For the sampling here, we use a commercial implementation of the Sobol’ sequence, SobolSeq8192 (developed by BRODA Ltd
and available as a C++ code), which can provide Sobol’ sequences up to dimension 8192. Since it is recommended for LSS designs [11], and to allow error estimation, we randomise each Sobol’ sample with a random shift, following [25]. Specifically, we generate a k-length vector, z, such that z → [0, 1]k , and for the jth point from the Sobol’ sequence X˜ j , we generate a shifted X j such that,
X j = X˜ j + z + X˜ j + z where α denotes the greatest integer less than or equal to α . This is equivalent to shifting the whole sample in a random direction, preserving points that move outside the unit hypercube. In this way, the structure of the sample is preserved, but a random element is introduced, which allows for variance estimation. 3.2. The sampling design Once the sampling strategy is decided, the second choice relates to the sampling design, the aim of which is to generate the triplet (i ) of matrices A, B, A B required to feed the model. 3.2.1. The Randomised Sobol’ QMC design In order to generate the matrices A and B, we generate a randomised Sobol’ sample of size N over the unit hypercube of di-
1064
S. Tarantola et al. / Computer Physics Communications 183 (2012) 1061–1072
Fig. 3. Convergence (top) and i-wise errors at N = 8192 (bottom) for A1-1, k = 30, approximate convergence rates in parentheses.
mension 2k. This creates an N × 2k matrix, T , that retains the equidistribution properties of the QMC method. The N × k matrix A is taken as the first k columns of T , with the remaining k columns given as B. The sample matrix X is defined by appending B be(i ) low A and subsequently the set of matrices A B for i = 1, . . . , k, below B. 3.2.2. The LSS design Caflisch, Morokoff and Owen introduced the concept of the “effective dimension” of a function [26], which is, roughly speaking, a measure of either the number of variables (truncation sense), or the orders of interactions (superposition sense) that are causing the majority of output variance. The point is that quite often a small subset of variables is responsible for most of the output variance, therefore the effective dimension of a function is often considerably smaller than k. LSS was introduced by Owen as a means of exploiting this kind of structure in high-dimensional integrals [11], by grouping influential variables into subsets and creating QMC designs in each subset. Owen explains that LSS is most effective when the variables in each subset interact strongly with one another. LSS can be viewed as a combination of Latin hypercube sampling (LHS) and Quasi-Monte Carlo sampling (QMC). In a basic LHS design, the support Ωi of each variable X i is divided into N
j
equally-sized intervals. For the jth interval, a sample X i is randomly drawn across the interval, giving a total of N points, spread somewhat uniformly over Ωi . This process is repeated for each X i , such that an N × k matrix is generated. Finally, the LHS design is created by randomly permuting the values in each column independently. This results in a design with N points that are roughly uniform over each dimension of Ω , and explores the full range of each variable. LSS extends LHS by dividing the k input dimensions into subd sets {s1 , s2 , . . . , sd }, such that r =1 |sr | = k, where |sr | denotes the cardinality of sr . For each subset sr a QMC design is constructed of N points, giving d matrices, each of size N × |sr |. Much like the LHS design, the full LSS design is constructed by randomly permuting the rows (samples) inside each sr , and concatenating them to make the full N × k matrix. In short, while the LHS randomises the order of samples for every dimension, LSS randomises over subsets of dimensions. LHS could therefore be viewed as a special case of LSS with |sr | = 1 ∀r, where sampling within subsets is achieved by random sampling inside equally-spaced bins (rather than QMC). Owen explains that in order to exploit the full potential of LSS, variables that interact strongly should be grouped into subsets together [11]. In practice this can be difficult, since variable interactions are generally unknown until a sensitivity analysis has been performed. To reflect this, in our investigation we examine first the case when
S. Tarantola et al. / Computer Physics Communications 183 (2012) 1061–1072
1065
Fig. 4. Convergence (top) and i-wise errors at N = 8192 (bottom) for A2, k = 10, approximate convergence rates in parentheses.
the ideal division of input variables is known, and second, where the division is unknown (this is explained further in Section 5). We therefore consider a simple case of LSS where there are only two subsets such that s1 + s2 = k. As before, we generate an N × 2k sample matrix using the Sobol’ sequence. This is split into two N × k matrices, A and B, each of which is altered to create the LSS design (randomising rows for each sr ), then the full sample matrix X is obtained in the same way as before. It is worth noting here that one drawback of the LSS design is that it does not allow sequential addition of new points. Extra points cannot be added to an existing sample since it requires a random permutation of coordinates over all sample points. Therefore if the sample size must be increased, a completely new sample must be designed.
to O ( N −1/2 ) (i.e. that of purely random Monte Carlo), only obtaining the higher QMC rate above a certain (high) value of N. This transition value appears to increase with dimensionality, but is also strongly affected by the nature of the integrand (the effective dimensionality) [12,13]. 4. Setting the experiment 4.1. The test function The test function employed in this work has been widely used as benchmark for sensitivity analysis (see [17]) and is defined as:
Y = V ( X 1 , X 2 , . . . , X k, a1 , a2 , . . . , ak ) =
k
vi
i =1
3.2.3. Properties of the sample designs It is worth considering briefly the properties of the designs here. The LSS design has been shown to exhibit “asymptoticallynegligible bias”, so long as randomised-QMC is used [11]. This comes at the expense of a small increase in the variance of the estimate due to the introduced randomisation, but was shown to be preferable. The Sobol’ sequence is a low-discrepancy design that has good uniform properties as N increases (discrepancy follows O ( N −1 logn−1 N ) for N = 2k ) [27]. However, it has been identified that for large dimensionality, convergence of integrals can reduce
with,
vi =
|4 X i − 2| + a i . 1 + ai
The nature of the function V is driven by the dimensionality k as well as by the values of the coefficients ai . These latter determine the relative importance of input variables (i.e., the lower ai the higher the dependence of Y on X i – see Fig. 1) and can therefore be used to generate a wide spectrum of test cases of different levels of complexity.
1066
S. Tarantola et al. / Computer Physics Communications 183 (2012) 1061–1072
Fig. 5. Convergence (top) and i-wise errors at N = 8192 (bottom) for A2, k = 30, approximate convergence rates in parentheses.
4.2. Variations of the test function The test cases used in this work can be grouped according to a classification proposed in [14] which defines three different types of functions (here labelled A, B and C) with different characteristics and different effective dimensionality: Type A includes functions with a few dominant input variables, which could also be involved in interactions (low effective dimension in the superposition and truncation sense). These are the most common functions in practice. Type B includes functions with a homogeneous effect of all input variables, which are involved in some interactions (effective dimension is roughly equal to k in the truncation sense). Type C is like Type B but the input variables are involved in strong interactions (effective dimension is similar to k in both truncation and superposition senses). In this work we present five test cases, classified according to the scheme above, that are representative of a wide range of model behaviours. For each of them we know the analytical values of the sensitivity indices (see Appendix A). We compare those exact values with the estimates obtained through the model evaluations. A1-1/A1-2: The first test case (called A1-1) uses the coefficients a1 = a2 = 0, ai >2 = 6.52; the second test case (A1-2) uses the co-
efficients ak−1 = ak = 0; ai
S i = 0.8268
i =1
from which it follows that 17.32% of the variance in the model output is due to interactions between the input variables. Note that these models have a low effective dimension, which could be exploited by the LSS design. However, there are no strong interactions. A2: The third test case, called A2, uses the coefficients a1 = a2 = · · · = a5 = 0; ai >5 = 6.52, which yield a different decomposition of the sensitivity in terms of first, second and higher-order effects. Here, the first-order effects account for 50.79% of the output variance (k = 10):
S. Tarantola et al. / Computer Physics Communications 183 (2012) 1061–1072
1067
Fig. 6. Convergence rate of Sobol’, LSS and ignorant LSS with model dimensionality for A2 function.
Fig. 7. Comparison of absolute MAES values for A2 function at increasing N and k.
k
S i = 0.5079
i =1
while interactions, at any order, account for the remaining 49.21% of the output variance. Given that the majority of interaction variance is caused by the set {x1 , . . . , x5 }, this should be an ideal case for the LSS design if the subsets are chosen appropriately. B: The fourth test case (called B) uses the coefficients ai = 6.52 for i = 1, . . . , k. The variance is therefore spread evenly over all inputs, mostly in first-order effects. Type B functions are quite rare:
models encountered in practice usually have a few dominant input variables and a larger number of less important ones. C: The final test case C uses the coefficients ai = 0 for i = 1, . . . , k. In contrast with B, the variables have important interactions, such that the sum of first-order effects for k = 10 is k
S i = 0.1989.
i =1
This means that around 80% of the variance in the output is driven by interactions. Again, this is not a common type of function, but is worth investigating for completion.
1068
S. Tarantola et al. / Computer Physics Communications 183 (2012) 1061–1072
Fig. 8. Convergence (top) and i-wise errors at N = 8192 (bottom) for B, k = 10, approximate convergence rates in parentheses.
As will be seen in the following sections, it becomes increasingly difficult to achieve a desired accuracy when estimating sensitivity indices as the sum of first-order interactions becomes smaller (i.e. in case C in particular). 4.3. Sample approach and dimensionality In order to assess the application of these designs to models of various dimensionality, we investigate the cases where k is equal to 10, 20, 30, 50 and 100. For each model, we run the experiments at increasing sample size N to check convergence of the sensitivity estimates to the true values. We set the sample size N at values 16, 32, 64, 128, 256, 512, 1024, 2048, 4096 and 8192. Finally, each experiment at a given N and k is replicated R = 50 times in order to quantify the statistical error of the sensitivity estimates, in terms of deviation between the estimates and the true values. 4.4. Measuring convergence We use four indicators to quantify the statistical error of the sensitivity estimates at given N and k – these are named AES, MAES, AEST and MAEST. Let us denote by S¯ i the analytic value of a first-order index for (r ) the ith input and by S i the corresponding estimated value of the
rth replication. We now define AES as the mean absolute error of each S i :
AESi =
R 1 (r ) S − S¯ i i R r =1
i.e., these help to investigate which model input causes most of the deviation from the analytical values. MAES is the average of the AES values over all inputs to give a single error measure,
1 1 (r ) S − S¯ i . i Rk k
MAES =
R
i =1 r =1
The definitions of AEST and MAEST are exactly equivalent to AES and MAES respectively, but measure errors based on estimates of S T i . For any test function and any k we regress both MAES and MAEST against N by a trend line c ∗ N −a (power regression) in a log–log plot. Values of a ∈ [0, 1] are given in brackets in the legend of the figures and show the rate of convergence of the sensitivity estimates for each implementation design at increasing sample size (larger a implies better convergence properties). For reference, standard Monte Carlo has a = 0.5. The coefficient c is related to the approximation error at small N, therefore small values of c mean that the implementation performs well with a small sample size.
S. Tarantola et al. / Computer Physics Communications 183 (2012) 1061–1072
1069
Fig. 9. Convergence (top) and i-wise errors at N = 8192 (bottom) for B, k = 30, approximate convergence rates in parentheses.
5. Results and analysis The following experiments compare convergence of sensitivity estimates using randomised Sobol’ against LSS designs, where the subsets are chosen with knowledge of the grouping of variables. Since, in practice, the interactions of variables are rarely known a priori, convergence is also calculated for an “ignorant” LSS sample. In this case, we assume no knowledge of the grouping of variables, and choose two subsets such that |s1 | = |s2 | = k/2. The intention is to see how LSS performs with no knowledge of groups, or if the groups are incorrectly specified. For the Type B and Type C functions, there is no “correct” grouping, since all variables are equally important. Therefore we adopt the “ignorant” LSS approach as default. 5.1. Type A1 functions Fig. 2 (top) shows the convergence of estimates with increasing N for both sample designs, in the A1-1 configuration, where k = 10. It is evident that in this case the Sobol’ sampling design is superior to the LSS implementation (with little difference between correctly-specified groups and ignorance): the values of MAES and MAEST are generally smaller and their slope at increasing N is steeper. Sobol’ has the largest values of a, however it should be
noted that at small sample sizes LSS tends to have a lower error. Moreover, AES and AEST show that the main contribution to the statistical error comes from the first two input variables. When k = 30 (Fig. 3, top) MAES values for Sobol’ are again almost always better than those of LSS. Interestingly however, the MAEST values are consistently superior for LSS sampling. This could be due to the deterioration of the Sobol’ design in certain dimensions as k becomes large. In contrast, LSS (if subsets are correctly specified) will not have this problem because each subset uses lower-dimensional Sobol’ sequences. It seems that on average, for this type of function, the Sobol’ design is the better choice. The results of the A1-2 function (where the driving variables are the last two, rather than the first two) were different to those of the A1-1 function. While the MAES curves were similar, the Sobol’ sampling outperformed LSS in estimating S T , which suggests that the ordering of variables is important. 5.2. Type A2 function The Type A2 function is similar to Type A1-1 but with five (instead of only two) important input variables, and a greater proportion of variance due to interactions. In the LSS implementation, these five model input constitute group s1 (this is therefore the
1070
S. Tarantola et al. / Computer Physics Communications 183 (2012) 1061–1072
Fig. 10. Convergence (top) and i-wise errors at N = 8192 (bottom) for C, k = 10, approximate convergence rates in parentheses.
ideal subset for the LSS design). In the case of k = 10, the results (Fig. 4) show that, while there is little difference in the estimation of first-order indices, Sobol’ sampling outperforms LSS in estimation of total effect indices. The results for k = 30 are basically identical (Fig. 5). LSS and Sobol’ designs perform very similarly in terms of MAES, but Sobol’ outperforms LSS at all sample sizes with respect to total effect indices (MAEST). Fig. 6 shows the result of measuring the convergence rate for the A2 function at increasing dimensionality. It is interesting to note here that the rate of convergence is nearly always superior for the Sobol’ design. However, besides the convergence rate, it is also important to consider the absolute error at a given sample value. To further investigate, Fig. 7 shows the difference between values of log(MAES) from the Sobol’ and LSS designs, for increasing values of k and N (essentially the pointwise difference between the convergence curves). The indication here is that Sobol’ sampling outperforms LSS at large sample sizes. For lower sample sizes (below 128) the lack of a clear pattern could be interpreted as noise since the number of sample points is too few to give meaningful estimates of sensitivity indices, particularly for higher-dimensional functions.
5.3. Type B function Examining the Type B function results in Fig. 8 (equal importance in all inputs, and little interaction), it is evident that for k = 10 the Sobol’ sampling design is superior to LSS in terms of MAES and MAEST for larger sample sizes. The convergence rate a for the Sobol’ design is very good (a ≈ 1), whereas the rate for LSS designs is around 0.50, typical of a standard Monte Carlo sampling strategy. This is to be expected – the strength of LSS lies in assembling inputs into interacting subsets, which cannot be achieved here. When k = 30 (see Fig. 9), in general, the convergence rate of the Sobol’ sampling decreases (compared to k = 10). For the estimation of S i , Sobol’ still outperforms LSS at large N, but its efficacy is decreased. However, LSS performs better in estimation of total effect indices. In the case where N = 8192, this is due to a large error in the Sobol’ estimation of S T for variable 20, which takes its MAEST value above that of LSS. Since the Type B function is symmetrical in all dimensions, this must be due to a quirk of the Sobol’ sequence. Upon further investigation, it was found that this “spike” was also present at lower sample sizes (although it appeared on different input variables), which caused the Sobol’ MAEST values to be consistently higher than those of LSS.
S. Tarantola et al. / Computer Physics Communications 183 (2012) 1061–1072
1071
Fig. 11. Convergence (top) and i-wise errors at N = 8192 (bottom) for C, k = 30, approximate convergence rates in parentheses.
5.4. Type C function Fig. 10 shows the performance of the designs on the C-type function, for the case where k = 10. The designs appear similar for MAES, although there is a considerable advantage for the Sobol’ design for MAEST. The results for k = 30 (Fig. 11) are very similar: for MAES the two designs show little difference, whereas for MAEST the Sobol’ design show a fairly clear advantage over LSS. In both designs there are error spikes over certain inputs – these are also found at other sample sizes, though not necessarily on the same variables. 6. Conclusion In this work we have presented five test cases, grouped according to the classification proposed by Kucherenko, that are representative of a wide range of model behaviours. The results here are based on the Sobol’ V-function (indeed, on particular implementations of this function), therefore it is important to stress that the conclusions here may not necessarily be extrapolated out to all functions in each group, though they can hopefully provide an indication. For A-type functions, it was found that the Sobol’ design generally performs better, particularly in the case where there is a larger subset of inputs interacting quite strongly (A2 functions). However,
the LSS design can offer some advantages, although these were observed to be restricted to particular cases such as the A1 function, and only ever with respect to the estimation of S T . A more detailed study of convergence in the A2 function revealed that for all dimensionalities investigated, the Sobol’ design was consistently superior. Furthermore, a comparison of absolute error for given dimensionality and sample size revealed that Sobol sampling almost always performed better. B-type and C-type functions were those that give equal importance to all inputs, with weak and strong interactions respectively. Since LSS is based on identifying a subset of strongly-interacting inputs, it was not expected that LSS would offer any advantage over a Sobol’ design, and indeed this was the case in general. However, one characteristic that was that the Sobol’ design tended to distribute the total estimation error more unevenly across inputs, i.e., low error for most indices, then a relatively high error for one or two inputs. This is thought to be due to a deterioration of the Sobol’ design in high dimensions. The LSS design did not generally suffer from this (despite being based on the Sobol’ design), which can be viewed as an advantage, since it would not be known which of the sensitivity indices were the accurate ones in most practical situations. In general, it is difficult to recommend LSS as a sampling strategy, unless the effective dimension of the function is known and is compatible with the strengths of LSS. It may be that in some cases,
1072
S. Tarantola et al. / Computer Physics Communications 183 (2012) 1061–1072
this information is available from some background knowledge of the function (the physical process behind a model perhaps), or from a preliminary screening test. However in general, little or nothing is known about the response of the model or function. In such cases the Sobol’ sequence provides a good general-purpose sampling approach. Finally, it should be noted that while the results here favour Sobol’ sampling over LSS, they are based on a particular test function which may not necessarily generalise to all model responses in practice. Appendix A. V -function and its analytical values The sensitivity indices for the V-function can be obtained analytically from the following,
V i = V X i E X i (Y | X i ) = VTi = Vi
1 /3 , (1 + ai )2
k (1 + V j ), j =i
V =
k
(1 + V i ) − 1
i =1
where S i = V i / V and S T i = V T i / V . Note that V here denotes variance, not to be confused with the variables of the V -function. The analytical values for the test cases are therefore: k = 30
k = 100
Type A1-1
sum( S i ) S 1 / S Tot 1 S k / S Tot k
k = 10 0.83 0.72 0.54
0.76 0.64 0.48
0.58 0.42 0.32
Type A1-2
sum( S i ) S 1 / S Tot 1 S k / S Tot k
0.83 0.54 0.72
0.76 0.48 0.64
0.58 0.32 0.42
Type A2
sum( S i ) S 1 / S Tot 1 S k / S Tot k
0.51 0.31 0.23
0.47 0.27 0.21
0.35 0.18 0.14
Type B
sum( S i ) S 1 / S Tot 1 S k / S Tot k
0.97 0.95 0.95
0.92 0.84 0.84
0.74 0.56 0.56
Type C
sum( S i ) S 1 / S Tot 1 S k / S Tot k
0.20 0.08 0.08
Functions
∼0 ∼0 ∼0
∼0 ∼0 ∼0
References [1] J.E. Oakley, A. O’Hagan, Probabilistic sensitivity analysis of complex models: a Bayesian approach, Journal of the Royal Statistical Society Series B 66 (3) (2004) 751–769. [2] T. Crestaux, I. Le Ma, O. Tre, J.-M. Martinez, Polynomial chaos expansion for sensitivity analysis, Reliability Engineering & System Safety 94 (7) (2009) 1161–1172.
[3] M. Ratto, A. Pagano, Using recursive algorithms for the efficient identification of smoothing spline ANOVA models, AStA Advances in Statistical Analysis 94 (4) (2010) 367–388. [4] R. Jin, W. Chen, T.W. Simpson, Comparative studies of metamodelling techniques under multiple modelling criteria, Structural and Multidisciplinary Optimization 23 (1) (2001) 1–13. [5] C.B. Storlie, L.P. Swiler, J.C. Helton, C.J. Sallaberry, Implementation and evaluation of nonparametric regression procedures for sensitivity analysis of computationally demanding models, Reliability Engineering & System Safety 94 (11) (2009) 1735–1763. [6] R.E. Caflisch, Monte Carlo and quasi-Monte Carlo methods, Acta Numerica 7 (1998) 1–49. [7] H. Niederreiter, Random Number Generation and Quasi-Monte Carlo Methods, Society for Industrial Mathematics, 1992. [8] 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, Computer Physics Communications 181 (2) (2010) 259–270. [9] I.M. Sobol’, On the distribution of points in a cube and the approximate evaluation of integrals, USSR Computational Mathematics and Mathematical Physics 7 (4) (1967) 86–112. [10] I.M. Sobol’, Uniformly distributed sequences with an additional uniform property, USSR Computational Mathematics and Mathematical Physics 16 (5) (1976) 236–242. [11] A.B. Owen, Latin supercube sampling for very high-dimensional simulations, ACM Transactions on Modeling and Computer Simulation 8 (1) (1998) 71–102. [12] W.J. Morokoff, R.E. Caflisch, Quasi-random sequences and their discrepancies, SIAM Journal on Scientific Computing 15 (6) (1994) 1251–1279. [13] W.J. Morokoff, R.E. Caflisch, Quasi-Monte Carlo integration, Journal of Computational Physics 122 (2) (1995) 218–230. [14] S. Kucherenko, B. Feil, N. Shah, W. Mauntz, The identification of model effective dimensions using global sensitivity analysis, Reliability Engineering & System Safety 96 (4) (2011) 440–449. [15] A. Saltelli, et al., Global Sensitivity Analysis: The Primer, John Wiley & Sons, Ltd., Chichester, 2008, x, 292 p. [16] R.I. Cukier, C.M. Fortuin, K.E. Shuler, A.G. Petschek, J.H. Schaibly, Study of the sensitivity of coupled reaction systems to uncertainties in rate coefficients. I. Theory, Journal of Chemical Physics 59 (8) (1973) 3873–3879. [17] I.M. Sobol’, Sensitivity estimates for nonlinear mathematical models, Mathematical Modeling and Computational Experiment 1 (4) (1993) 407–414. [18] T. Homma, A. Saltelli, Importance measures in global sensitivity analysis of nonlinear models, Reliability Engineering & System Safety 52 (1) (1996) 1–17. [19] A. Saltelli, S. Tarantola, F. Campolongo, M. Ratto, Sensitivity Analysis in Practice: A Guide to Assessing Scientific Models, John Wiley & Sons Inc., 2004. [20] I.M. Sobol’, S. Tarantola, D. Gatelli, S.S. Kucherenko, W. Mauntz, Estimating the approximation error when fixing unessential factors in global sensitivity analysis, Reliability Engineering & System Safety 92 (7) (2007) 957–960. [21] M.J.W. Jansen, Analysis of variance designs for model output, Computer Physics Communications 117 (1–2) (1999) 35–43. [22] M.J.W. Jansen, W.A.H. Rossing, R.A. Daamen, Monte Carlo estimation of uncertainty contributions from several independent multivariate sources, in: G. van Straaten, J. Grasman (Eds.), Predictability and Nonlinear Modelling in Natural Sciences and Economics, Kluwer Academic Publishing, Dordrecht, 1994, pp. 335–343. [23] I.M. Sobol’, S.S. Kucherenko, On global sensitivity analysis of quasi-Monte Carlo algorithms, Monte Carlo Methods and Applications 11 (1) (2005) 83–92. [24] P. Bratley, B.L. Fox, ALGORITHM 659: implementing Sobol’s quasirandom sequence generator, ACM Transactions on Mathematical Software 14 (1) (1988) 88–100. [25] R. Cranley, T.N.L. Patterson, Randomization of number theoretic methods for multiple integration, SIAM Journal on Numerical Analysis 13 (6) (1976) 904– 914. [26] R.E. Caflisch, W. Morokoff, A. Owen, Valuation of mortgage backed securities using Brownian bridges to reduce effective dimension, Journal of Computational Finance 1 (1) (1997) 27–46. [27] I.M. Sobol, On the systematic search in a hypercube, SIAM Journal on Numerical Analysis 16 (5) (1979) 790–793.
Contents lists available at SciVerse ScienceDirect
Computer Physics Communications www.elsevier.com/locate/cpc
A comparison of two sampling methods for global sensitivity analysis Stefano Tarantola, William Becker ∗ , Dirk Zeitz Joint Research Centre of the European Commission, Institute for the Protection and Security of the Citizen, Via E. Fermi 2749, 21027 Ispra, Italy
a r t i c l e
i n f o
a b s t r a c t
Article history: Received 12 May 2011 Received in revised form 5 December 2011 Accepted 21 December 2011 Available online 24 December 2011 Keywords: Latin supercube Sobol’ sequence Monte Carlo Quasi-random sequence Effective dimension
We compare the convergence properties of two different quasi-random sampling designs – Sobol’s quasiMonte Carlo, and Latin supercube sampling in variance-based global sensitivity analysis. We use the non-monotonic V-function of Sobol’ as base case-study, and compare the performance of both sampling strategies at increasing sample size and dimensionality against analytical values. The results indicate that in almost all cases investigated here, the Sobol’ design performs better. This, coupled with the fact that effective Latin supercube sampling requires a priori knowledge of the interaction properties of the function, leads us to recommend Sobol’ sampling in most practical cases. © 2012 Elsevier B.V. All rights reserved.
1. Introduction The approach used when performing a sensitivity analysis is very dependent on the problem at hand. In the simplest cases, analytical solutions may be available. In very large models, emulator approaches are often the only option – recent popular approaches include Gaussian processes [1], polynomial chaos expansions [2], and smoothing splines [3], with comparisons and summaries in [4, 5]. However, where computational cost permits, Monte Carlo (MC) and Quasi-Monte Carlo (QMC) approaches are preferred since they make the fewest assumptions about the model, and can estimate sensitivity indices of any order (unlike HDMR approaches, for instance). When implementing such (Q)MC analyses, practitioners are required to make a number of decisions. The first concerns the sampling strategy, which could be either simple random sampling, or quasi-random (low-discrepancy) sampling, including stratification approaches (which involve division of the input space into strata and random sampling within each one [6]). The former can be implemented by various available random number generators, while the latter can be realised by Sobol’, Halton, Faure or Niederreiter sequences, to name but a few – see [7] for a summary. Once the sampling strategy is decided, the second choice concerns the sampling design, which aims to prepare a specific sample to feed the model and generate outputs. Finally, the resulting sample data is used to estimate the various sensitivity
*
Corresponding author. E-mail address: [email protected] (W. Becker).
0010-4655/$ – see front matter doi:10.1016/j.cpc.2011.12.015
©
2012 Elsevier B.V. All rights reserved.
indices, which requires a choice of estimator. Of the many sensitivity estimators available, we apply here the most efficient Monte Carlo estimators [8] to calculate first and total order sensitivity indices. In this paper, we test two alternative sampling designs and investigate their performance in terms of convergence rate for global sensitivity analysis. The two designs under consideration are the randomised Sobol’ Quasi-Monte Carlo sampling routine [9,10] and Latin Supercube Sampling (LSS) [11]. LSS was developed because the efficacy of QMC has been observed to deteriorate for highdimensional problems [12,13]. The LSS design has proven to be more efficient than standard Monte Carlo [11] and the Sobol’ QMC design has proven to be more efficient than Latin hypercube sampling [14]. However, to date, no comparison between LSS and Sobol’ QMC designs has been made in the context of sensitivity analysis. The test model used in this performance assessment is the V-function of Sobol’ (also sometimes known as the G-function), which is a base case study for which analytical values of first and total order sensitivity indices are available. Different settings of the V-function are tested, varying the number of input variables and their relative importance. The paper is organised as follows. In Section 2 we give an overview of variance-based sensitivity analysis, introducing the sensitivity estimators that will be used in the rest of the paper. Section 3 describes the sampling designs used, while Section 4 contains a description of the test models and the performance measurements used. We present our results in Section 5 and we conclude in Section 6.
1062
S. Tarantola et al. / Computer Physics Communications 183 (2012) 1061–1072
2. Variance-based sensitivity analysis 2.1. Background Global sensitivity analysis is the study of how the uncertainty in the output of a computational model can be decomposed according to uncertainty in individual inputs. Recently, variance-based methods have become very popular among practitioners of global sensitivity analysis [15]. Variance-based methods have a long history in sensitivity analysis, starting with a Fourier implementation in the seventies [16], having a milestone in the work of Sobol’ [17], followed by the introduction of total sensitivity indices by Homma and Saltelli [18]. Variance-based methods have a number of interesting features such as (see [19]):
• Model independence: the sensitivity analysis can be performed for any kind of model;
• The capacity to capture the influence of the full range of variation of each input variable;
• The quantification of overall interaction effects between input variables via total sensitivity indices; • Sensitivity analysis can be performed on groups of input variables. In this section we offer an overview of variance decomposition methods (full details can be found in [17]). Given a squareintegrable function Y = f ( X 1 , X 2 , . . . , X k ) defined over Ω , the kdimensional unit hypercube:
Ω = ( X | 0 x i 1; i = 1, . . . , k )
Vi +
i
i
j >i
V i j = V X i X j E X ∼i j ( Y | X i , X j ) − V X i E X ∼i ( Y | X i )
− V X j E X∼ j ( V | X j )
(2)
and so on for higher-order terms. If we divide both sides of (1) by V (Y ), we obtain:
Si +
i
i
N 1
N
(i )
f (B) j f A B
j =1
j
− f ( A) j
which is a slight modification of the estimator proposed in [20], more suited to Sobol’ sequences. The same matrices may be used to estimate S T i . Here we use the estimator provided by Jansen [21, 22] which expresses the numerator in Eq. (3):
E X ∼i V X i ( Y | X ∼ i ) ≈
N 1
2N
j =1
(i )
f ( A) j − f A B
2 j
.
S i j + · · · + S 12...k = 1,
while the outputs f ( A B ) j are excluded from the computation as they are not independent.
A number of quasi-random sampling procedures are available in the literature. In the following, we adopt the randomised Sobol’ quasi-random sequences for both sample designs, since the LSS design requires a quasi-random number generator. We make here a distinction between “sampling strategy”, which refers to the practice of generating uniformly-distributed sample points across Ω , and the “sample design”, which we regard as the assembly of the matrices for use in the estimation of sensitivity indices. 3.1. The sampling strategy The sampling strategy refers to the selection of an appropriate random number generator. In this paper, Sobol’ quasi-random
j >i
E X ∼i ( V X i (Y | X ∼i )) V (Y )
=1−
(5)
Thus, all that is needed to compute both sets of S i and S T i for the k inputs using formulas (4) and (5) is the triplet of matrices A, B (i ) and A B . The convergence rate and bias of these estimators is of course dependent on the ability to generate a uniform sample – this is addressed in the following section. Note that the total output variance V (Y ) is estimated using only the outputs f ( A ) j and f ( B ) j (since they are independent),
where the S i , S i j , etc., are the widely-known sensitivity indices for first and higher orders. Since this variance decomposition has 2k − 1 terms, it is typical to compute only two sets of k indices: the k so-called ‘first-order’ effects ( S i ) and the k so-called ‘total effects’ S T i , which are defined as:
STi =
(4)
(1)
where
3. Sampling
V i j + · · · + V 12...k
V i = V X i E X ∼i ( Y | X i ) ,
V X i E X ∼i ( Y | X i ) ≈
(i )
the decomposition of the model output variance V (Y ), which holds for independent input variables, is:
V (Y ) =
Let us generate a sample of N points within Ω , the input space of dimension k (we ignore for the moment the issue of the sample design). The coordinates of the N points can be summarised by a matrix A with N rows and k columns. Let us now generate a further set of N points within Ω , independently from the previous set, and summarise the coordinates in a second N × k matrix, B. (i ) We now define the following matrices { A B : i = 1, . . . , k}, where all columns are from A except the ith column (i.e. the ith input variable) which is from B. Saltelli et al. suggest the following estimator for V i [8]:
V X ∼i ( E X i (Y | X ∼i )) V (Y )
.
(3)
Details on total effects can be found in [18]. In the following investigation, we will use the estimation of both the first-order and total sensitivity indices to test the two sampling designs. 2.2. On the implementation of sensitivity estimators First-order and total sensitivity estimates are based on generated samples, subsequent model evaluations, and the use of specific estimators for Eqs. (2) and (3). Here we describe the best available practice to compute S i and S T i .
Fig. 1. The Sobol’ V -function for various coefficient values.
S. Tarantola et al. / Computer Physics Communications 183 (2012) 1061–1072
1063
Fig. 2. Convergence (top) and i-wise errors at N = 8192 (bottom) for A1-1, k = 10, approximate convergence rates in parentheses.
sequences have been chosen. Quasi-random sequences are also known as “low-discrepancy” sequences, having the property of near-uniformity over Ω (strictly speaking a low-discrepancy sequence is defined by a star discrepancy below a certain measure – see [6] for further information). Unlike random numbers, quasirandom points are deterministic sequences that know about the position of previously sampled points and are constructed to avoid the presence of gaps and clusters (discrepancies) as much as possible. They generally outperform crude Monte Carlo sampling in the estimation of multidimensional integrals [23]. The Sobol’ quasi-random sequences satisfy, in particular, the following three main requirements [10]:
• Best uniformity of distribution as N → ∞. • Good distribution for fairly small initial sets. • A very fast computational algorithm. A full description of the Sobol’ sequence will not be presented here since it would require a lengthy discourse. However, the underpinnings of the theory can be found in [9], and an implementation of the algorithm and further practical details can also be found in [24]. For the sampling here, we use a commercial implementation of the Sobol’ sequence, SobolSeq8192 (developed by BRODA Ltd
and available as a C++ code), which can provide Sobol’ sequences up to dimension 8192. Since it is recommended for LSS designs [11], and to allow error estimation, we randomise each Sobol’ sample with a random shift, following [25]. Specifically, we generate a k-length vector, z, such that z → [0, 1]k , and for the jth point from the Sobol’ sequence X˜ j , we generate a shifted X j such that,
X j = X˜ j + z + X˜ j + z where α denotes the greatest integer less than or equal to α . This is equivalent to shifting the whole sample in a random direction, preserving points that move outside the unit hypercube. In this way, the structure of the sample is preserved, but a random element is introduced, which allows for variance estimation. 3.2. The sampling design Once the sampling strategy is decided, the second choice relates to the sampling design, the aim of which is to generate the triplet (i ) of matrices A, B, A B required to feed the model. 3.2.1. The Randomised Sobol’ QMC design In order to generate the matrices A and B, we generate a randomised Sobol’ sample of size N over the unit hypercube of di-
1064
S. Tarantola et al. / Computer Physics Communications 183 (2012) 1061–1072
Fig. 3. Convergence (top) and i-wise errors at N = 8192 (bottom) for A1-1, k = 30, approximate convergence rates in parentheses.
mension 2k. This creates an N × 2k matrix, T , that retains the equidistribution properties of the QMC method. The N × k matrix A is taken as the first k columns of T , with the remaining k columns given as B. The sample matrix X is defined by appending B be(i ) low A and subsequently the set of matrices A B for i = 1, . . . , k, below B. 3.2.2. The LSS design Caflisch, Morokoff and Owen introduced the concept of the “effective dimension” of a function [26], which is, roughly speaking, a measure of either the number of variables (truncation sense), or the orders of interactions (superposition sense) that are causing the majority of output variance. The point is that quite often a small subset of variables is responsible for most of the output variance, therefore the effective dimension of a function is often considerably smaller than k. LSS was introduced by Owen as a means of exploiting this kind of structure in high-dimensional integrals [11], by grouping influential variables into subsets and creating QMC designs in each subset. Owen explains that LSS is most effective when the variables in each subset interact strongly with one another. LSS can be viewed as a combination of Latin hypercube sampling (LHS) and Quasi-Monte Carlo sampling (QMC). In a basic LHS design, the support Ωi of each variable X i is divided into N
j
equally-sized intervals. For the jth interval, a sample X i is randomly drawn across the interval, giving a total of N points, spread somewhat uniformly over Ωi . This process is repeated for each X i , such that an N × k matrix is generated. Finally, the LHS design is created by randomly permuting the values in each column independently. This results in a design with N points that are roughly uniform over each dimension of Ω , and explores the full range of each variable. LSS extends LHS by dividing the k input dimensions into subd sets {s1 , s2 , . . . , sd }, such that r =1 |sr | = k, where |sr | denotes the cardinality of sr . For each subset sr a QMC design is constructed of N points, giving d matrices, each of size N × |sr |. Much like the LHS design, the full LSS design is constructed by randomly permuting the rows (samples) inside each sr , and concatenating them to make the full N × k matrix. In short, while the LHS randomises the order of samples for every dimension, LSS randomises over subsets of dimensions. LHS could therefore be viewed as a special case of LSS with |sr | = 1 ∀r, where sampling within subsets is achieved by random sampling inside equally-spaced bins (rather than QMC). Owen explains that in order to exploit the full potential of LSS, variables that interact strongly should be grouped into subsets together [11]. In practice this can be difficult, since variable interactions are generally unknown until a sensitivity analysis has been performed. To reflect this, in our investigation we examine first the case when
S. Tarantola et al. / Computer Physics Communications 183 (2012) 1061–1072
1065
Fig. 4. Convergence (top) and i-wise errors at N = 8192 (bottom) for A2, k = 10, approximate convergence rates in parentheses.
the ideal division of input variables is known, and second, where the division is unknown (this is explained further in Section 5). We therefore consider a simple case of LSS where there are only two subsets such that s1 + s2 = k. As before, we generate an N × 2k sample matrix using the Sobol’ sequence. This is split into two N × k matrices, A and B, each of which is altered to create the LSS design (randomising rows for each sr ), then the full sample matrix X is obtained in the same way as before. It is worth noting here that one drawback of the LSS design is that it does not allow sequential addition of new points. Extra points cannot be added to an existing sample since it requires a random permutation of coordinates over all sample points. Therefore if the sample size must be increased, a completely new sample must be designed.
to O ( N −1/2 ) (i.e. that of purely random Monte Carlo), only obtaining the higher QMC rate above a certain (high) value of N. This transition value appears to increase with dimensionality, but is also strongly affected by the nature of the integrand (the effective dimensionality) [12,13]. 4. Setting the experiment 4.1. The test function The test function employed in this work has been widely used as benchmark for sensitivity analysis (see [17]) and is defined as:
Y = V ( X 1 , X 2 , . . . , X k, a1 , a2 , . . . , ak ) =
k
vi
i =1
3.2.3. Properties of the sample designs It is worth considering briefly the properties of the designs here. The LSS design has been shown to exhibit “asymptoticallynegligible bias”, so long as randomised-QMC is used [11]. This comes at the expense of a small increase in the variance of the estimate due to the introduced randomisation, but was shown to be preferable. The Sobol’ sequence is a low-discrepancy design that has good uniform properties as N increases (discrepancy follows O ( N −1 logn−1 N ) for N = 2k ) [27]. However, it has been identified that for large dimensionality, convergence of integrals can reduce
with,
vi =
|4 X i − 2| + a i . 1 + ai
The nature of the function V is driven by the dimensionality k as well as by the values of the coefficients ai . These latter determine the relative importance of input variables (i.e., the lower ai the higher the dependence of Y on X i – see Fig. 1) and can therefore be used to generate a wide spectrum of test cases of different levels of complexity.
1066
S. Tarantola et al. / Computer Physics Communications 183 (2012) 1061–1072
Fig. 5. Convergence (top) and i-wise errors at N = 8192 (bottom) for A2, k = 30, approximate convergence rates in parentheses.
4.2. Variations of the test function The test cases used in this work can be grouped according to a classification proposed in [14] which defines three different types of functions (here labelled A, B and C) with different characteristics and different effective dimensionality: Type A includes functions with a few dominant input variables, which could also be involved in interactions (low effective dimension in the superposition and truncation sense). These are the most common functions in practice. Type B includes functions with a homogeneous effect of all input variables, which are involved in some interactions (effective dimension is roughly equal to k in the truncation sense). Type C is like Type B but the input variables are involved in strong interactions (effective dimension is similar to k in both truncation and superposition senses). In this work we present five test cases, classified according to the scheme above, that are representative of a wide range of model behaviours. For each of them we know the analytical values of the sensitivity indices (see Appendix A). We compare those exact values with the estimates obtained through the model evaluations. A1-1/A1-2: The first test case (called A1-1) uses the coefficients a1 = a2 = 0, ai >2 = 6.52; the second test case (A1-2) uses the co-
efficients ak−1 = ak = 0; ai
S i = 0.8268
i =1
from which it follows that 17.32% of the variance in the model output is due to interactions between the input variables. Note that these models have a low effective dimension, which could be exploited by the LSS design. However, there are no strong interactions. A2: The third test case, called A2, uses the coefficients a1 = a2 = · · · = a5 = 0; ai >5 = 6.52, which yield a different decomposition of the sensitivity in terms of first, second and higher-order effects. Here, the first-order effects account for 50.79% of the output variance (k = 10):
S. Tarantola et al. / Computer Physics Communications 183 (2012) 1061–1072
1067
Fig. 6. Convergence rate of Sobol’, LSS and ignorant LSS with model dimensionality for A2 function.
Fig. 7. Comparison of absolute MAES values for A2 function at increasing N and k.
k
S i = 0.5079
i =1
while interactions, at any order, account for the remaining 49.21% of the output variance. Given that the majority of interaction variance is caused by the set {x1 , . . . , x5 }, this should be an ideal case for the LSS design if the subsets are chosen appropriately. B: The fourth test case (called B) uses the coefficients ai = 6.52 for i = 1, . . . , k. The variance is therefore spread evenly over all inputs, mostly in first-order effects. Type B functions are quite rare:
models encountered in practice usually have a few dominant input variables and a larger number of less important ones. C: The final test case C uses the coefficients ai = 0 for i = 1, . . . , k. In contrast with B, the variables have important interactions, such that the sum of first-order effects for k = 10 is k
S i = 0.1989.
i =1
This means that around 80% of the variance in the output is driven by interactions. Again, this is not a common type of function, but is worth investigating for completion.
1068
S. Tarantola et al. / Computer Physics Communications 183 (2012) 1061–1072
Fig. 8. Convergence (top) and i-wise errors at N = 8192 (bottom) for B, k = 10, approximate convergence rates in parentheses.
As will be seen in the following sections, it becomes increasingly difficult to achieve a desired accuracy when estimating sensitivity indices as the sum of first-order interactions becomes smaller (i.e. in case C in particular). 4.3. Sample approach and dimensionality In order to assess the application of these designs to models of various dimensionality, we investigate the cases where k is equal to 10, 20, 30, 50 and 100. For each model, we run the experiments at increasing sample size N to check convergence of the sensitivity estimates to the true values. We set the sample size N at values 16, 32, 64, 128, 256, 512, 1024, 2048, 4096 and 8192. Finally, each experiment at a given N and k is replicated R = 50 times in order to quantify the statistical error of the sensitivity estimates, in terms of deviation between the estimates and the true values. 4.4. Measuring convergence We use four indicators to quantify the statistical error of the sensitivity estimates at given N and k – these are named AES, MAES, AEST and MAEST. Let us denote by S¯ i the analytic value of a first-order index for (r ) the ith input and by S i the corresponding estimated value of the
rth replication. We now define AES as the mean absolute error of each S i :
AESi =
R 1 (r ) S − S¯ i i R r =1
i.e., these help to investigate which model input causes most of the deviation from the analytical values. MAES is the average of the AES values over all inputs to give a single error measure,
1 1 (r ) S − S¯ i . i Rk k
MAES =
R
i =1 r =1
The definitions of AEST and MAEST are exactly equivalent to AES and MAES respectively, but measure errors based on estimates of S T i . For any test function and any k we regress both MAES and MAEST against N by a trend line c ∗ N −a (power regression) in a log–log plot. Values of a ∈ [0, 1] are given in brackets in the legend of the figures and show the rate of convergence of the sensitivity estimates for each implementation design at increasing sample size (larger a implies better convergence properties). For reference, standard Monte Carlo has a = 0.5. The coefficient c is related to the approximation error at small N, therefore small values of c mean that the implementation performs well with a small sample size.
S. Tarantola et al. / Computer Physics Communications 183 (2012) 1061–1072
1069
Fig. 9. Convergence (top) and i-wise errors at N = 8192 (bottom) for B, k = 30, approximate convergence rates in parentheses.
5. Results and analysis The following experiments compare convergence of sensitivity estimates using randomised Sobol’ against LSS designs, where the subsets are chosen with knowledge of the grouping of variables. Since, in practice, the interactions of variables are rarely known a priori, convergence is also calculated for an “ignorant” LSS sample. In this case, we assume no knowledge of the grouping of variables, and choose two subsets such that |s1 | = |s2 | = k/2. The intention is to see how LSS performs with no knowledge of groups, or if the groups are incorrectly specified. For the Type B and Type C functions, there is no “correct” grouping, since all variables are equally important. Therefore we adopt the “ignorant” LSS approach as default. 5.1. Type A1 functions Fig. 2 (top) shows the convergence of estimates with increasing N for both sample designs, in the A1-1 configuration, where k = 10. It is evident that in this case the Sobol’ sampling design is superior to the LSS implementation (with little difference between correctly-specified groups and ignorance): the values of MAES and MAEST are generally smaller and their slope at increasing N is steeper. Sobol’ has the largest values of a, however it should be
noted that at small sample sizes LSS tends to have a lower error. Moreover, AES and AEST show that the main contribution to the statistical error comes from the first two input variables. When k = 30 (Fig. 3, top) MAES values for Sobol’ are again almost always better than those of LSS. Interestingly however, the MAEST values are consistently superior for LSS sampling. This could be due to the deterioration of the Sobol’ design in certain dimensions as k becomes large. In contrast, LSS (if subsets are correctly specified) will not have this problem because each subset uses lower-dimensional Sobol’ sequences. It seems that on average, for this type of function, the Sobol’ design is the better choice. The results of the A1-2 function (where the driving variables are the last two, rather than the first two) were different to those of the A1-1 function. While the MAES curves were similar, the Sobol’ sampling outperformed LSS in estimating S T , which suggests that the ordering of variables is important. 5.2. Type A2 function The Type A2 function is similar to Type A1-1 but with five (instead of only two) important input variables, and a greater proportion of variance due to interactions. In the LSS implementation, these five model input constitute group s1 (this is therefore the
1070
S. Tarantola et al. / Computer Physics Communications 183 (2012) 1061–1072
Fig. 10. Convergence (top) and i-wise errors at N = 8192 (bottom) for C, k = 10, approximate convergence rates in parentheses.
ideal subset for the LSS design). In the case of k = 10, the results (Fig. 4) show that, while there is little difference in the estimation of first-order indices, Sobol’ sampling outperforms LSS in estimation of total effect indices. The results for k = 30 are basically identical (Fig. 5). LSS and Sobol’ designs perform very similarly in terms of MAES, but Sobol’ outperforms LSS at all sample sizes with respect to total effect indices (MAEST). Fig. 6 shows the result of measuring the convergence rate for the A2 function at increasing dimensionality. It is interesting to note here that the rate of convergence is nearly always superior for the Sobol’ design. However, besides the convergence rate, it is also important to consider the absolute error at a given sample value. To further investigate, Fig. 7 shows the difference between values of log(MAES) from the Sobol’ and LSS designs, for increasing values of k and N (essentially the pointwise difference between the convergence curves). The indication here is that Sobol’ sampling outperforms LSS at large sample sizes. For lower sample sizes (below 128) the lack of a clear pattern could be interpreted as noise since the number of sample points is too few to give meaningful estimates of sensitivity indices, particularly for higher-dimensional functions.
5.3. Type B function Examining the Type B function results in Fig. 8 (equal importance in all inputs, and little interaction), it is evident that for k = 10 the Sobol’ sampling design is superior to LSS in terms of MAES and MAEST for larger sample sizes. The convergence rate a for the Sobol’ design is very good (a ≈ 1), whereas the rate for LSS designs is around 0.50, typical of a standard Monte Carlo sampling strategy. This is to be expected – the strength of LSS lies in assembling inputs into interacting subsets, which cannot be achieved here. When k = 30 (see Fig. 9), in general, the convergence rate of the Sobol’ sampling decreases (compared to k = 10). For the estimation of S i , Sobol’ still outperforms LSS at large N, but its efficacy is decreased. However, LSS performs better in estimation of total effect indices. In the case where N = 8192, this is due to a large error in the Sobol’ estimation of S T for variable 20, which takes its MAEST value above that of LSS. Since the Type B function is symmetrical in all dimensions, this must be due to a quirk of the Sobol’ sequence. Upon further investigation, it was found that this “spike” was also present at lower sample sizes (although it appeared on different input variables), which caused the Sobol’ MAEST values to be consistently higher than those of LSS.
S. Tarantola et al. / Computer Physics Communications 183 (2012) 1061–1072
1071
Fig. 11. Convergence (top) and i-wise errors at N = 8192 (bottom) for C, k = 30, approximate convergence rates in parentheses.
5.4. Type C function Fig. 10 shows the performance of the designs on the C-type function, for the case where k = 10. The designs appear similar for MAES, although there is a considerable advantage for the Sobol’ design for MAEST. The results for k = 30 (Fig. 11) are very similar: for MAES the two designs show little difference, whereas for MAEST the Sobol’ design show a fairly clear advantage over LSS. In both designs there are error spikes over certain inputs – these are also found at other sample sizes, though not necessarily on the same variables. 6. Conclusion In this work we have presented five test cases, grouped according to the classification proposed by Kucherenko, that are representative of a wide range of model behaviours. The results here are based on the Sobol’ V-function (indeed, on particular implementations of this function), therefore it is important to stress that the conclusions here may not necessarily be extrapolated out to all functions in each group, though they can hopefully provide an indication. For A-type functions, it was found that the Sobol’ design generally performs better, particularly in the case where there is a larger subset of inputs interacting quite strongly (A2 functions). However,
the LSS design can offer some advantages, although these were observed to be restricted to particular cases such as the A1 function, and only ever with respect to the estimation of S T . A more detailed study of convergence in the A2 function revealed that for all dimensionalities investigated, the Sobol’ design was consistently superior. Furthermore, a comparison of absolute error for given dimensionality and sample size revealed that Sobol sampling almost always performed better. B-type and C-type functions were those that give equal importance to all inputs, with weak and strong interactions respectively. Since LSS is based on identifying a subset of strongly-interacting inputs, it was not expected that LSS would offer any advantage over a Sobol’ design, and indeed this was the case in general. However, one characteristic that was that the Sobol’ design tended to distribute the total estimation error more unevenly across inputs, i.e., low error for most indices, then a relatively high error for one or two inputs. This is thought to be due to a deterioration of the Sobol’ design in high dimensions. The LSS design did not generally suffer from this (despite being based on the Sobol’ design), which can be viewed as an advantage, since it would not be known which of the sensitivity indices were the accurate ones in most practical situations. In general, it is difficult to recommend LSS as a sampling strategy, unless the effective dimension of the function is known and is compatible with the strengths of LSS. It may be that in some cases,
1072
S. Tarantola et al. / Computer Physics Communications 183 (2012) 1061–1072
this information is available from some background knowledge of the function (the physical process behind a model perhaps), or from a preliminary screening test. However in general, little or nothing is known about the response of the model or function. In such cases the Sobol’ sequence provides a good general-purpose sampling approach. Finally, it should be noted that while the results here favour Sobol’ sampling over LSS, they are based on a particular test function which may not necessarily generalise to all model responses in practice. Appendix A. V -function and its analytical values The sensitivity indices for the V-function can be obtained analytically from the following,
V i = V X i E X i (Y | X i ) = VTi = Vi
1 /3 , (1 + ai )2
k (1 + V j ), j =i
V =
k
(1 + V i ) − 1
i =1
where S i = V i / V and S T i = V T i / V . Note that V here denotes variance, not to be confused with the variables of the V -function. The analytical values for the test cases are therefore: k = 30
k = 100
Type A1-1
sum( S i ) S 1 / S Tot 1 S k / S Tot k
k = 10 0.83 0.72 0.54
0.76 0.64 0.48
0.58 0.42 0.32
Type A1-2
sum( S i ) S 1 / S Tot 1 S k / S Tot k
0.83 0.54 0.72
0.76 0.48 0.64
0.58 0.32 0.42
Type A2
sum( S i ) S 1 / S Tot 1 S k / S Tot k
0.51 0.31 0.23
0.47 0.27 0.21
0.35 0.18 0.14
Type B
sum( S i ) S 1 / S Tot 1 S k / S Tot k
0.97 0.95 0.95
0.92 0.84 0.84
0.74 0.56 0.56
Type C
sum( S i ) S 1 / S Tot 1 S k / S Tot k
0.20 0.08 0.08
Functions
∼0 ∼0 ∼0
∼0 ∼0 ∼0
References [1] J.E. Oakley, A. O’Hagan, Probabilistic sensitivity analysis of complex models: a Bayesian approach, Journal of the Royal Statistical Society Series B 66 (3) (2004) 751–769. [2] T. Crestaux, I. Le Ma, O. Tre, J.-M. Martinez, Polynomial chaos expansion for sensitivity analysis, Reliability Engineering & System Safety 94 (7) (2009) 1161–1172.
[3] M. Ratto, A. Pagano, Using recursive algorithms for the efficient identification of smoothing spline ANOVA models, AStA Advances in Statistical Analysis 94 (4) (2010) 367–388. [4] R. Jin, W. Chen, T.W. Simpson, Comparative studies of metamodelling techniques under multiple modelling criteria, Structural and Multidisciplinary Optimization 23 (1) (2001) 1–13. [5] C.B. Storlie, L.P. Swiler, J.C. Helton, C.J. Sallaberry, Implementation and evaluation of nonparametric regression procedures for sensitivity analysis of computationally demanding models, Reliability Engineering & System Safety 94 (11) (2009) 1735–1763. [6] R.E. Caflisch, Monte Carlo and quasi-Monte Carlo methods, Acta Numerica 7 (1998) 1–49. [7] H. Niederreiter, Random Number Generation and Quasi-Monte Carlo Methods, Society for Industrial Mathematics, 1992. [8] 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, Computer Physics Communications 181 (2) (2010) 259–270. [9] I.M. Sobol’, On the distribution of points in a cube and the approximate evaluation of integrals, USSR Computational Mathematics and Mathematical Physics 7 (4) (1967) 86–112. [10] I.M. Sobol’, Uniformly distributed sequences with an additional uniform property, USSR Computational Mathematics and Mathematical Physics 16 (5) (1976) 236–242. [11] A.B. Owen, Latin supercube sampling for very high-dimensional simulations, ACM Transactions on Modeling and Computer Simulation 8 (1) (1998) 71–102. [12] W.J. Morokoff, R.E. Caflisch, Quasi-random sequences and their discrepancies, SIAM Journal on Scientific Computing 15 (6) (1994) 1251–1279. [13] W.J. Morokoff, R.E. Caflisch, Quasi-Monte Carlo integration, Journal of Computational Physics 122 (2) (1995) 218–230. [14] S. Kucherenko, B. Feil, N. Shah, W. Mauntz, The identification of model effective dimensions using global sensitivity analysis, Reliability Engineering & System Safety 96 (4) (2011) 440–449. [15] A. Saltelli, et al., Global Sensitivity Analysis: The Primer, John Wiley & Sons, Ltd., Chichester, 2008, x, 292 p. [16] R.I. Cukier, C.M. Fortuin, K.E. Shuler, A.G. Petschek, J.H. Schaibly, Study of the sensitivity of coupled reaction systems to uncertainties in rate coefficients. I. Theory, Journal of Chemical Physics 59 (8) (1973) 3873–3879. [17] I.M. Sobol’, Sensitivity estimates for nonlinear mathematical models, Mathematical Modeling and Computational Experiment 1 (4) (1993) 407–414. [18] T. Homma, A. Saltelli, Importance measures in global sensitivity analysis of nonlinear models, Reliability Engineering & System Safety 52 (1) (1996) 1–17. [19] A. Saltelli, S. Tarantola, F. Campolongo, M. Ratto, Sensitivity Analysis in Practice: A Guide to Assessing Scientific Models, John Wiley & Sons Inc., 2004. [20] I.M. Sobol’, S. Tarantola, D. Gatelli, S.S. Kucherenko, W. Mauntz, Estimating the approximation error when fixing unessential factors in global sensitivity analysis, Reliability Engineering & System Safety 92 (7) (2007) 957–960. [21] M.J.W. Jansen, Analysis of variance designs for model output, Computer Physics Communications 117 (1–2) (1999) 35–43. [22] M.J.W. Jansen, W.A.H. Rossing, R.A. Daamen, Monte Carlo estimation of uncertainty contributions from several independent multivariate sources, in: G. van Straaten, J. Grasman (Eds.), Predictability and Nonlinear Modelling in Natural Sciences and Economics, Kluwer Academic Publishing, Dordrecht, 1994, pp. 335–343. [23] I.M. Sobol’, S.S. Kucherenko, On global sensitivity analysis of quasi-Monte Carlo algorithms, Monte Carlo Methods and Applications 11 (1) (2005) 83–92. [24] P. Bratley, B.L. Fox, ALGORITHM 659: implementing Sobol’s quasirandom sequence generator, ACM Transactions on Mathematical Software 14 (1) (1988) 88–100. [25] R. Cranley, T.N.L. Patterson, Randomization of number theoretic methods for multiple integration, SIAM Journal on Numerical Analysis 13 (6) (1976) 904– 914. [26] R.E. Caflisch, W. Morokoff, A. Owen, Valuation of mortgage backed securities using Brownian bridges to reduce effective dimension, Journal of Computational Finance 1 (1) (1997) 27–46. [27] I.M. Sobol, On the systematic search in a hypercube, SIAM Journal on Numerical Analysis 16 (5) (1979) 790–793.