Statistics and Computing (2000) 10, 187–196

Winding Stairs: A sampling tool to compute sensitivity indices KAREN CHAN, ANDREA SALTELLI and STEFANO TARANTOLA Institute for Systems, Informatics and Safety, European Commission Joint Research Centre, Italy [email protected]

Accepted October 1999

Sensitivity analysis aims to ascertain how each model input factor influences the variation in the model output. In performing global sensitivity analysis, we often encounter the problem of selecting the required number of runs in order to estimate the first order and/or the total indices accurately at a reasonable computational cost. The Winding Stairs sampling scheme (Jansen M.J.W., Rossing W.A.H., and Daamen R.A. 1994. In: Gasman J. and van Straten G. (Eds.), Predictability and Nonlinear Modelling in Natural Sciences and Economics. pp. 334–343.) is designed to provide an economic way to compute these indices. The main advantage of it is the multiple use of model evaluations, hence reducing the total number of model evaluations by more than half. The scheme is used in three simulation studies to compare its performance with the classic Sobol’ LPτ . Results suggest that the Jansen Winding Stairs method provides better estimates of the Total Sensitivity Indices at small sample sizes. Keywords: first order indices, Sensitivity Analysis, Sobol’ LPτ sequences, Total Sensitivity Indices, Winding Stairs sampling scheme

1. Introduction This paper describes a sampling scheme used to perform a Monte-Carlo based sensitivity analysis (SA) of model output. In general, SA is conducted by the following steps: (i) defining a model which approximates economic, engineering, environmental, physical or social phenomenon of various levels of complexity, and its input factors and output variable(s); (ii) assigning probability distributions to the input factors; (iii) generating input values via an appropriate random sampling method and evaluating the output; and (iv) assessing the influence or relative importance of each input factor on the output variable. The focus of this article is on steps (iii) and (iv). First we will introduce some standard notation. We assume a model output Y is dependent on k input factors, say x1 , x2 , . . . , xk , namely y = f (x1 , x2 , . . . , xk ). Difficulty can arise from within a modelling process when the model factors vary about nominal values in some manner which is unknown to us. In this situation we model uncertainties about C 2000 Kluwer Academic Publishers 0960-3174 °

the input factors by treating them as random variables. Hence, we often face the problem of what values to use for the input factors when we investigate real world phenomena with numerical experiments using a mathematical model. In SA, several values of x = (x1 , x2 , . . . , xk ), say x1 , x2 , . . . , xn , are generated as successive sets of inputs in order to obtain the desired and accurate information concerning Y . For some complex models, the computational cost can be expensive in that a single set of input values may require several minutes in order to yield one output. In practice, the sample size, n, should be large enough to provide accurate information about the relationships between the input factors and the model output, while it should be small enough to minimise the computational cost. Several techniques to generate sample input points, such as Crude Monte Carlo, Latin hypercube, Sobol’ LPτ , have been proposed and compared (Homma and Saltelli 1995, McKay, Beckman and Conover 1979). The next section describes two sampling methods for selecting (sampling) input values of the input factors. In Section 3, methods for assessing the influence or relative importance of each input factor on the output variable are given. Three simulation studies and their results are described and presented in Section 4. Finally, we summarise our conclusions and discussion in Section 5.

188

Chan, Saltelli and Tarantola

2. Sampling schemes 2.1. The Winding Stairs sampling method In ordinary Monte Carlo (MC) sampling a new realisation of the model output y is obtained by drawing values for the inputs according to their joint probability distribution, and calculating y for each multivariate realisation x. Such a sample contains no information about the role of the individual input factors. The winding stairs method (Jansen, Rossing and Daamen 1994) consists of calculating y after each drawing of a new value of an individual factor, xi for i = 1, 2, . . . , k. Each new value is drawn at random from the marginal distribution of X i . A sequence of sample points is generated as follows: we generate a value from each input factor space to obtain a sample point in the k dimensional space, e.g. {x11 , x21 , x31 , . . . , xk1 }. (Note that the first suffix denotes parameter number and the second suffix denotes sample number.); we then sample a value from the sample space of the second parameter to obtain a second sample point, namely {x11 , x22 , x31 , . . . , xk1 }. Similarly, a third sample point is obtained by sampling a value from the sample space of the third parameter giving {x11 , x22 , x32 , . . . , xk1 }, and so on. Hence, the kth sample point is {x11 , x22 , x32 , . . . , xk2 } and the next sampling point is obtained by changing the value of the first parameter giving, {x12 , x22 , x32 , . . . , xk2 }. We repeat these steps until we obtain the desired number of sample points. This means that new input values are sampled in a fixed cyclic order. The output is evaluated after each sample input point is generated, yielding a sequence of output yl , for l = 1, 2, . . . , N , where N is the total number of model evaluations. We arrange the sequence of N outputs into k columns and r rows, where r (= N /k) is the number of turns to repeat the cyclic order. For example, the following output matrix shows a set of outputs y1 , y2 , . . . , y15 for k = 3 and r = 4,   y1 y2 y3    y4 y5 y6     y7 y8 y9       y10 y11 y12  y13 y14 y15   f (x11 , x21 , x31 ) f (x11 , x22 , x31 ) f (x11 , x22 , x32 )    f (x12 , x22 , x32 ) f (x12 , x23 , x32 ) f (x12 , x23 , x33 )     =  f (x13 , x23 , x33 ) f (x13 , x24 , x33 ) f (x13 , x24 , x34 )  .    f (x14 , x24 , x34 ) f (x14 , x25 , x34 ) f (x14 , x25 , x35 )  f (x15 , x25 , x35 ) f (x15 , x26 , x35 ) f (x15 , x26 , x36 ) The entries within each column are independent of each other (e.g. y1 = f (x11 , x21 , x31 ) and y4 = f (x12 , x22 , x32 )), since the values of each input factor are all different. However, the consecutive points within each row are not independent in the

Fig. 1. A plot of 63 Winding Stairs sample points and the trajectory of first 26 points for three factors

sense that the points differ in a single input factor value (e.g. the output y3 and y4 are computed from f (x11 , x22 , x32 ) and f (x12 , x22 , x32 ), respectively). In total, we generate k × (r + 1) input points and the sequence of these input points forms a trajectory in the input factor space. Figure 1 shows an example of a trajectory of the first 26 sample points generated by the WS sampling scheme in a three factor case; the total number of points plotted is 63 (i.e. r = 20). Note that a pseudo-random number generator is used to obtain each input factor value. 2.2. Sobol’ LP τ sampling scheme Sequences of LPτ vectors represent a strategy to produce sample points uniformly distributed in a unit cube, and are essentially quasi-random sequences which are defined as sequences of points that have no intrinsic random properties. Generally, sequences of quasi-random vectors, Q 1 , Q 2 , . . . , Q n , should fulfill the following requirements (Sobol’ 1994): (i) The uniformity of the distribution is optimal when the length of the sequence tends to infinity. (ii) The uniformity of vectors Q 1 , Q 2 , . . . , Q n should be observed for fairly small n. (iii) The algorithm used for the computation of the vectors should be simple. Quasi-random sequences are used in place of random points to guarantee convergence of estimation in the classical sense. The use of points of LPτ sequences usually results in better convergence when employed instead of random points in a Monte Carlo algorithm with finite constructive dimension,1 n. The LPτ sequences were introduced by Sobol’ in 1966 and obey the requirements listed above. The theory underlying the Sobol’s LPτ sequences is given in Sobol’ (1967, 1976) and the algorithm to generate the sequences has been coded in FORTRAN77 and C (Bratley and Fox 1988, Sobol’ et al. 1992). These programs are widely available hence we will not describe the algorithm here. For readers who wish to know more about the theoretical development of LPτ sequences,

Winding Stairs: A sampling tool

189

Fig. 2. Filling-in points generated by LP τ sequences with (a) 512 points, and (c) 1024 points, and (d) 4096 points and (b) 512 points generated by pseudo-random algorithm

a good summary description can be found in Bratley and Fox (1988). With as little as 512 points, the uniformly-regular filling-in feature of the LPτ sequence can clearly be seen (Fig. 2(a)), compared with the purely random placement of pseudo-random points (Fig. 2(b)). The structure of the image of quasi-random points changes as the number of points increases but the property of regularity remains (Fig. 2(c) and (d)). However, sequences that are intended for numerical computations must satisfy the following additional property Z 1 Z 1 n 1X ··· g(x1 , x2 , . . . , xk ) d x1 · · · d xk = lim g(Q j ), n→∞ n 0 0 j=1 where g(x1 , x2 , . . . , xk ) is an arbitrary integrable function.

The main reason that investigators favour the use of quasirandom points is the fast rate of convergence; for large n, the approximate error can be of the order 1/n compared with or√ der 1/ n for standard Monte Carlo methods. Thus, without changing the computation algorithm, but merely by replacing the random numbers with coordinates of quasi-random points, we can improve our results considerably.

3. Sensitivity Analysis The ultimate goal in performing SA is to investigate the relative importance of each input factor with respect to the model output variation. This could be done by measuring the main effect and/or total effect of each individual input factor on the model

190

Chan, Saltelli and Tarantola

output. There is a large literature on these topics, the most recent are Sobol’ (1993), Jansen, Rossing and Daamen (1994), Homma and Saltelli (1996), McKay (1997), Saltelli and Bolado (1998), Saltelli, Tarantola and Chan (1999), and Jansen (1999). In the next two sections, we will discuss how to measure the main effects and the total effects of the input factors.

integrals, namely Z

1

f i (xi ) = − f 0 +

Z

1

···

0

f (x) dx∼i 0

Z

1

f i j (xi , x j ) = − f 0 − f i (xi ) − f j (x j ) + 0

3.1. First order Sensitivity Indices, SI The main effects can be measured by the so-called partial variance (Sobol’ 1993) or correlation ratios (McKay 1997), or Top Marginal Variance (TMV) (Jansen Rossing and Daamen 1994), which is defined to be the expected variance reduction due to fixing factor X i while the remaining X∼i vary. Here, x∼i denotes a vector of input values x excluding the input value for factor X i . Here we will describe briefly the Sobol’ and Jansen’ methods.

Äk = {x | 0 ≤ xi ≤ 1; i = 1, . . . , k}. The function f (x) is decomposed into summands of increasing dimensionality, namely f (x1 , . . . , xk ) = f 0 +

k X i=1

f i (xi ) +

X

f i j (xi , x j ) (1)

A more general representation of this decomposition using multiple integrals can be summarized as follows. For (1) to hold f 0 must be a constant, and the integrals of every summand over any of its own variables must be zero, i.e. Z

1

¡

¢

f i1 ,...,is xi1 , . . . , xis d xik = 0, if 1 ≤ k ≤ s.

f (x) dx∼(i j) 0

Äk

while the second order moments of the terms in (1) are given by Z 1 Z 1 ··· f i21 ,...,is (x1 , . . . , xs ) d xi1 , . . . , d xis Di1 ,...,is = 0

where 1 ≤ i 1 < · · · < i s ≤ k and s = 1, . . . , k. Since we treat the point x as random and uniformly distributed in Äk , then f (x) and f i1 ,...,is (x1 , . . . , xs ) are also random, hence D and Di1 ,...,is are the measures of their variances. The partial variance for factor X i is given by ¸2 Z ·Z Z ··· f (x) d x∼i d xi − f 02 . (5) Di = The sensitivity indices develop very naturally, that is by squaring and integrating (1) over Äk and using (3) we obtain

1≤i< j≤k

+ · · · + f 1,2,...,k (x1 , . . . , xk ).

1

with the convention that dx∼i , dx∼(i j) denote integration over all variables except xi , and xi and x j , respectively. Analogous formulae can be obtained for the higher order terms. The total variance D of f (x) is defined to be Z D= f 2 (x) dx − f 02 (4)

0

3.1.1. The Sobol’ method In order to describe the Sobol’ method (Sobol’ (1990 in Russian and 1993 in English)), let us define the input factor space Äk as a k-dimensional unit cube, i.e. the region

Z ···

(2)

0

A consequence of (1) and (2) is that all the summands in (1) are orthogonal, i.e. if (i 1 , . . . , i s ) 6= ( j1 , . . . , jl ), then Z f i1 ,...,is f j1 ,..., jl dx = 0. (3) Äk

D=

Äk

Sobol’ (1993) showed that decomposition (1) is unique and that all the terms in (1) can be evaluated via multidimensional

X

Di +

i=1

Di j + · · · + D1,2,...,k

1≤i< j≤k

Hence, the sensitivity measure Si1 ,...,is are given by Si1 ,...,is =

Di1 ,...,is for 1 ≤ i 1 < · · · < i s ≤ k D

(6)

where Si is called the first order sensitivity index for factor xi , which measures the main effect of xi on the output (the fractional contribution of xi to the variance of f (x)), Si j , for i 6= j, is called the second order sensitivity index which measures the interaction effect (the part of the variation in f (x) due to xi and x j which cannot be explained by the sum of the individual effects of xi and x j ), and so on. It can easily be seen that all the terms in (6), sum to 1, namely k X

Since at least one of the indices will not be repeated, the corresponding integral will vanish due to (2). Another consequence is that Z f (x) dx. f0 =

k X

i=1

Si +

X

Si j + · · · + S1,2,...,k = 1.

1≤i< j≤k

3.1.2. The Jansen method Consider the analysis of variance decomposition of f with two groups of input factors, X = {U, V } say. Then, by equation (1), we have f (x) = f 0 + f 1 (u) + f 2 (v) + f 12 (u, v)

Winding Stairs: A sampling tool with

Z

Z f 1 du =

Let

Z f 2 dv =

Z

Du =

191

f 12 du =

f 12 dv = 0.

Z f 12 du,

Dv =

where S∼i is the sum of all the indices which do not include factor X i , and is computed by

Z

Dˆ ∼i . Dˆ

Z f 22 dv,

Duv =

2 f 12 du dv.

The Jansen’ method can also be applied to compute the TSI. It can be shown (see Appendix) that 1 E[ f (X i , X∼i ) − f (X i0 , X∼i )]2 2 = D − Cov [ f (X i , X∼i ), f (X i0 , X∼i )].

Then, the total variance of the output is given by D = Du + Dv + Duv .

(7)

Jansen, Rossing and Daamen (1994) proposed to measure the main effect of a group of factors u, by Du = D −

1 E[ f (U, V) − f (U, V0 )]2 , 2

The right hand side of equation (11) can be written in terms of an MC integral, namely Cov [ f (X i , X∼i ), f (X i0 , X∼i )] Z = f (xi , x∼i ) f (xi0 , x∼i ) d xi d xi0 dx∼i − f 02 .

(8)

where D is the output variance. This squared difference mea˘ sure has been described in Saltenis and Dzemyda (1982). Jansen (1996) showed that D−

1 E[ f (U, V) − f (U, V0 )]2 = Cov [ f (U, V), f (U, V0 )] 2 (9)

We shall see later that, from the computational point of view, the right hand of equation (9) is equivalent to that of the Sobol’ first order indices. By partitioning X into X i and X ∼i , we can measure the main effect of the individual factor, X i . For example, the first order index of factor X i , denoted by DiJ here, is given by DiJ = D −

1 E[ f (X i , X∼i ) − f (X i , X0∼i )]2 . 2

= D∼i .

(12)

Hence, the total effect of factor X i can be measured by DiJtot =

1 E[ f (X i , X∼i ) − f (X i0 , X∼i )]2 . 2

(13)

Again from the computational point of view, Jansen’s method to compute TSI is equivalent to that of Sobol’. In summary, for both indices Jansen’s method uses the squared differences of two sets of model outputs; while the Sobol’ method uses the product. In the next section, we will discuss how these indices are computed.

(10)

In the other words, the expected squared difference on the right hand side of equation (10) is computed by fixing factor X i and varying the remaining factors X∼i . 3.2. Total Sensitivity Indices, TSI The Total Sensitivity Index (TSI) is a measure of the total influence of an input factor on the model output variation, and is defined as the sum of all the sensitivity indices (including all the interaction effects) involving the factor of interest (Homma and Saltelli 1996). From equation (7), the total effect of factor X i is given by Ditot = Di + Di(∼i) = D − D∼i , where D∼i is the total fractional variance complementary to factor X i , and is given by Z f (xi , x∼i ) f (xi0 , x∼i ) d xi d xi0 dx∼i − f 02 . D∼i = Hence, from equation (6), the TSI of factor X i , denoted by STi , is given by TS i = 1 − S∼i ,

(11)

3.3. Computational issues The computation of the Sobol’ sensitivity indices is usually associated with the LPτ quasi-random numbers, described in Section 2.2. The partial variance for factor X i given in equation (5) can be computed by performing Monte Carlo integrals. Two sampling data matrices M1 and M2 both of dimension n × k for the input factors X are generated. The MC estimates of quantity Di and Ditot are given by n ¡ 1X (M1) ¢ ¡ (M2) ¢ f xim , x∼im f xim , x∼im − fˆ02 , Dˆ i = n m=1

(14)

and

" # n ¡ (M1) ¢ ¡ (M2) ¢ 1X tot 2 ˆ ˆ f xim , x∼im f xim , x∼im − f 0 , Di = D − n m=1 (15)

where n is the sample size. For example, Dˆ i is computed by multiplying the output corresponding to x from matrix M1 with the output computed from a different matrix M2 with the ith column taken from matrix M1. The matrix M1 is usually called the data base matrix and M2 the resampling matrix. The sample variance of the output, D given in equation (4), is estimated

192

Chan, Saltelli and Tarantola

using the base matrix, and is given by 1 Dˆ = n

n X

f 2 (xm ) − fˆ02 ,

(16)

m=1

P where fˆ0 = n1 nm=1 f (xm ) is an estimate of the mean of the model output. So far, Jansen’s method has been associated with the Winding Stairs sampling scheme. Recall that the outputs within each column of the WS output matrix are independent of each other, so D can be estimated by pooling all the k sample variances of the r independent outputs. Hence, the WS sample estimate of D is given by Dˆ WS

k X 1 = k(r − 1) i=1

(

r X

"

#2 )

r 1X y 2 (m, i) − y(m, i) r m=1 m=1

, (17)

where {y(m, i)} denotes the (m, i)th element of the output matrix Y. The expected squared difference on the right hand side of equation (10) is computed by fixing factor X i and varying the remaining factors X∼i . Hence, it can be estimated by using the (r + 1) × k WS output matrix as described in Section 2.1, namely r £ ¤2 1 X Dˆ iWS = Dˆ WS − yk(l−1)+i − ykl+i−1 , 2r l=i

(18)

where Dˆ iWS denotes the estimate of Di using the Jansen’s method with the Winding Stairs sampling scheme. For example, with k = 3 and r = 4, to compute the SI of factor X 1 , we calculate the sum of the squared differences between the entries of the first and the third column of the first four rows of the output matrix; however, for factor X 2 , the squared differences are computed using the values in column 1 and column 2 with a row shifted. The corresponding WS sample estimate of the SI is given by Dˆ WS SˆiWS = 1 − i for i = 1, 2, . . . , k. Dˆ WS Similarly, the sum of all the partial variances which are complementary to factor X i , i.e. D∼i is computed by varying factor X i and fixing the remaining factors X ∼i . By using the WS output matrix and equation (13), the WS estimates of Ditot for factor X i is given by  r 1 X   [ylk − ylk+1 ]2 if i = 1     2r l=1 (19) DiWStot =  r  X £ ¤ 1  2   if i 6= 1. yk(l−1)+i−1 − yk(l−1)+i  2r l=1 For example, with k = 3 and r = 4, the TSI of Factor X 1 is estimated by computing the sum of the squared differences between the entries of the first and the last column of the output matrix with a row shifted; however, the TSI of factor X 2 is estimated

by taking the squared differences between the entries of the first and the second column of the output matrix. It is easy to see that both the LPτ and the Winding Stairs sampling schemes can be used to compute the Sobol’ and Jansen’s sensitivity indices. 3.4. Sobol’-LP τ vs Jansen-WS A drawback of the Sobol’ method is that a separate computation is needed for each factor. To compute both sets of first-order and total indices, the Sobol’ method requires N = n(2k + 1) model evaluations, while the number of model evaluations required for the WS method is N = r k, where r denotes the ‘sample size’ for the WS method. The efficiency of the winding stairs method lies mainly in the multiple use of model evaluations in that the first order and total sensitivity indices of each factor can be computed using a single set of model evaluations, each model evaluation is used approximately four times. With the same sample sizes (i.e. r = n) and computing all the first order and total sensitivity indices, the WS sampling method reduces the total number of model evaluations by more than a half. In the next section we will present results of three simulation studies to compare the Sobol’ indices with LPτ sampling scheme and Jansen’s sensitivity measures with Winding Stairs.

4. Simulation studies In this section we describe three simulation studies using two analytical test functions to illustrate Jansen’s method with the WS sampling scheme (referred to hereafter as the Jansen-WS method or simply the WS method) to estimate the first order and total sensitivity indices, and to compare the method with the Sobol’-LPτ . 4.1. Example 1 This first example uses Legendre polynomials which have been used to compare different SA techniques (McKay 1997, Saltelli, Tarantola and Chan 1999). The Legendre polynomial of order d has two input factors, namely, x (a uniformly distributed random variable taking values over the range [−1, +1]) and d (a discrete uniformly distributed random variable, assumed to take values from 1 to 5). The analytical values of the partial variances due to the factors d and x (Dd and Dx ) are given in Table 3 of McKay (1997). These analytical values are given as conditional expectations and variances, e.g. Var [E(Y | x)] = Dx , which is the first-order partial variance of Y due to factor x and since Var [Y ] = Var [E(Y | x)] + E(Var [Y | x]), we have that E(Var [Y | x]) is the partial variance complementary to x, namely, D∼x + Dx,∼x .

Winding Stairs: A sampling tool

193

Fig. 3. Plot of Jansen-WS and Sobol’-LP τ estimates of (a) first order and (b) total sensitivity indices against number of model evaluations. Analytical values of the respective indices for factors x and d are Sx = 0.2 & Sd = 0.0, and TS x = 1.0 & TS d = 0.8.

The Jansen-WS and Sobol’-LPτ estimates of the first order and total sensitivity indices are plotted against the number of model evaluations (correspond to sample sizes of 32, 64, 128, 256, 512, 1024 and 2048) in Fig. 3(a) and (b), respectively. With the same number of model evaluations at small sample sizes, the Jansen-WS provides better estimates for both first and total indices than the Sobol’ method. To examine robustness of the Jansen-WS method, the experiment is repeated 100 times, the average, minimum and maximum of the 100WS estimates of the first order and total indices are plotted against the number of model evaluations. For both indices, as the number of model evaluations increases the ranges decrease and the averages quickly converge to the analytical values.

reduces to a function used by Davis and Rabinowitz (1984) as a test case for performing multidimensional integration. Note R1 that, for each gi (xi ), 0 gi (xi ) d xi = 1 for i = 1, 2, . . . , k, and R1 R1 hence 0 · · · 0 f d x1 · · · d xk = 1. Furthermore, it can be shown that the range of the function gi (xi ) is given by 1−

1 1 ≤ gi (xi ) ≤ 1 + . 1 + ai 1 + ai

Hence, the values of the ai ’s determine the relative importance of the xi ’s since the range of uncertainty of gi (xi ) depends exclusively on the value of ai . The g-functions for four ai values are plotted in Fig. 4. The corresponding relative important of the

4.2. Examples 2 and 3 These two examples are illustrated by the Sobol’ g-function, which has been used by several investigators (Saltelli and Sobol’ 1995, Archer, Saltelli and Sobol’ 1997, and Saltelli and Bolado 1998) to compare different techniques for computing sensitivity indices. The Sobol’ g-function is defined in the k-dimensional unit cube, as follows: f =

k Y

gi (xi )

i=1

where gi (xi ) is given by gi (xi ) =

|4xi − 2| + ai , for 0 ≤ xi ≤ 1 and ai ≥ 0. 1 + ai

The factor ai is chosen to specify the role of the corresponding input factor xi . For ai = 0, the function, displayed in Fig. 4,

Fig. 4. Plot of g-function vs the input factor x for ai = 0 (solid), ai = 1 (dotted), ai = 9 (dashed) and ai = 99 (long dashed). The relative importance of the corresponding input factors are “Very important”, “Important”, “Non-important”, and “Non-significant”, respectively

194

Chan, Saltelli and Tarantola

Table 1. Analytical values of the first order and total indices for the g-function with six and eight factors Sensitivity Index Parameter (ai ) First order Total 0.2 0.8 0.0 0.1 1.6 0.4

0.1261 0.0560 0.1816 0.1550 0.0269 0.0926

0.2904 0.1441 0.3862 0.3336 0.0726 0.2245

Sensitivity Index Parameter (ai )

First order Total

0.0 1.0 4.5 9.9 99.0 99.0 99.0 99.0

0.7162 0.1790 0.0237 0.0072 0.0001 0.0001 0.0001 0.0001

0.7871 0.0242 0.0034 0.0105 0.0001 0.0001 0.0001 0.0001

input factor can clearly be seen. Note that gi is non-monotonic in all its input variables and its derivative changes sign within the interval of variation, hence it is undefined at its midpoint. The exact (analytical) value of the sensitivity index for each input factor can be computed (Saltelli and Sobol’ 1995). In this example, the numbers of input factors chosen are eight and six. The respective sets of ai values are chosen to be {0, 1, 4.5, 9, 99, 99, 99, 99} and {0.2, 0.8, 0.0, 0.1, 1.6, 0.4}. The first set represents the input factors in order of descending importance, while the second set represents the factors as more or less equally important. Saltelli and Sobol’ (1995) showed that the case where all factors are equally important or unimportant is the most difficult case to deal with. The analytical values of the first order and of the total sensitivity indices for these two examples are given in Table 1. The experiment is performed repeatedly so that the accuracy and precision of the two methods can be assessed and compared. The number of replicates is 100. For the Sobol’ LPτ method, the

replicates are obtained by permuting the columns of the data matrix. Comparisons are based on the Total Absolute Error criterion (which we abbreviate here as TAE), namely TAE =

k X

|Il − Iˆl |

l=1

where Il and Iˆl are the analytical value of either the first-order or the total sensitivity index and its estimate for the factor l, respectively. This absolute difference measure is particularly useful if one wants to compare methods resulting from a wide range of models, e.g., linear, non-linear or additive, etc. In order to make a fair comparison between the two methods, both sets of sensitivity indices are computed with the same number of model evaluations. The average and the standard deviation (expressed as average ± standard deviation) of the TAEs of the first order and the total sensitivity indices are plotted against the number of model evaluations for the eight factor case in Fig. 5(a) and (b), respectively. As expected, for both methods and both indices, the averages and the ranges of TAE’s decrease as the number of model evaluations increases, indicating that the estimates converge to the analytical values and their precision increases, as the sample size increases. The Sobol’-LPτ method gives a better estimate of SI than the Jansen-WS method. However, for the TSI, the Sobol’-LPτ method performs worse than the Jansen-WS method at small sample size, and seems to perform marginally better than the Jansen-WS method when the sample size is large, but the differences are relatively small (see Fig 5(b)). Results of the six factor case are plotted in Fig. 6(a) for the first order indices and in Fig. 6(b) for the total indices. As in the eight factor example, for both methods and both set of indices, the TAE’s decrease as the number of model evaluations increases

Fig. 5. Comparisons between the Jansen-WS and Sobol’-LP τ methods using TAE criterion at different numbers of model evaluations: (a) First order indices and (b) Total indices; g-function with eight factors

Winding Stairs: A sampling tool

195

Fig. 6. Comparisons between the Jansen-WS and Sobol’-LP τ methods using TAE criterion at different numbers of model evaluations: (a) First order indices and (b) Total indices; g-function with six factors

and the Sobol’-LPτ method gives a better estimate of SI than the Jansen-WS method. However, for TSI the Jansen-WS method performs better than the Sobol’-LPτ method at all sample sizes, and is considerably better at much lower sample size (n = 32 and 64). Recall that this example represents a difficult case as all the input factors are more or less equally important. The Jansen-WS method outperforms the Sobol’-LPτ method in computing the Total sensitivity indices at very small sample sizes.

5. Conclusion and discussion In this paper, we have described the Winding Stairs sampling scheme and Jansen’s sensitivity measures. Jansen’s methods can be shown to be closely related to those of Sobol’. From the computational viewpoint these methods are equivalent. However, it can be shown (see Chan et al. (2000) for example) that ¤ £ ¤ £ (20) Var DiJtot ≤ Var Ditot , while

£ ¤ Var DiJ ≥ Var [Di ].

It should be noted that the winding Stairs sampling scheme is similar to that of the One-At-Time Morris method (Morris 1991). The Morris method consists of generating r independent trajectories, each of (k + 1) successive points. The successive points are generated by changing a value of one factor at a time with equal steps so that the resulting points are equally spaced. Hence, the Morris method involves dividing the input factors space into a k-dimensional p-level equally spaced grid. The sample points are used to compute two sensitivity measures for each factor: the main effect and a measure which is the sum of all the curvatures and interactions involving that factor. Like the WS scheme, the Morris method makes multiple use of the model evaluations, reducing the number of model evaluations by half.

Appendix Define E X [ f (X)] = f 0 , Var X [ f (X)] = E[ f (X) − f 0 ]2 . Then 1 E[ f (X i , X∼i ) − f (X i0 , X∼i )]2 2 = Cov [ f (X i , X∼i ), f (X i0 , X∼i )]

D− (21)

In the other words, Jansen’s method provides better estimates for measuring the total effects while the Sobol’ method gives better estimates for the main effects. Our numerical examples have confirmed this. The advantage of the WS sampling scheme is the multiple use of model evaluations, giving a saving of more than half and hence, for large problems, reducing the computational burden. However, equations (20) and (21) might also suggest the superiority of the Jansen-WS method could be due to the method of computing the sensitivity measures and not to the WS sampling scheme. What remains, however, is the investigation of the Jansen’s method with the LPτ sequences and the Sobol’ method with the WS sampling scheme.

Proof: 1 E[ f (X i , X∼i ) − f (X i0 , X∼i )]2 2 1 = E[ f (X i , X∼i ) − f 0 − f (X i0 , X∼i ) + f 0 ]2 2 1 = E[( f (X i , X∼i ) − f 0 )2 − 2( f (X i , X∼i ) − f 0 ) 2 × ( f (X i0 , X∼i ) − f 0 ) + ( f (X i0 , X∼i ) − f 0 )2 ] 1 = [D − 2E[( f (X i , X∼i ) − f 0 )( f (X i0 , X∼i ) − f 0 ) + D] 2

196

Chan, Saltelli and Tarantola (since E[ f (X i , X∼i ) − f 0 ]2 = D) = D − Cov [ f (X i , X∼i ), f (X i0 , X∼i )]

¤

Acknowledgment The authors wish to thank Prof. Sobol’, two referees and the associate editor for their helpful comments and suggestions. Also, the authors would like to acknowledge the financial support provided by Eurostat through the SUPCOM Lot 14 (Contract No. 13592).

Note 1. The constructive dimension of a modelling algorithm is the maximum number of points used to compute one realisation of the random variable 3, for instance (Sobol’ 1994).

References Abramowitz M. and Stegun I.A. 1970. Handbook of Mathematical Functions. New York, Dover Publications. Archer G., Saltelli A., and Sobol’ I.M. 1997. Sensitivity measures, ANOVA like techniques and the use of bootstrap. Journal of Statistical Computation and Simulation 58: 99–120. Bratley P. and Fox B.L. 1988. Algorithm 659 implementing Sobol’ quasirandom sequence generator. ACM Transactions on Mathematical Sofware 14: 88–100. Chan K., Saltelli A., Tarantola S., and Sobol’ I.M. 2000. Variance based methods. In: Saltelli A., Chan K., and Scott E.M. (Eds.), Mathematical and Statistical Methods for Sensitivity Analysis. London, John Wiley & Sons, Forthcoming. Homma T. and Saltelli A. 1995. Sensitivity analysis of model output. Performance of the Sobol’ quasi-random sequence generator for the integration of the modified Hora and Iman importance measure. Journal Nuclear Science and Technology 32: 1164– 1173. Homma T. and Saltelli A. 1996. Importance measures in global sensitivity analysis of nonlinear models. Reliability Engineering and System Safety 52: 1–17.

Jansen M.J.W. 1996. Winding stairs sample analysis program WINDINGS 2.0. Technical Report, Private communication. Jansen M.J.W. 1999. Analysis of variance designs for model output. Computer Physics Communications 117: 35–43. Jansen M.J.W., Rossing W.A.H., and Daamen R.A. 1994. Monte Carlo estimation of uncertainty contributions from several independent multivariate sources. In: Gasman J. and van Straten G. (Eds.), Predictability and Nonlinear Modelling in Natural Sciences and Economics, Kluwer Academic Publishers, pp. 334–343. McKay M.D. 1997. Non-parametric variance-based methods for assessing uncertainty importance. Reliability Engineering and System Safety 57: 267–279. McKay M.D., Beckman R.J., and Conover W.J. 1979. A comparison of three methods for selecting values of input variables in the analysis of output from a computer code. Technometrics 21: 239–245. Morris M.D. 1991. Factorial sampling plans for preliminary computational experiments. Technometrics 33: 161–174. Saltelli A. and Bolado R. 1998. An alternative way to compute Fourier Amplitude Sensitivity Test (FAST). Computational Statistics and Data Analysis 26: 445–460. Saltelli A. and Sobol’ I.M. 1995. About the use of rank transformation in sensitivity analysis of model output. Reliability Engineering and System Safety 50: 225–239. Saltelli A., Tarantola S., and Chan K. 1999. A quantitative, model independent method for global sensitivity analysis of model output. Technometrics 41: 39–56. ˘ Saltenis V. and Dzemyda G. 1982. Structure analysis of extremal problems using an approximation of characteristics. Optimal Decision Theory 8: 124–138. Sobol’ I.M. 1967. On the distribution of points in a cube and the approximate evaluation of integrals. USSR Computational Mathematics and Mathematical Physics 7: 86–112. Sobol’ I.M. 1976. Uniformly distributed sequences with an additional uniform property. USSR Computational Mathematics and Mathematical Physics 16: 236–242. Sobol’ I.M. 1993. Sensitivity analysis for nonlinear mathematical models. Mathematical Modeling & Computational Experiment. 1: 407–414. Sobol’ I.M. 1994. A Primer for the Monte Carlo Method. Boca Raton, CRC Press, USA. Sobol’ I.M., Turchaninov V.I., Levitan Y.L., and Shukhman B.V. 1992. Quasirandom sequence generators. Preprint, Keldysh Institute of Applied Mathematics, Russian Academy of Sciences.