AVaN Pack: An Analytical/Numerical Solution for Variance-Based Sensitivity Analysis
AAV A N P
ACK : A N A NALYTICAL /N UMERICAL S OLUTION FOR V ARIANCE -B ASED S ENSITIVITY A NALYSIS
A P
REPRINT
Eduardo Vasconcelos
Academic Departament of Electronic Systems ControlFederal Institute of PernambucoRecife, Pernambuco, Brazil [email protected]
Adriano Souza
Academic Departament of Electronic Systems ControlFederal Institute of PernambucoRecife, Pernambuco, Brazil [email protected]
Kelvin Dias
Center of InformaticFederal University of PernambucoRecife, Pernambuco, Brazil [email protected]
December 25, 2019 A BSTRACT
Sensitivity analysis is an important concept to analyze the influences of parameters in a system,an equation or a collection of data. The methods used for sensitivity analysis are divided intodeterministic and statistical techniques. Generally, deterministic techniques analyze fixed points of amodel whilst stochastic techniques analyze a range of values. Deterministic methods fail in analyzethe entire range of input values and stochastic methods generate outcomes with random errors. In thismanuscript, we are interested in stochastic methods, mainly in variance-based techniques such asVariance and Sobol indices, since this class of techniques is largely used on literature. The objectiveof this manuscript is to present an analytical solution for variance based sensitive analysis. As a resultof this research, two small programs were developed in Javascript named as AVaN Pack ( A nalysis of Va riance through N umerical solution). These programs allow users to find the contribution of eachindividual parameter in any function by means of a mathematical solution, instead of sampling-basedones. K eywords Sensitivity Analysis, Analytical Solution, Analysis of Variance, Parameter Contribution, Sobol Indices
The concept of sensitivity analysis is well known in the literature and it is used in several scientific areas to determinethe impact of parameters in complex systems. There are a large number of techniques on literature, that can bedivided into statistical and deterministic methods [1]. Among these techniques, we can mention the Partial DifferentialApproach, the Rank of Correlation, the Regression method [2] and Variance-Based Methods such as ANOVA (Analysisof Variance), and Sobol indices [3]. The sensitivity analysis concept is widely used in many areas such as finance [4],[5], [6], engineering [7], Biology [8] [9] and many other areas.Deterministic approaches consist in mathematically find the contribution of parameters in a function. So In this concept,the analysis is performed considering specific points over the range of inputs. Partial Derivative analysis [2] is adeterministic technique that determine the contribution of parameters by analyzing the differential of parameter x i over the function f ( x , x , ..., x n ) for specific values of X . The main advantage of this approach is that results do notpresent random errors. The main disadvantage, however, is that it fails in analyze the global scope of complex functionsdue to the fact that this technique’s class analyzes deterministic points (local analysis). The main disadvantage, however, a r X i v : . [ s t a t . O T ] D ec PREPRINT - D
ECEMBER
25, 2019is that due to the fact that this class of technique analyzes deterministic points (local analysis) it fails in analyze theglobal scope of a complex function.Statistic approaches, on the other hand, consist of creating a set of data series through the generation of random values tobe passed as parameters to function f ( X ) . With the data series, it’s possible to generate statistics such as the Correlationof Pearson between each parameter and the function output. The main advantage of this concept is that it can analyzethe range of values of each parameter. The problem with this concept is due to its random nature, since it generatesresults with random error, what can decreases the overall quality of the analysis, as discussed in [1].Besides the problems presented on both concepts of sensitivity analysis, for complex multi-variable functions f ( X ) ,both concepts can fail to analyze certain regions of function f . For deterministic analysis, the person responsible forplanning the analyzes can empirically choose values that cover only some regions of f . For statistic concept, the randomnature of analyzes can bypass some regions of function f , even with a great number of samples.To solve the problems presented, we propose a method to solve statistical sensitivity analysis of continuous functionsthrough analytical solution. In this paper, our proposal focuses on the variance-based class of statistic techniques suchas the indices of Sobol and Variance, since they are largely studied on literature [2][10].As a result of this research, we have developed two scripts, coded in Javascript, that can find the contributions ofparameters in a continuous function through of numerical solution. To this set of scripts, we have given the name ofAVaN Pack ( A nalysis of Va riance through N umerical solution).This work is organized as following: section 2 presents the method used for analytically solve the Variance and theSobol indices; section 3 presents and discusses a practical example commonly used on literature comparing the resultsproduced by the proposed method with those obtained on [2]; section 4 presents the scripts of AVaN Pack; and section 5presents the considerations and future works. Let f ( X ) be a function and X = x , x , x , ..., x n an array of input parameters. S i is the first level sensitivitycontribution of parameter x i over f ( X ) (named as "First Order index") that, in analysis of variance, can be obtained by[2]: S i = V ar ( f ( x i )) V ar ( f ( X )) (1)where V ar ( f ( x i )) corresponds to the variance of function f ( X ) by varying the parameter x i and keeping the otherparameters unchanged, and V ar ( f ( X )) is the variance of function f ( X ) .The method used for calculating the Sobol index So i is different from S i since it considers the variance decompositionof f ( X ) [3]. So, the index So i is defined as: So i = V i V ( Y ) = V ar x i ( E − x i [ f ( X ) | x i ]) V ( Y ) (2)with V ( Y ) = n (cid:88) i V i + n (cid:88) i,j V ij + . . . + V , , ,...,n = V ar x i ( E − x i [ f ( X ) | x i ]) + E x i [ V ar − x i ( f ( X ) | x i )] (3) Equations S i..j = V ar ( f ( x i ...x j )) V ar ( f ( X )) (4)and So i..j = V ar x i...j ( E − x i...j [ f ( X ) | x i...j ]) − V ar x i ( E − x i [ f ( X ) | x i ]) − . . . − V ar x j ( E − x j [ f ( X ) | x j ]) V ( Y ) (5)In general, several authors have used sampling techniques to analyze functions using variance based methods in orderto determine the variances of Equations
1, 2, 4 and 5. To solve the variance of any function analytically, the followingexpression can be used:
V ar ( f ( x )) = E [ f ( x ) ] − E [ f ( x )] (6)2 PREPRINT - D
ECEMBER
25, 2019The expected value E [ f ( x )] of a function can be obtained as: E [ f ( x )] = (cid:90) f ( x ) ϕ ( f ( x )) dx (7)where ϕ ( x ) is the Probability Distribution Function of f ( x ) . It is not trivial and some times impossible to find theProbability Distribution Function of complex equations. So, we can use the mathematical Expected Value of a functiondenoted by: E [ f ( x )] = ( b − a ) − (cid:90) ba f ( x ) dx, b > a (8)The problem of Equation f ( x ) is necessary a range of values. However, asall studies of the sensitivity analysis are performed under a range of values, the solution presented in Equation E [ f ( x ) ] , we can use the following expression: E [ f ( x ) ] = ( b − a ) − (cid:90) ba f ( x ) dx, b > a (9)Based on Equations f ( x ) . But, these equations just calculate the firstdimension variance, and to find the first order contribution we need the n dimension variance of f ( X ) ; in other words,we need to find V ar ( f ( X )) . So, we can make use of the Fubini’s Theorem for any integrable and continuous equationto find this answer: (cid:90) y max y min (cid:90) x max x min f ( x, y ) dxdy = (cid:90) x max x min (cid:90) y max y min f ( x, y ) dydx (10)So, from property presented in Equation (10) we have: E [ f ( X )] = ( n (cid:89) k b k − a k ) − (cid:90) b a (cid:90) b a (cid:90) b a ... (cid:90) b n a n f ( X ) dx n ...dx dx dx (11)Thus, from Equations
6, 8, 9 and 11 we can calculate the indices presented in
Equations
1, 2, 4 and 5.
The problem of calculating
V ar ( f ( x )) and V ( Y ) is the number of integrations that are required, since computationally,these values are generally only obtained through numerical solutions. But, to compute the first order indices for equation f ( X ) we can consider these indices in percentage terms. So, we can write the percentage variance indices as: S % i = S i (cid:32) n (cid:88) k S k (cid:33) − (12)Substituting the indices S i by Equation S % i = V ar ( f ( x i )) V ar ( f ( X )) (cid:32) n (cid:88) k V ar ( f ( x k )) V ar ( f ( X )) (cid:33) − (13)As V ar ( f ( X )) is fixed and behaves as a constant in Equation
13, we have: S % i = V ar ( f ( x i )) V ar ( f ( X )) (cid:32) V ar ( f ( X )) n (cid:88) k V ar ( f ( x k )) (cid:33) − (14)Finally, the percentage index of variance is defined by Equation S % i = V ar ( f ( x i )) (cid:32) n (cid:88) k V ar ( f ( x k )) (cid:33) − (15)Following the same idea for Sobol indices formula, we can derive the percentage contribution of Sobol method asfollow: So % i = So i (cid:32) n (cid:88) k So k (cid:33) − (16)3 PREPRINT - D
ECEMBER
25, 2019 = V ar x i ( E − x i [ f ( X ) | x i ]) V ( Y ) (cid:32) n (cid:88) k V ar x k ( E − x k [ f ( X ) | x k ]) V ( Y ) (cid:33) − (17) = V ar x i ( E − x i [ f ( X ) | x i ]) V ( Y ) (cid:32) V ( Y ) n (cid:88) k V ar x k ( E − x k [ f ( X ) | x k ]) (cid:33) − (18) = V ar x i ( E − x i [ f ( X ) | x i ]) (cid:32) n (cid:88) k V ar x k ( E − x k [ f ( X ) | x k ]) (cid:33) − (19)The Equations
15 and 16 represent the first order contribution of parameters of f ( X ) both for the variance analysisand for the Sobol analysis and can be efficiently calculated, since, we do not have to calculate n integrals to find V ar ( f ( X )) . It’s important to mention that S % i and So % i do not represent the index of variance neither the Sobol indexof f ( X ) . It represents just the individual percentage contribution of parameter x i for variance and Sobol sensitivitytechniques. To demonstrate the use of the analytical solution proposed in this manuscript for both Variance and Sobol indices, wereproduce the experiment presented in [2] by extracting the variance indices of Ishigami Function: Y = sin ( x ) + 7( sin ( y )) + 0 . z ) sin ( x ) (20) To obtain the percentage first order variance of
Equation
20 we need to calculate the Expected Value presented in
Equation
6, considering that the range used in [2] is ( − π/ , π/ for x, y and z. So, the integrations for functions V ar ( f ( x )) , V ar ( f ( y )) and V ar ( f ( z )) are: V ar ( Y ( x )) = E [ Y ( x ) ] − E [ Y ( x )] (21)where: E [ Y ( x ) ] = ( π/ − ( − π/ − (cid:90) π/ − π/ ( sin ( x ) + 7( sin ( y )) + 0 . z ) sin ( x )) dx = 4900 π ( cos (4 y ) − cos (2 y )) − (cid:112) − √ z + 10) + 4 π ( z + 20 z + 2775)4000( π/ − ( − π/ and E [ Y ( x )] = (cid:32) ( π/ − ( − π/ − (cid:90) π/ − π/ ( sin ( x ) + 7( sin ( y )) + 0 . z ) sin ( x )) dx (cid:33) = (cid:18) πsin ( y )5( π/ − ( − π/ (cid:19) For the variances of Y ( y ) and Y ( z ) we have: V ar ( Y ( y )) = E [ Y ( y ) ] − E [ Y ( y )] (22) E [ Y ( y ) ] = ( π/ − ( − π/ − (cid:90) π/ − π/ ( sin ( x ) + 7( sin ( y )) + 0 . z ) sin ( x )) dy = 5 π − (0 . . z + 0 . z +( − . − . z − . z ) cos (2 x ) +(0 . . z ) sin ( x )) PREPRINT - D
ECEMBER
25, 2019 E [ Y ( y )] = (cid:32) ( π/ − ( − π/ − (cid:90) π/ − π/ ( sin ( x ) + 7( sin ( y )) + 0 . z ) sin ( x )) dy (cid:33) = π − (25((0 . z + 0 . sin ( x ) + 0 . and V ar ( Y ( z )) = E [ Y ( z ) ] − E [ Y ( z )] (23) E [ Y ( z ) ] = ( π/ − ( − π/ − (cid:90) π/ − π/ ( sin ( x ) + 7( sin ( y )) + 0 . z ) sin ( x )) dy = ( π/ − (cid:18) π (500000 + π ) sin ( x ) sin ( y )1 . × (cid:19) +( π/ − (cid:18)(cid:2) π/ π / (1 . × ) + ( π / . × ) (cid:3) sin ( x ) + 49 πsin ( x )5 (cid:19) E [ Y ( z )] = (cid:32) ( π/ − ( − π/ − (cid:90) π/ − π/ ( sin ( x ) + 7( sin ( y )) + 0 . z ) sin ( x )) dz (cid:33) = ( π/ sin ( x ) + 7 sin ( y )) + 0 . sin ( x )) As we are interested in first-order influences, we need to define the values assumed to unchanged variables on equations,since in [2] authors did not present these values explicitly. So, In this manuscript, we have fixed the values of x, y and zin 0, since it is the median value of the range − π/ and π/ After solving
Equations
21, 22 and 23, we have obtained the following values:
V ar ( Y ( x )) = 0 . , V ar ( Y ( y )) = 0 . and V ar ( Y ( z )) = 0 . From the variances, we can obtain the percentage contributionusing Equation
15 obtaining the following results: S % x = 44 . , S % y = 55 . and S % z = 0% . Now, we are interested in obtaining the first order Sobol indices of Ishigami Function presented in
Equation
20. Inorder to obtain the Sobol index So % i , we have to follow a simple process: To calculate V ar x i ( E − x i [ f ( X ) | x i ]) we haveto find E [ f ( X )] fixing x i and varying the other parameters; after, we obtain the variance, by varying only x i . Whenusing the sampling method, it is necessary to choose a random value to x i then obtain the Expected Value of f ( X ) byvarying the others x − i parameters; this process must be repeated until we have the desired number of E [ f ( X )] , then,we extract the variance of these Expected Values [2].For analytical solution, the process is more intuitive: we just need to obtain the E − x i [ f ( X ) | x i ] equation and then,obtain its variance over x i using Equations
6, 8, 9 and 11.So, the Sobol percentage indices of
Equation
20 are calculated as follow:
V ar x ( E − x [ Y | x ]) = E x [ E − x [ Y | x ] ] − E x [ E − x [ Y | x ]] (24)with E − x [ Y | x ] = ( π/ − ( − π/ − (cid:90) π/ − π/ (cid:90) π/ − π/ sin ( x ) + 7 sin ( y ) + 0 . z sin ( x ) dydz = 25 π − (0 . sin ( x ) + 0 . and E x [ E − x [ Y | x ] ] = ( π/ − ( − π/ − (cid:90) π/ − π/ (25 π − (0 . sin ( x ) + 0 . dx = 0 . PREPRINT - D
ECEMBER
25, 2019 E x [ E − x [ Y | x ]] = (cid:32) ( π/ − ( − π/ − (cid:90) π/ − π/ π − (0 . sin ( x ) + 0 . dx (cid:33) = 0 . So,
V ar x ( E − x [ Y | x ]) = (0 . − . . . Calculating the other parameters, we have: V ar y ( E − y [ Y | y ]) = E y [ E − y [ Y | y ] ] − E y [ E − y [ Y | y ]] (25) E − y [ Y | y ] = ( π/ − ( − π/ − (cid:90) π/ − π/ (cid:90) π/ − π/ sin ( x ) + 7 sin ( y ) + 0 . z sin ( x ) dxdz = 7 sin ( y ) E y [ E − y [ Y | y ] ] = ( π/ − ( − π/ − (cid:90) π/ − π/ (7 sin ( y )) dy = 0 . E y [ E − y [ Y | y ]] = (cid:32) ( π/ − ( − π/ − (cid:90) π/ − π/ Sin ( y ) dx (cid:33) = 0 . and finally for z , we have: V ar z ( E − z [ Y | z ]) = E z [ E − z [ Y | z ] ] − E z [ E − z [ Y | z ]] (26) E − z [ Y | z ] = ( π/ − ( − π/ − (cid:90) π/ − π/ (cid:90) π/ − π/ sin ( x ) + 7 sin ( y ) + 0 . z sin ( x ) dxdy = 7(4 π − (cid:112) − √ πE z [ E − z [ Y | z ] ] = ( π/ − ( − π/ − (cid:90) π/ − π/ (cid:32) π − (cid:112) − √ π (cid:33) dz = 0 . E z [ E − z [ Y | z ]] = (cid:32) ( π/ − ( − π/ − (cid:90) π/ − π/ π − (cid:112) − √ π dz (cid:33) = 0 . In summary, the
V ar y ( E − y [ Y | y ]) = (0 . − . . and V ar z ( E − z [ Y | z ]) = (0 . − . . Results obtained from Sections 3.1 and 3.2 are similar and may be used to demonstrate the influences of parameters x, y, z . Table
Table
1, the results for variance and Sobol indices analytically computed have a difference ofjust . . This difference can be explained by the precision of integrations values. Comparing with results presentedon [2] we can see the differences between the analytical solutions and the sampling-based solutions, mainly if wecompare the Sobol sampling-based results. These differences can be explained by the nature of the analysis, since, theanalysis presented in [2] has been performed by sampling methods, that by nature have a random error.6 PREPRINT - D
ECEMBER
25, 2019Table 1: Results of ExamplesIndices Analytical Var Var Indices[2] Analitycal Sobol Sobol R[2] Sobol MatLab[2] S % x .
58% 44 .
8% 44 .
59% 41 .
55% 41 . S % y .
42% 55 .
2% 55 .
41% 58 .
46% 58 . S % z
0% 0% 0% 0% 0%
In this section, we present the main features of AVaN Pack . This package is composed by two scripts (Variance Scriptand Sobol Script) that compute the sensitivity contribution of a function through a numerical solution. We decided tosplit the contribution of this manuscript into two programs to increase the code readability, since, these solutions arecomplex. Variance and Sobol scripts are small programs developed in Javascript (JS). These programs can find theindividual percentage contributions of parameters of one equation. We have chosen JS to code these programs for somereasons: • JS is possibly the most used programming language throughout the world; • The pure and old JS code can be executed practically in every smart device, provided it supports a web browser; • Users don’t need of any further structure such as servers or frameworks to execute a JS code, only the textsource file; • JS is a flexible language and codes can easily be modified and extended; • With JS eval function, it’s possible to convert symbolic equations into JS native functions;So, two single-file codes have been developed to compute variance-based sensitivity analysis of any solvable andintegrable equation. To keep these programs as simple as possible, we do not add any JS library or framework such asJQuery or Angular in our programs. With this, we can ensure that this program can execute even in old browsers, oncethey give support to JSON objects.Following, we present the main features of each program.
The variance script has a simple user interface composed by two HTML text areas and one button as showed in
Figure
1. The first text area is used for parameter descriptions, the second one is used for the equation input. Parameters mustbe defined in JSON format. This decision allowed us to keep the program as simple as possible, decreasing the numberof HTML elements and codes that should be necessary to treat these elements. Parameters inserted in the first text areamust be described following 4 JSON parameters: • "param" - represents the name of the parameter; • "min" - contains the minimum value assumed for this parameter; • "max" - contains the maximum value; • "fixed" - The fixed value for these parameters.For equation sin ( x ) + cos ( y ) the text inserted into the first text area must be: { ” param ” : ” x ” , ” min ” : ”1” , ” max ” :”10” , ” f ixed ” : ”5” } & { ” param ” : ” y ” , ” min ” : ”2” , ” max ” : ”5” , ” f ixed ” : ”3” } . Note that we have used "&" toseparate each parameter. The equation must be inputted into the second text area using the JS math format. In thiscase, users must write their equations as if they were programming in Javascript. For example, let’s suppose that anuser wants to analyze equation sin ( x ) + e cos ( y + x ) . Hence, the text that must be inputted into the second text area is: M ath.sin ( x ) + M ath.exp ( M ath.cos ( y + x )) . AVaN Pack is a GPL-3.0 licensed software and is available in https://github.com/dr-eduardovasconcelos/AVaNPack PREPRINT - D
ECEMBER
25, 2019Figure 1: Graphical User Interface of Variance ScriptTo solve the analytical variance, the variance script makes use of numerical integration to find the expected values. TheJS integration algorithm used is presented as follows:
Algorithm 1:
Numerical Integration function numericalIntegration( f, a, b ) { var delta = 0 . ∗ ( b − a ) ; var area = 0 ; for ( var i = a ; i < b − delta ; i + = delta ) { area + = ( delta ) ∗ (( f ( i ) + f ( i + delta )) / } return area ; }Function numericalIntegration () uses the trapezoidal method of integration to find the area of a function between a and b . This method of integration consists in multiplying the height ( b − a ) with the average of two bases f ( a ) and f ( b ) . The variable delta stores the partial heights of each trapezium that composes the integration. The value . represents the granularity of the process and, with this value, we guarantee the existence of trapeziums perintegration. We assume that the numerical integration imposes an error to results, but, we argue that different from thesampling-based solution, the error produced is not random, actually, this is a precision error.With the method of integration defined, it is necessary to define the method to compose the JS functions. As we cansee in numericalIntegration() function, parameter f is a single-parameter function, and contains the equation passedby the user. But, as a sensitivity analysis is made over an equation that has a set of parameters, define one singleparameter function is unfeasible. So, to solve this problem, we used the flexibility of JS to create a set of texts thatcontain the partial functions for integration. After, we make use of code eval() in order to turn these partial functionsinto JS functions. The definition of these function has been split into two phases: in the first, we need to create a8 PREPRINT - D
ECEMBER
25, 2019multi-parameter function to execute the user equation; in the second phase, we define single parameter functions thatwill be used for numerical integration. The code used to create the first function is showed as following:
Algorithm 2:
String Equation Contruction var equationInJSF mt = ” f unction equation (” ; for ( var i = 0; i < parameters.length ; i + + ){ equationInJSF mt + = parameter [ i ] .param ; equationInJSF mt + = ( i == parameters.length −
1) ? ”) { return ” + equationString + ”; } ” : ” , ” ;} eval ( equationInJSF mt ); The variable parameters contains an array of JS objects, each of these objects contains the parameters defined into thefirst text area (See
Figure equationString contains the equation defined into the second text area. At theend of the definition of equation() , we will have one function such as " function equation(x,y){returns Math. sin(x)+ Math.cos(y);};" for the equation
Math.sin(x)+Math.cos(y) . After all, the function eval() converts the text into a JSfunction.To find the variance of f ( x ) we need to calculate f ( x ) , thus we have to create a second function to be used to compute E [ f ( x ) ] . This second function is defined as the following: Algorithm 3:
Creating Equation Moment var momentEquationString = equationInJSF mt.replace (” equation ” , ” equationM oment ”) .replace (” return ” , ” returnM ath.pow (”) .replace (”; ” , ” , eval ( momentEquationString ); This code basically uses the replace function to change equation into equationMoment . As a result, the String momentEquationString will contain " function equationMoment(x,y) {return Math.pow(Math.sin(x)+Math.cos(y),2);}; ".Functions equation and equationMoment are the first equations used in our solution. Now, we have to define thefunctions that will be directly used by numericalIntegration() . So the definition of these single-parameter functions isdefined as follow:
Algorithm 4:
Definition of _fun and _funMoment var _ f un ; var _ f unM oment ; var transientF unctionEqString = ” _ f un = f unction ( _ param ) { return equation (”; for ( var j = 0; j < parameters.length ; j + + ) { transientF unctionEqString + = ( j == i ) ? ” _ param ” : parameters [ j ] .f ixed + ””; transientF unctionEqString + = ( j ! = parameters.length −
1) ? ” , ” : ”); } ”; } transientF unctionM oString = transientF unctionEqString .replace (” equation ” , ” equationM oment ”) .replace (” _ f un ” , ” _ f unM oment ”); eval ( transientF unctionEqString ); eval ( transientF unctionM oString ); This code repeats for each parameter on analysis. Functions _fun and _funMoment contain JS native functions and willbe used to compute the expected values E [ f ( x )] and E [ f ( x ) ] . The idea of the above code is create a function thatreceives one parameter and changes the others for their respective fixed values. For the function function equation(x,y){... the resultant function for the first iteration is " _fun = function(_param){return equation(_param, 3);}; " if the fixedvalue of y is 3. 9 PREPRINT - D
ECEMBER
25, 2019Once the functions have been defined the rest of the process consists in integrate each function using the following code:
Algorithm 5:
Calculation of Variance Indices var variances = new Array ( parameters.length ); for ( var i = 0; i < parameters.length ; i + + ){ var min = parseF loat ( parameters [ i ] .min ); var max = parseF loat ( parameters [ i ] .max ); var squareEsperance = numericalIntegration ( _ f un, min, max ) / ( max − min ); squareEsperance = M ath.pow ( squareEsperance, var moment = numericalIntegration ( _ f unM oment, min, max ) / ( max − min ); variances [ i ] = moment − squareEsperance ; }After this process, we obtain the analytical variances of parameters, and by using the Equation
15, we can obtain thefirst order percentage influences.
The user interface of Sobol script is similar to that showed in
Figure E − x i [ f ( X ) | x i ] . This expected value is computed by function nMinusXiAverageIntegration() that is showed asfollow: Algorithm 6:
First part of nMinusXiAverageIntegration function function nMinusXiAverageIntegration( equation, xi, X, delta_base ) { var delta _ X = newArray ( X.length ); var delta _ inc X = newArray ( X.length ); var baseP roduct = 1 . var aM inusbP rod = 1 . for ( var i = 0; i < X.length ; i + + ) { delta _ X [ i ] = delta _ base ∗ ( X [ i ] .max − X [ i ] .min ); delta _ inc _ X [ i ] = X [ i ] .min ; baseP roduct ∗ = delta _ X [ i ]; aM inusbP rod ∗ = X [ i ] .max − X [ i ] .min ; delta _ X [ i ] = delta _ base ∗ ( X [ i ] .max − X [ i ] .min ); delta _ inc _ X [ i ] = X [ i ] .min ; . . .Function nMinusXiAverageIntegration() receives 4 arguments: The first is a literal that contains the equation that mustbe analyzed; the second is an object that receives the name and the value of x i ; X is an array of objects representingthe other parameters of equation, and delta_base is a float number that represents the granularity of our integration.The variable delta_X is an array that contains the heights of parameters of the array X and delta_inc_X is an array ofnumbers used for controlling the process of integration.The following code shows the process of defining the literal functions used on analyzes. Algorithm 7:
Second part of nMinusXiAverageIntegration function var transBaseF un = ” _ baseF un = f unction (” + xi.param + ” , ” ; var _ basef un ; for ( var i = 0; i < X.length ; i + + ) { transBaseF un = transBaseF un + X [ i ] .param + (( i ! = X.length −
1) ? ” , ” :”) { return ” + equation + ”; } ; ”); } eval ( transBaseF un ); . . .Variable transBaseFun has the base function of the integration, which computes the equation passed as parameters. Thefirst argument of _basefun is x i , the other arguments are inserted into the function according to the sequence of array X .10 PREPRINT - D
ECEMBER
25, 2019If the equation passed as argument is x + y + z then the resultant funtion is " _basefun(y,x,z){return x+y+x;}; ". Notethat, in this case, the variable y is the x i in this process.Function _basefun contains only the analyzed equation. To make the multiple integrations, we need another functionthat makes the sums necessary for the trapezoidal integration. This second function is assembled as follow: Algorithm 8:
Third part of nMinusXiAverageIntegration function var transIntF un = ” _ intF un = f unction ( delta _ inc _ X ) { return _ baseF un (” + xi.value + ” , ”; var _ intF un ; for ( var i = 0; i < M ath.pow (2 , X.length ); i + + ) { for ( var j = 0; j < X.length ; j + + ) { transIntF un + = ” delta _ inc _ X [” + j + ”]”; if ( i & M ath.pow (2 , j ))! = 0) { transIntF un + = ” + ” + delta _ X [ j ]; } transIntF un + = ( j == X.length −
1) ? ”) ” : ” , ”; } transIntF un + = ( i == M ath.pow (2 , X.length ) −
1) ? ”; } ; ” : ” + _ baseF un (” + xi.value + ” , ”; } eval ( transIntF un ); . . .The base of function _intFun represents the generalization of trapezoidal method for integration of multiple variables.To explain the resultant format of _intFun , let’s consider that _basefun contains 3 arguments x, y and z, and yis the analyzed parameter. Let’s consider also that the fixed value of y is 3 and the other parameters are beingintegrated with the intervals a, b and c, d. So, the general form of this trapezoidal integration is " (b-a)*(d-c) *(_basefun(3,a,c)+_basefun(3,a+b,c)+_basefun(3,a,c+d)+ _basefun(3,a+b,c+d))/4 ". Function intFun computes thesums of _basefun and the product among (b-a) and (d-c) is represented by the variable baseProduct that is introducedinto the next piece of code. Basically, in function _intFun the intervals of each parameter are given by array delta_X .Once we have the two functions _basefun and _intFun , the next step is to compute the expected value E − x i [ f ( X ) | x i ] through the following piece of code: Algorithm 9:
Last part of nMinusXiAverageIntegration function var partialArea = 0 . for ( var i = 0; i < M ath.pow (1 . /delta _ base, X.length ); i + + ) { partialArea + = baseP roduct ∗ ( _ intF un ( delta _ inc _ X ) /M ath.pow (2 , X.length )); delta _ inc _ X [0] = ( delta _ inc _ X [0] + delta _ X [0]); if ( delta _ inc _ X [0] > X [0] .max − delta _ X [0] ) { delta _ inc _ X [0] = X [0] .min ; for ( var j = 1; j < X.length ; j + + ) { if ( delta _ inc _ X [ j −
1] == X [ j − .min ) { delta _ inc _ X [ j ] = ( delta _ inc _ X [ j ] + delta _ X [ j ]); if ( delta _ inc _ X [ j ] > X [ j ] .max − delta _ X [ j ] ){ delta _ inc _ X [ j ] = X [ j ] .min ; }}}}} return partialArea/aM inusbP rod ; Since for computing an n-variable numerical integration we need of n loops, in Sobol Script, to adapt the algorithmfor any number of parameters, we have decided to create a bigger single loop. This loop will have (1 /delta _ base ) n iterations, that is the same complexity of having n loops. After calculating the sum of trapeziums, we need to updatethe values of delta_inc_X array. As delta_inc_X[i] has a value between X[i].min and
X[i].max , this value must beincremented until it reaches the maximum value of parameter
X[i] . Once delta_inc_X[i] reaches its maximum value,the program resets it and increment delta_inc_X[i+1] . So, in the above piece of code, delta_inc_X[0] is incrementedin every iteration, and always which this variable reaches its maximum value, the code scrolls through the array toupdate the other parameters. Then, after computing the main sum, the function nMinusXiAverageIntegration returns theratio between the sum result and the product of maximum and minimum values of parameters, this ratio represents themulti-variable expected value.Once we have computed E − x i [ f ( X ) | x i ] it’s necessary to compute V ar x i ( E − x i [ f ( X ) | x i ]) . The process to compute V ar x i for symbolic integration is known, we need to calculate the resultant equation E − x i [ f ( X ) | x i ] and then integrate11 PREPRINT - D
ECEMBER
25, 2019it to obtain the variance. However, for numerical integration it cannot be done, since it’s not possible to integrate amulti-variable function in parts. Then, the following piece of code shows the solution.
Algorithm 10:
Sobol Indices Calculation var delta _ base = 0 . for ( var i = 0; i < parameters.length ; i + + ){ var minV al = parseF loat ( parameters [ i ] .min ); var maxV al = parseF loat ( parameters [ i ] .max ); var delta _ xi = ( maxV al − M inV al ) ∗ delta _ base ; var areaM omento = 0 . var areaSquare = 0 . for ( var j = minV al ; j < maxV al ; j + = delta _ xi ) { xi.value = j + delta _ xi/ . var partialInteg = nM inusXiAverageIntegration ( equationString, xi, X, delta _ base ); areaM omento + = delta _ xi ∗ M ath.pow ( partialInteg, areaSquare + = delta _ xi ∗ partialInteg ; } V [ i ] = areaM omento/ ( maxV al − minV al ) − M ath.pow ( areaSquare/ ( maxV al − minV al ) , sumV + = V [ i ]; }In this solution, we basically scroll through the minimum and maximum values in delta_xi steps, compute the expectedvalue of − x i parameters and then, sum these expected values to obtain E x i [ E − x i [ f ( x ) | x i ] ] and ( E x i [ E − x i [ f ( x ) | x i ]]) .Note that the value of x i is defined as j + delta_xi/2.0 , it’s because, for this integration, we have used the rectangle rule.This measure allowed us to reduce the complexity of our solution, since, to compute the trapezoidal rule we need to call nMinusXiAverageIntegration() twice for each value of x i . The code " areaMomento += delta_xi * Math.pow(partialInteg,2) " represents the partial calculation of E x i [ E − x i [ f ( x ) | x i ] ] as presented in Equation
9. Finally, all variances arestored into array V . From V , it’s possible to compute the Sobol percentage indices of user equation. The complexity ofSobol script is O ( n ) = n × ( delta _ base − n ) for delta _ base < and n the number of parameters. The objective of this section is to compare the outcomes of Variance and Sobol scripts with results presented in
Table
1. The object of analysis is the Ishigami Function presented on
Equation
20. The JS format of Ishigami Function ispresented as following:
M ath.sin ( x ) + 7 ∗ M ath.pow ( M ath.sin ( y ) ,
2) + 0 . ∗ ( M ath.pow ( z, ∗ M ath.sin ( x ) .The above text is set into the second text area of both programs. The parameters in JSON format are described asfollowing:Variance Script{"param":"x","min":"-0.31415926535897932384626433", "max":"0.31415926535897932384626433","fixed":"0.0"} &{"param":"y","min":"-0.31415926535897932384626433", "max":"0.31415926535897932384626433","fixed":"0.0"} &{"param":"z","min":"-0.31415926535897932384626433", "max":"0.31415926535897932384626433","fixed":"0.0"}Sobol Script{"param":"x","min":"-0.31415926535897932384626433", "max":"0.31415926535897932384626433"} &{"param":"y","min":"-0.31415926535897932384626433", "max":"0.31415926535897932384626433"} &{"param":"z","min":"-0.31415926535897932384626433", "max":"0.31415926535897932384626433"}The values of delta and delta_base for both Variance Script and Sobol Script, has been set as and ,respectively. The reason for setting delta_base with a greater value then delta is the complexity of Sobol Script code. Table
Table delta and delta_base doesn’t mean that these results will repeat for any experiment performed on AVaN Pack.So, whenever questions remain about the outcomes, it will be necessary to decrease the values of delta or delta_base .12 PREPRINT - D
ECEMBER
25, 2019Table 2: Program outcomesIndices Variance Script Analytical Var Sobol Script Analitycal Sobol S % x .
58% 44 .
58% 44 .
59% 44 . S % y .
42% 55 .
42% 55 .
41% 55 . S % z
0% 0% 0% 0%
References [1] Robert G. Hendrickson. A survey of sensitivity analysis methodology. national bureau of standards department ofcommerce washington. 1984.[2] ALEXANDRE ALLARD and NICOLAS FISCHER.
SENSITIVITY ANALYSIS IN METROLOGY: STUDY ANDCOMPARISON ON DIFFERENT INDICES FOR MEASUREMENT UNCERTAINTY , pages 1–6. 2009.[3] G. E. B. Archer, A. Saltelli, and I. M. Sobol. Sensitivity measures,anova-like techniques and the use of bootstrap.
Journal of Statistical Computation and Simulation , 58(2):99–120, 1997.[4] Xiaodong Xu and John R. Birge. Joint production and financing decisions: Modeling and analysis. 2005.[5] E. Borgonovo, S. Gatti, and L. Peccati. What drives value creation in investment projects? an application ofsensitivity analysis to project finance transactions.
European Journal of Operational Research , 205(1):227 – 236,2010.[6] Jae-Il Yoo, Eul-Bum Lee, and Jin-Woo Choi. Balancing project financing and mezzanine project financing withoption value to mitigate sponsor’s risks for overseas investment projects.
Sustainability , 10(5), 2018.[7] Alessandra Oppio, Stefano Corsi, Francesca Torrieri, and Sergio Mattia.
Infrastructure Development andTerritorial Vulnerability. The Role of Composite Indicators for Addressing Siting Decisions , pages 277–290.Springer International Publishing, Cham, 2017.[8] N G Cogan, Atanaska Dobreva, and Ralf Paus. Analysing the dynamics of a model for alopecia areata asan autoimmune disorder of hair follicle cycling.
Mathematical Medicine and Biology: A Journal of the IMA ,35(3):387–407, 08 2017.[9] Carl M. Gay, Pan Tong, Robert J. Cardnell, Xiao Su, Nene N. Kalu, Upasana Banerjee, Rasha O. Bara, Faye M.Johnson, John V. Heymach, Jing Wang, and Lauren A. Byers. Abstract 1560: Differential sensitivity analysis forresistant malignancies (disarm), a novel approach for drug screen analysis, identifies common candidate drugsacross platinum-resistant cancer types.
Cancer Research , 77(13 Supplement):1560–1560, 2017.[10] Bertrand Iooss and Paul Lemaître.