Home

NISP Toolbox Manual - Forge

1. 3 a o o 3 a Figure 5 2 Monte Carlo Sampling Normal random variable Empirical histogram of X2 0 0 2 0 2 1 2 2 2 3 24 2 5 2 6 2 7 2 8 2 9 3 0 Figure 5 3 Monte Carlo Sampling Uniform random variable 30 In order to install this module we can run the atomsInstall function as in the following script atomsInstall stixbox The following script compares the empirical and theoretical distributions scf histplot 50 sampling 1 xtitle Empirical histogramuofuX1 x linspace 15 15 1000 y dnorm x 1 3 plot x y r legend Empirical Exact The previous script produces the figure 5 4 E Empirical Exact TN MIA 0 05 4 a 3 n od o 3 a Figure 5 4 Monte Carlo Sampling Histogram and exact distribution functions for the first variable The following script performs the same comparison for the second variable scf histplot 50 sampling 2 xtitle Empirical histogram of X2 x linspace 2 3 1000 y ones 1000 1 plot x y r The previous script produces the figure 5 5 5 2 2 A Monte
2. passed ref created Summary tests 4 100 passed 0 o failed 0 o skipped 0 o length 3 84 sec Tests ending the 2009 11 18 at 12 47 48 end verbatim 15 Chapter 3 Configuration functions In this section we present functions which allow to configure the NISP toolbox The nisp_ functions allows to configure the global behaviour of the toolbox These func tions allows to startup and shutdown the toolbox and initialize the seed of the random number generator They are presented in the figure 3 1 nisp_startup Starts up the NISP toolbox nisp_shutdown Shuts down the NISP toolbox level nisp_verboselevelget Returns the current verbose level nisp_verboselevelset level Sets the value of the verbose level nisp_initseed seed Sets the seed of the uniform random number generator nisp_destroyall Destroy all current objects nisp_getpath Returns the path to the current module nisp_printall Prints all current objects Figure 3 1 Outline of the configuration methods The user has no need to explicitely call the nisp_startup and nisp_shutdown func tions Indeed these functions are called automatically by the etc NISP start and etc NISP quit scripts located in the toolbox directory structure The nisp_initseed seed is especially useful when we want to have reproductible re sults It allows to set the seed of the generator at a particular value so that the sequence of
3. f i cos x sin z da 7 115 T TT The previous equation leads to T cos z sin z da m 7 116 T which implies n cos z sin x dx a 7 117 T We finally plug the equation 7 117 into 7 110 and get the equation 7 108 7 11 2 Expectation By assumption the three random variables X X and X3 are independent so that the joint distribution function is the product of the three distribution functions fi i e maan 22 23 g1 21 92 12 ga 23 7 118 By definition the expectation of the random variable sin X4 is E sin X f f H sin t1 91 23 T1 22 73 dzydzodzg 7 119 H i d f sin 21 gi 21 g gt 22 ga 2a dzydzydzs 7 120 1 re x sin 11 dx1 7 121 0 7 122 69 By definition the expectation of the random variable sin X2 is E sin X d f f sin z3 gi 3 21 12 23 d1 dx2dx3 J sin a2 g2 r2 dx2 The equality 7 100 then implies that l d ede Be o 2m 1 y By definition the expectation of the random variable XJ is T T T 4 J d d T391 2 3 T1 La 23 daydxod2g TY TR YT E X3 vimine f x393 x3 dx3 1 T 7 123 7 124 7 125 7 126 7 127 7 128 7 129 7 130 7 131 7 132 7 133 7 134 7 135 We are now going to use the expectations 7 122 7 127 and 7 135 in order to compute the expectation of the output Y The model 7 97 is a sum of functions Since the expectation of a sum of two rand
4. Name Parameter 1 a Parameter 2 b Conditions Normale Ii 0 g 1 Bion Uniforme a 0 b a b Exponentielle 1 LogNormale w 0 1 o 1 0 nc LogUniforme a 0 1 b 1 0 a b gt 0 a lt b Figure 4 2 Default parameters for distributions functions 4 1 2 Parameters of the Log normal distribution For the LogNormale law the distribution function is usually defined by the expected value y and the standard deviation o of the underlying Normal random variable But when we create a LogNormale randvar the parameters to pass to the constructor are the expected value of the LogNormal random variable E X and the standard deviation of the underlying Normale random variable o The expected value and the variance of the Log Normal law are given by E X exp 50 4 1 V X exp o 1 exp 2u 07 4 2 It is possible to invert these formulas in the situation where the given parameters are the expected value and the variance of the Log Normal random variable We can invert completely the previous equations and get u ELCH 5n rk 4 3 ei 1 G 4 4 In particular the expected value u of with the Normal random variable satisfies the equation p In E X oi 4 5 4 1 3 Uniform random number generation In this section we present the generation of uniform random numbers The goal of this section is to warn users about a current limitation of the
5. 01413 H103 We divide the previous inequality by V Y and get ei mob SEN OM Therefore the sum of the first order sensitivity indices satisfies the inequality Si 5 lt 1 62 7 55 7 56 7 57 7 58 7 59 7 60 7 61 7 62 7 63 7 64 7 65 7 66 7 67 Hence in this example one part of the variance V Y cannot be explained neither by Xi alone or by X alone because it is caused by the interactions between X and X We define by 512 the sensitivity index associated with the group of variables X X5 as 010 vw 1 1 1 S9 7 68 The following Scilab script performs the sensitivity analysis on the previous example We consider two normal variables where the first variable has mean 1 5 and standard deviation 0 5 while the second variable has mean 3 5 and standard deviation 2 5 function y Exemple x y 1 x 1 x 2 endfunction First variable Normal mui 1 5 Sigmal 0 5 Second variable Normal mu2 3 5 Sigma2 2 5 1 Two stochastic normalized variables vxi randvar new Normale vx2 randvar new Normale 2 A collection of stoch variables srvx setrandvar_new setrandvar_addrandvar srvx vxl setrandvar_addrandvar srvx vx2 3 Two uncertain parameters vul randvar_new Normale mul sigmal vu2 randvar_new Normale mu2 sigma2 4 A collection of uncertai
6. setrandvar_new srv setrandvar_new n srv setrandvar_new file Methods setrandvar_setsample setrandvar_setsample setrandvar_setsample srv k value setrandvar setsample srv value setrandvar save srv file np setrandvar getsize srv sample setrandvar getsample srv k sample setrandvar getsample srv k sample setrandvar getsample srv setrandvar getlog srv nx setrandvar getdimension srv setrandvar freememory srv setrandvar buildsample srv srv2 setrandvar buildsample srv name np setrandvar buildsample srv name np ne setrandvar addrandvar srv rv Destructor setrandvar destroy srv Static methods tokenmatrix setrandvar tokens nb setrandvar size srv name np srv k i value i Figure 5 1 Outline of the methods of the setrandvar class 28 variable is associated with a Uniform distribution function The simulation is based on 1000 experiments The function nisp_initseed is used to set the value of the seed to zero so that the re sults can be reproduced The setrandvar_new function is used to create a new set of ran dom variables Then we create two new random variables with the randvar_new function These two variables are added to the set with the setrandvar_addrandvar function The setrandvar_buildsample allows to build the design of experiments which can be retrieved as matri
7. 1 Then we create a Petras DOE for the stochastic collection srvx and transform it into a DOE for the uncertain parameters srvu nx setrandvar_getdimension srvu srvx setrandvar_new nx degre 9 setrandvar buildsample srvx Petras degre AT setrandvar_buildsample srvu srvx We use the polychaos_new function and create the new polynomial pc with 3 inputs and 1 output noutput 1 pc polychaos_new srvx noutput The next step allows to perform the simulations associated with the DOE prescribed by the collection srvu Here we must perform np 751 experiments np setrandvar getsize srvu polychaos setsizetarget pc np This is slow for k 1 np inputdata setrandvar getsample srvu k outputdata Exemple inputdata polychaos settarget pc k outputdata end The previous loop works but is slow when np is large Instead the following script is fast because it uses vectorization This is fast inputdata setrandvar getsample srvu outputdata Exemple inputdata polychaos settarget pc outputdata We can now compute the polynomial expansion by integration polychaos setdegree pc degre polychaos computeexp pc srvx Integration Everything is now ready so that we can do the sensitivy analysis as in the following script average polychaos getmean pc var polychaos getvariance pc mprintf Meanuuuuuuuu uibf n average mprintf Varianceyyuu sf n var m
8. 1 750000 Variance 1 000000 Indice de sensibilite du 1er ordre Variable X1 0 765625 Variable X2 0 187500 Indice de sensibilite Totale Variable X1 0 812500 Variable X2 0 234375 In order to free the memory required for the computation it is necessary to delete all the objects created so far polychaos destroy pc randvar destroy vu1 randvar destroy vu2 randvar destroy vx1 randvar destroy vx2 setrandvar_destroy srvu setrandvar_destroy srvx Prior to destroying the objects we can inquire a little more about the density of the output of the chaos polynomial In the following script we create a Latin Hypercube Sampling made of 10 000 points Then get the output of the polynomial on these inputs and plot the histogram of the output polychaos_buildsample pc Lhs 10000 0 sample_output polychaos_getsample pc scf histplot 50 sample output xtitle Product function Empirical Histogram X P X The previous script produces the figure 6 3 We may explore the following topics e Perform the same analysis where the variable X5 is a normal variable with mean 2 and standard deviation 2 e Check that the development in polynomial chaos on a Hermite Hermite basis does not allow to get exact results See that the convergence can be obtained by increasing the degree e Check that the development on a basis Hermite Legendre allows to get exact results with degree 2 6 2 2 The Ishigami test c
9. E 5 gt ist i d d 4d 2 8 S eg 131 2 i fr a_i i i i Ml 0 2 0 3 0 4 0 5 0 6 0 7 0 8 0 Variable 1 Uniforme 0 0 1 0 Figure 5 14 Petras sampling Two uniform variables in the interval 1 1 39 Quadrature Sampling 1 0 TT 3 T T T T T T FF rf T Er e e E E H E e o ee ese o e D e e o e e ee 0 9 ho e e e e e e e e o e ee ho D D D CJ D D D D D D D D ee pP EE e o e ee Fe e o e A Q9 boe o e e e o e e ee e S 06 g po o o H D D D D D D D D D D oe Se ar ew dm Yen tas Ue Ce Ce Cu een leg S pee e o eee 047 2 S ooo e e o o e ee S 7 0 30 0 0 e e c H e e e ee eens e e a e e ee 027 boss e o e D D e o e o ee ke s 0 e e o o c ee 0 1 ee e D D D LU D D D D D D D D D eee o D D e o o o o Mes C ii oo I 0 0 0 1 0 2 0 3 04 0 5 0 6 0 7 0 8 0 9 Variable 1 Uniforme 0 0 1 0 Figure 5 15 Quadrature sampling Two uniform variables in the interval 1 1 Sobol Sampling Variable 2 Uniforme 1 0 1 0 Variable 1 Uniforme 1 0 1 0 Figure 5 16 Sobol sampling Two uniform variables in the interval 1 1 40 Chapter 6 The polychaos class 6 1 Introduction The polychaos class allows to manage a polynomial chaos expansion The coefficients of the
10. This set may be read from a data file The object is linked with a collection of random variables Then the coefficients of the polynomial can be computed by integration quadra ture Once done the mean the variance and the Sobol indices can be directly computed from the coefficients The figure 1 1 presents the NISP methodology The process requires that the user has a numerical solver which has the form Y f X where X are input uncertain parameters and Y are output random variables The method is based on the following steps e We begin by defining normalized random variables For example we may use a random variables in the interval 0 1 or a Normal random variable with mean 0 and variance 1 This choice allows to define the basis for the polynomial chaos denoted by Yx x gt o Depending on the type of random variable the polynomials Yx x gt o are based on Hermite Legendre or Laguerre polynomials e We can now define a Design Of Experiments DOE and with random variable transforma tions rules we get the physical uncertain parameters X Several types of DOE are available Monte Carlo Latin Hypercube Sampling etc If N experiments are required the DOE define the collection of normalized random variables amp 1 w Transformation rules allows to compute the uncertain parameters X 1 v which are the input of the numerical solver f e We can now perform the simulations that is compute the collection of outputs
11. Yi 1 w where Y f X e The variables Y are then projected on the polynomial basis and the coefficients y are computed by integration or regression Random Uncertain Numerical Spectral Variable Parameter Solver Projection 6 X Y f X Y 2 y w 5 Figure 1 1 The NISP methodology 13 The NISP module The NISP toolbox is available under the following operating systems e Linux 32 bits e Linux 64 bits e Windows 32 bits e Mac OS X The following list presents the features provided by the NISP toolbox 4 e Manage various types of random variables uniform normal exponential log normal Generate random numbers from a given random variable Transform an outcome from a given random variable into another e Manage various Design of Experiments for sets of random variables Monte Carlo Sobol Latin Hypercube Sampling various samplings based on Smolyak designs Manage polynomial chaos expansion and get specific outputs including mean variance quantile correlation etc e Generate the C source code which computes the output of the polynomial chaos expansion This User s Manual completes the online help provided with the toolbox but does not replace it The goal of this document is to provide both a global overview of the toolbox and to give some details about its implementation
12. randvar_size Figure 4 3 Outline of the methods of the randvar class 4 2 2 The Oriented Object system In this section we present the token system which allows to emulate an oriented object program ming with Scilab We also present the naming convention we used to create the names of the functions The randvar class provides the following functions e The constructor function randvar_new allows to create a new random variable and returns a token rv e The method randvar_getvalue takes the token rv as its first argument In fact all methods takes as their first argument the object on which they apply 19 e The destructor randvar_destroy allows to delete the current object from the memory of the library e The static methods randvar_tokens and randvar_size allows to quiery the current object which are in use More specifically the randvar_size function returns the number of current randvar objects and the randvar_tokens returns the list of current randvar objects In the following Scilab sessions we present these ideas with practical uses of the toolbox Assume that we start Scilab and that the toolbox is automatically loaded At startup there are no objects so that the randvar_size function returns 0 and the randvar_tokens function returns an empty matrix gt nb randvar_size nb 0 gt tokenmatrix randvar tokens tokenmatrix We now create 3 new random variables based on the Unif
13. s2 sin x 2 x34 x 3 4 72 y 1 st a s2 2 b x34 s1 endfunction function C nisp_cov x y Returns the empirical covariance matrix of x and y x x y y n size x x x mean x y y mean y C 1 1 x x n 1 C 1 2 x y n 1 C 2 1 C 1 2 C 2 2 y y n 1 endfunction function s sensitivityindex ya yc Returns the sensitivity index associated with experiments ya and yc C nisp_cov ya yc s C 1 2 st deviation ya st deviation yc endfunction Create the uncertain parameters rvui randvar_new Uniforme pi pi rvu2 randvar_new Uniforme pi pi rvu3 randvar_new Uniforme pi pi srvu setrandvar_new setrandvar_addrandvar srvu rvul setrandvar_addrandvar srvu rvu2 setrandvar_addrandvar srvu rvu3 The number of uncertain parameters is nx setrandvar_getdimension srvu np 10000 Create a first sampling A Setrandvar buildsample srvu Lhs np A setrandvar getsample srvu Create a first sampling B setrandvar_buildsample srvu Lhs np B setrandvar getsample srvu Perform the experiments in A ya ishigami A Compute the first order sensitivity index for X1 C B C 1 np 1 A 1 np 1 yc ishigami C si sensitivityindex ya yc mprintf Siw u4fu expectedu u4f Nn s1 0 3139 Compute the first order sensitivity index for X2 73 C B C 1 np 2 A 1 np 2
14. srvx Quadrature degre setrandvar_buildsample srvu srvx The next steps will be to create the polynomial and actually perform the DOE But before doing this we can take a look at the DOE associated with the stochastic and uncertain collection of random variables We can use the setrandvar_getsample twice and get the following output gt setrandvar_getsample srvx ans 1 7320508 0 1127017 1 7320508 0 5 1 7320508 0 8872983 0 0 1127017 0 0 5 0 0 8872983 1 7320508 0 1127017 1 7320508 0 5 1 7320508 0 8872983 gt setrandvar_getsample srvu ans 0 1339746 1 1690525 0 1339746 1 75 0 1339746 2 3309475 1 1 1690525 1 1 75 1 2 3309475 1 8660254 1 1690525 1 8660254 1 75 1 8660254 2 3309475 These two matrices are a 9 x 2 matrices where each line represents an experiment and each column represents an input random variable The stochastic normalized srvx DOE has been created first then the srvu has been deduced from srvx based on random variable transformations We now use the polychaos_new function and create a new polynomial pc The number of input variables corresponds to the number of variables in the stochastic collection srvx that is 2 and the number of output variables is given as the input argument ny In this particular case the number of experiments to perform is equal to np 9 as returned by the setrandvar_getsize function This parameter is passed to the polynomial pc with the polychaos_set
15. uniform pseudo random numbers is deterministic When the toolbox is started up the seed is automatically set to 0 which allows to get the same results from session to session 16 Chapter 4 The randvar class In this section we present the randvar class which allows to define a random variable and to generate random numbers from a given distribution function 4 1 The distribution functions In this section we present the distribution functions provided by the randvar class We especially present the Log normal distribution function 4 1 1 Overview The table 4 1 gives the list of distribution functions which are available with the randvar class 3 Each distribution functions have zero one or two parameters One random variable can be specified by giving explicitely its parameters or by using default parameters The parameters for all distribution function are presented in the figure 4 2 which also presents the conditions which must be satisfied by the parameters Name F a EX V X Normale VE exp 1 ep u c Uniforme ie i E i bta Bes Exponentielle pem FAR y S S g i LogNormale ox am O P 2 ar p eed exp u 30 exp o 1 exp 2u e 0 fx lt 0 LogUniforme TI e E i i MORBO soi E x Figure 4 1 Distributions functions of the randvar class The expected value is denoted by E X and the variance is denoted by V X 17
16. yc ishigami C s2 sensitivityindex ya yc mprintf S2 4fuCexpected f n s2 0 4424 Compute the first order sensitivity index for X3 C Bj C 1 np 3 A 1 np 3 yc ishigami C s3 sensitivityindex ya yc mprintf S3 f expected y ff in s3 0 0 Compute the first order sensitivity index for X1 X3 C A C 1 np 2 B 1 np 2 yc ishigami C s13 sensitivityindex ya yc mprintf S13 14fu expected 04f n s13 0 5576 77 Clean up randvar_destroy rvul randvar_destroy rvu2 randvar_destroy rvu3 Setrandvar destroy srvu The previous script produces the following output S1 0 317135 expected 0 313900 S2 0 449154 expected 0 442400 3 0 010450 expected 0 000000 S13 0 561985 expected 0 557600 In the following script we create the histogram of the output of the Ishigami function scf histplot 50 ya xtitle Ishigami function X P X The previous script produces the figure 7 6 74 0 12 0 10 7 0 08 7 0 06 7 0 04 0 02 7 Ishigami function 0 00 Figure 7 6 Histogram of the output of the Ishigami function 75 Chapter 8 Thanks Many thanks to Allan Cornet
17. 76 Bibliography 1 T Homma and A Saltelli Importance measures in global sensitivity analysis of non linear models Reliability Engineering and System Safety 52 1 17 1996 2 Julien Jacques Contributions l analyse de sensibilit et l analyse discriminante g n ral is e 2005 3 Didier Pelat Bases et m thodes pour le traitement des donn es Bruits et Signaux Master M2 Recherche Astronomie astrophysique 2006 4 I M Sobol Sensitivity estimates for nonlinear mathematical models Mathematical Modelling and Computational Experiments 1 407 414 2003 TT
18. Variable transformations In this section we present the transformation of uniform random variables into other types of variables The transformations which are available in the randvar class are presented in figure 4 5 We begin the analysis by a presentation of the theory required to perform transformations Then we present some of the many the transformations which are provided by the library We now present some additionnal details for the function randvar_getvalue rv rv2 value2 This method allows to transform a random variable sample from one law to another The statement 22 Histogram of X 0 9 0 8 0 77 0 6 7 0 3 7 0 277 0 1 0 0 T T T T Figure 4 4 The histogram of a Normal random variable with 1000 samples Source Target Source Target Normale LogNormale Normale Normale Uniforme Uniforme Exponentielle Exponentielle LogNormale LogNormale LogUniforme LogUniforme Source Target Source Target Uniforme LogUniforme Uniforme Uniforme Normale Normale Exponentielle Exponentielle LogNormale LogNormale LogUniforme LogUniforme Source Target Exponentielle Exponentielle Figure 4 5 Variable transfo
19. alli 1 2 p the equation 7 80 can be written as filz fo E Y X 7 88 66 We now plug the equation 7 86 into the previous equality and get fai E Y X E Y 7 89 Similarily we can compute the first order decomposition functions f by computing the conditional expectation with respect to X and X Indeed since the conditional expectation with respect to X is 1 0 for all 2 1 2 p the equation 7 81 can be written as fislei 2 fo fixi fj E Y X Ech 7 91 E Y X Xj E Y E Y X E Y X 7 92 Similarily the equation 7 83 leads to Jug Ung Zj de fo file f 25 Pill Jam Fik Ti Tk HKL Tx E Y X X5 Xx 7 94 E Y X Xj Xx E Y E Y X E Y X E Y Xx 7 95 E Y Xy E Y X X4 Y Xy Xa 7 96 TODO 7 9 Decomposition of the variance TODO 7 10 Total sensitivity indices TODO 7 11 Ishigami function In this section we consider the model Y f Xi Xo X3 sin X1 asin X2 bX3 sin X 7 97 where X X X3 are three random variables uniform in 7 7 This implies that the distribution function of the variable X f satisfies the equation 1 7 98 F X for i 1 2 3 We are going to compute the expectation the variance and the sensitivity indices of this function Before this we need auxialiary results which are presented first 67 7 11 1 Elementary integration We
20. my handle children x label text Variable 1 Normale 1 0 0 5 my _handle children y_label text Variable 2 Uniforme 2 0 3 0 The following script allows to plot the histogram of the two variables which are presented in figures 5 9 and 5 10 Plot Var 1 34 Latin Hypercube Sampling Variable 2 Uniforme 2 0 3 0 1 0 0 5 0 0 0 5 1 0 1 5 2 0 2 5 3 0 Variable 1 Normale 1 0 0 5 Figure 5 8 Latin Hypercube Sampling The first variable is Normal the second variable is Uniform 39 my handle scf clf my handle reset histplot 50 sampling 1 my handle children title text Variable 1 Normale 1 0 0 5 Plot Var 2 my_handle scf clf my handle reset histplot 50 sampling 2 my handle children title text Variable 2 Uniforme 2 0 3 0 Variable 1 Normale 1 0 0 5 Figure 5 9 Latin Hypercube Sampling Normal random variable We can use the mean and variance on each random variable and check that the expected result is computed We insist on the fact that the mean and variance functions are not provided by the NISP library these are pre defined functions which are available in the Scilab library That means that any Scilab function can be now used to process the data generated by the
21. of nisp ort cpp Compilation of nisp pc cpp Compilation of nisp_polyrule cpp 12 Compilation Compilation Compilation Compilation Compilation Compilation of of of of of of nisp_qua cpp nisp_random cpp nisp_smo cpp nisp_util cpp nisp_va cpp nisp_smolyak cpp Building shared library be patient Generate a cleaner file Building gateway Generate a gateway file Generate a loader file Generate a Makefile Running the Compilation Compilation Compilation Compilation Compilation Compilation Compilation Compilation Compilation Compilation Compilation Compilation Compilation Compilation Compilation Compilation Compilation Compilation Compilation Compilation Compilation Compilation Compilation Compilation Compilation Compilation Compilation Compilation Compilation Compilation Compilation Compilation Compilation Compilation Compilation Compilation Makelib makefile of of of of of of of of of of of of of of of of of of of of of of of of of of of of of of of of of of of of nisp_gettoken cpp nisp_gwsupport cpp nisp_PolynomialChaos_map cpp nisp_RandomVariable_map cpp nisp_SetRandomVariable_map cpp sci_nisp_startup cpp sci_nisp_shutdown cpp sci_nisp_verboselevelset cpp Sci nisp verboselevelget cpp Sci nisp initseed cpp Sci randvar new cpp Sci randvar destroy cpp Sci randvar size cpp Sci randvar tokens cpp Sci randvar getlog cpp Sci randvar getval
22. the following Scilab session we use the atomsList function to print the list of all ATOMS toolboxes atomsList ANN_Toolbox ANN Toolbox dde_toolbox Dynamic Data Exchange client for Scilab module_lycee Scilab pour les lyc l es NISP Non Intrusive Spectral Projection plotlib Matlab like Plotting library for Scilab scipad Scipad 7 20 sndfile_toolbox Read amp write sound files stixbox Statistics toolbox for Scilab 5 2 In the following Scilab session we use the atomsShow function to print the details about the NISP toolbox atomsShow NISP Package NISP Title NISP Summary Non Intrusive Spectral Projection Version 2 1 Depend Category ies Optimization Maintainer s Pierre Marechal lt pierre marechal scilab org gt Michael Baudin lt michael baudin scilab org gt 10 Entity CEA DIGITEO WebSite License LGPL Scilab Version gt 5 2 0 Status Not installed Description This toolbox allows to approximate a given model which is associated with input random variables This toolbox has been created in the context of the OPUS project http opus project fr within the workpackage 2 1 1 Construction de m l ta mod les This project has received funding by Agence Nationale de la recherche http www agence nationale recherche fr See in the help provided in the help en_US directory of the toolbox for more information about its use Use cases are presented in the de
23. the relationship 7 4 The linear correlation coefficient between Y and X is Cow Y X PY X L 7 10 VV VV Xi for i 1 2 p In the particular case of the affine model 7 4 we have Cov Y Xi Cov fo Xi amp Cov Xi Xi B2Cov X2 Xi 7 11 BiCov Xi Xi suu By Cov Xp Xi 7 12 7 13 Since the random variables X are independent we have Cov X X 0 for any j i There fore Cov Y X BiCov Xi Xi 7 14 BV 7 15 Hence the correlation coefficient can be simplified into az iV Xi 7 16 PT WOE in _ Bi V X 7 17 V Y 54 We square the previous equality and get 2 2 BV Xi d Sa Za 7 18 Py X V Y Therefore the square of the linear correlation coefficient is equal to the first order sensitivity index i e Pix SRC 7 19 7 4 Using scatter plots In this section we present an example of an affine model where the difference between local and global sensitivity is made clearer by the use of scatter plots Assume four independent random variables X for i 1 2 3 4 We assume that the variables X are normally distributed with zero mean and i variance Let us consider the affine model Notice that the derivative of Y with respect to any of its input is equal to one i e OY 1 7 21 for 1 2 3 4 This means that locally the inputs all have the same effect on the output Let us compute the standardized regression coefficients of this model B
24. toolbox for ivar 1 2 m mean sampling ivar mprintf Variable 4d LMeany 04f n ivar m v variance sampling ivar mprintf Variable 4d Variance Af n ivar v end The previous script produces the following output Variable 1 Mean 1 000000 Variable 1 Variance 0 249925 Variable 2 Mean 2 500000 Variable 2 Variance 0 083417 Our numerical simulation is now finished but we must destroy the objects so that the memory managed by the toolbox is deleted 36 Variable 2 Uniforme 2 0 3 0 Figure 5 10 Latin Hypercube Sampling Uniform random variable randvar_destroy vul randvar_destroy vu2 setrandvar_destroy srv 5 2 4 Other types of DOEs The following Scilab session allows to generate a Monte Carlo sampling with two uniform variables in the interval 1 1 The figure 5 11 presents this sampling and the figures 5 12 and 5 13 present the histograms of the two uniform random variables vul randvar new Uniforme 1 0 1 0 vu2 randvar new Uniforme 1 0 1 0 Srv setrandvar new Setrandvar addrandvar srv vul Setrandvar addrandvar srv vu2 setrandvar_buildsample srv MonteCarlo 1000 sampling setrandvar getsample srv randvar destroy vu1 randvar destroy vu2 setrandvar_destroy s
25. 2 330948 output 4 349607 There is actually a much faster way of computing the output Indeed using vectorisation we can compute all the outputs in one single call to the Exemple function This is fast inputdata setrandvar getsample srvu outputdata Exemple inputdata polychaos settarget pc outputdata We can compute the polynomial expansion based on numerical integration so that the coeffi cients of the polynomial are determined This is done with the polychaos_computeexp function which stands for compute the expansion polychaos setdegree pc degre polychaos computeexp pc srvx Integration Everything is now ready for the sensitivity analysis Indeed the polychaos_getmean returns the mean while the polychaos getvariance returns the variance average polychaos getmean pc var polychaos getvariance pc mprintf Mean yuuu u fAn average mprintf Varianceyyuu sf n var mprintf Indice de sensibilite du ler ordre n mprintf uuuuVariable X10 14f n polychaos_getindexfirst pc 1 mprintf uuuuVariable X20 14f n polychaos_getindexfirst pc 2 mprintf Indice de sensibilite Totale n 45 mprintf mprintf uVariable X1 4f n polychaos_getindextotal pc 1 uVariable X2 4f n polychaos_getindextotal pc 2 The previous script produces the following output Mean
26. Carlo design with 2 variables In this section we create a Monte Carlo design with 2 variables We are going to use the exponential distribution function which is not defined in Scilab The following exppdf function computes the probability distribution function of the exponential distribution function 31 Empirical histogram of X2 Figure 5 5 Monte Carlo Sampling Histogram and exact distribution functions for the second variable function p exppdf x lambda p lambda exp lambda x endfunction The following script creates a Monte Carlo sampling where the first variable is Normal and the second variable is Exponential Then we compare the empirical histogram and the exact distribution function We use the dnorm function defined in the Stixbox module _initseed 0 rvi randvar_new Normale 1 0 0 5 rv2 randvar new Exponentielle 5 Definition d un groupe de variables aleatoires srv setrandvar new Setrandvar addrandvar srv rvi setrandvar_addrandvar srv rv2 nisp np 1000 setrandvar_buildsample srv MonteCarlo np sampling setrandvar_getsample srv Check sampling of random variable 1
27. It is now necessary to setup your Scilab system so that the toolbox is loaded automatically at startup The way to do this is to configure the Scilab startup configuration file The directory where this file is located is stored in the Scilab variable SCIHOME In the following Scilab session we use Scilab v5 2 0 beta 1 in order to know the value of the SCIHOME global variable gt SCIHOME SCIHOME C Users baudin AppData Roaming Scilab scilab 5 2 0 beta 1 On my Linux system the Scilab 5 1 startup file is located in home myname Scilab scilab 5 1 scilab On my Windows system the Scilab 5 1 startup file is located in C Users myname AppData Roaming Scilab scilab 5 1 scilab This file is a regular Scilab script which is automatically loaded at Scilab s startup If that file does not already exist create it Copy the following lines into the scilab file and configure the path to the toolboxes stored in the SCILABTBX variable exec C tbxnisp loader sce The following script presents the messages which are generated when the unit tests script of the toolbox is launched gt exec C tbxnisp runtests sce Tests beginning the 2009 11 18 at 12 47 45 TMPDIR C Users baudin AppData Local Temp SCI_TMP_6372_ 001 004 tbxnisp nisp passed ref created 002 004 tbxnisp polychaosi passed ref created 003 004 tbxnisp randvarl passed ref created 004 004 tbxnisp setrandvarl
28. Methods polychaos_setsizetarget pc np polychaos_settarget pc output polychaos_setinput pc invalue polychaos_setdimoutput pc ny polychaos_setdegree pc no polychaos_getvariance pc polychaos_getmean pc Destructor polychaos_destroy pc Static methods tokenmatrix polychaos_tokens nb polychaos_size Figure 6 1 Outline of the methods of the polychaos class tivity analysis This script is based on the NISP methodology which has been presented in the Introduction chapter We will use the figure 1 1 as a framework and will follow the steps in order In the following Scilab script we define the function Example which takes a vector of size 2 as input and returns a scalar as output function y Exemple x y 1 x 1 x 2 endfunction We now create a collection of two stochastic normalized random variables Since the ran dom variables are normalized we use the default parameters of the randvar_new function The normalized collection is stored in the variable srvx vx1 randvar_new Normale vx2 randvar_new Uniforme srvx setrandvar_new setrandvar_addrandvar srvx vxl setrandvar_addrandvar srvx vx2 We create a collection of two uncertain parameters We explicitely set the parameters of each random variable that is the first Normal variable is associated with a mean equal to 1 0 and a standard deviation equal to 0 5 while the seco
29. NISP Toolbox Manual Michael Baudin INRIA Jean Marc Martinez CEA Version 0 2 January 2011 Abstract This document is a brief introduction to the NISP module We present the installation process of the module in binary from ATOMS or from the sources We present the configuration functions and the randvar setrandvar and polychaos classes Contents 1 Introduction LI The OPUS DESEE siii EES EP RES e L2 The NISP DDE eccosc Si RAS Me TETTO RITO TETTE lo The NISP MON lt gt A pl RA ROR SO ORDERS dU Ea eS as 2 Installation dl here 2 2 Installing the toolbox from ATOMS 2 3 Installing the toolbox from the sources 3 Configuration functions 4 The randvar class dl The distribution MEN r A de bars dos ans deb RE de bane ds H ane a All TUBOS Run v Rh Eodcm GU GE Oe Aon de oH doni de wd EG 4 1 2 Parameters of the Log normal distribution 4 1 8 Uniform random number generation da a ce au xe x Rode EERE REPS e d e e e Ee ael ANSE eL bo Pw bed hed ei hee ee URDDO UD DD 42 2 The Oriented Object system A s La lus du hoe RO R4 e REDE nat diem die ERLE e RG DA ABT AREA lt lt LD oce dU KEES HED aus RD ae RS 4 3 2 Variable transformations ER or eee ee ee GR 5 The setrandvar class SB o xod ee a Re RS eR SG ZR CL e Be A uuu ando quao Rae ey ey sa ri A 5 2 1 A Monte Carlo design with 2 v
30. The detailed calling sequence of each function is provided by the online help and will not be reproduced in this document The inline help is presented in the figure 1 2 For example in order to access to the help associated with the randvar class we type the following statements in the Scilab console help randvar The previous statements opens the Help Browser and displays the helps page presented in figure Several demonstration scripts are provided with the toolbox and are presented in the figure 1 4 These demonstrations are available under the question mark in the menu of the Scilab console Finally the unit tests provided with the toolbox cover all the features of the toolbox When we want to know how to use a particular feature and do not find the information we can search in the unit tests which often provide the answer Help Browser NISP Toolbox Table of Contents nisp Functions to configure the NISP library overview An overview of the NISP toolbox polychaos A class to manage a Polynomial Chaos expansion randvar A class to manage a random variable setrandvar A class to manage a set of random variables Figure 1 2 The NISP inline help Help Browser Nom randvar class to manage a random variable SYNOPSIS tokens randvar tokens randvar size randvar new name a b randvar new name a randvar new name randvar destroy rv value ran
31. andardized regression coefficients of affine models In this section we present the standardized regression coefficients of an affine model Assume that the random variables X are independent with mean E X and finite variances V X for i 1 2 p Let us consider the random variable Y as an affine function of the variables X Y bo p BiXi 7 4 i 1 2 p where 68 are real parameters for i 1 2 p The expectation of the random variable Y is E Y B0 BE X 7 5 i 1 2 p Since the variables X are independent the variance of the sum of variables is the sum of the variances Hence V Y V B gt gt V GXI 7 6 i 1 2 p 53 which leads to the equality V Y gt EV 7 7 1 1 2 Hence each term 82V X is the part of the total variance V Y which is caused by the variable Xi The standardized regression coefficient is defined as BV X SRC e 7 8 VY 7 8 for i 1 2 p Hence the sum of the standardized regression coefficients is one SRC SRC2 SRC 1 7 9 7 3 Link with the linear correlation coefficients In this section we present the link between the linear correlation coefficients of an affine model and the standardized regression coefficients Assume that the random variables X are independent with mean E X and finite variances V X for i 1 2 p Let us consider the random variable Y which depends linearily on the variables X by
32. ariables 5 2 2 A Monte Carlo design with 2 variables dl A LHS ISA ooe es e RR A A RR no Othertypes of DOES o eg de EE a E 329 79 4 x 6 The polychaos class 01 Tobin 4 ce een D a hes de Od 0 Dada Dom de YI Sede g Gaede Roe bi oi Ze rg Zeie eee ed Peed VR eee Pa bed Peed da 6 2 1 Product of two random variables Dad The lenge test Case cia a c SED NEN H Co Co D 7 A tutorial introduction to sensitivity analysis 7 1 Sensitivity analysis 7 2 Standardized regression coefficients of affine models 7 3 Link with the linear correlation coefficients 7 4 Using scatter plots 7 5 Sensitivity analysis for nonlinear models 7 6 The effect of the interactions 7 7 Sobol decomposition 7 8 Decomposition of the expectation 7 9 Decomposition of the variance 7 10 Total sensitivity indices 7 11 Ishigami function Elementary integration Expectation s sc s nm ee ee ee a MT 3 og dd Ped Sobol decomposition The Sobol method for sensitivity analysis TILI 711 2 TALA 7 11 4 7 11 5 7 11 6 8 Thanks Bibliography The Ishigami function by the Sobol method 52 92 93 94 99 57 60 65 66 67 67 67 68 69 71 72 72 72 76 77 Chapter 1 Introduction 1 1 The OPUS project The goal of this toolbox is to provide a tool to manage uncertainties in simulated models This toolbox is based on the NISP library where NISP stand
33. ase In this section we present the Ishigami test case The function Exemple is the model that we consider in this numerical experiment This function takes a vector of size 3 in input and returns a scalar output 46 Product function Empirical Histogram 0 45 0 40 0 35 0 30 0 25 7 0 20 0 15 0 10 7 0 05 0 00 y T g T z T 3 T T y T T T li TE 1 Figure 6 3 Product function Histogram of the output on a LHS design with 10000 experiments function y Exemple x a 7 b 0 1 si sin x 1 s2 sin x 2 y 1 st a s2 s2 b x 3 x 3 x 3 x 3 s1 endfunction We create 3 uncertain parameters which are uniform in the interval r m and put these random variables into the collection srvu rvui randvar new Uniforme pi Api rvu2 randvar_new Uniforme pi pi rvu3 randvar_new Uniforme pi pi srvu setrandvar_new setrandvar_addrandvar srvu rvul setrandvar_addrandvar srvu rvu2 setrandvar_addrandvar srvu rvu3 The collection of stochastic variables is created with the function setrandvar_new The calling sequence srvx setrandvar_new nx allows to create a collection of nx 3 random variables uniform in the interval 0
34. cript produces the following output SRC_1 0 03538 SRC_2 0 12570 SRC_3 0 31817 SRC_4 0 54314 expected 0 expected 0 expected 0 expected 0 03333 13333 30000 53333 SUM 1 00066 expected 1 00000 The previous script also produces the plots 7 2 7 3 7 4 and 7 5 7 5 Sensitivity analysis for nonlinear models Let us focus on one particular input X of the model f with 1 2 particular value say x for example then the variance of the output Y surely decreases because the variable X is not randome anymore We can then measure the conditionnal variance given Xi denoted by V Y X z Since X is a random variable the conditionnal variance V Y X x is a random variable We are interested in the average value of this variance that is we are interested in E V Y X If X has a large weight in the variance V Y then E V Y X is small The theorem of the total variance states that if V Y is finite then V Y V E Y X E V Y X If X has a large weight in the variance V Y then V E Y X is large 57 p If we set X toa Scatter plot for X1 25 Figure 7 2 Scatter plot for an affine model Variable X4 Scatter plot for X2 25 Figure 7 3 Scatter plot for an affine model Variable X5 58 Scatter plot for X3 15 20 15 10 5 o 5 10 15 20 25 Figure 7 4 Scatter plot for an affine model Variable X3 Scatter
35. dvar getvalue rv value randvar getvalue rv rv2 value2 randvar getlog rv Figure 1 3 The online help of the randvar function GUI Genetic Algorithms Simulated Annealing Graphics Signal Processing CACSD Optimization and Simulation Polynomials Simulation Scicos Metanet Tcl Tk Sound file handling Random Spreadsheet Nisp scilE 2c Usecase 1 Usecase 2 Randvar 1 Randvar 2 Setrandvar 1 Setrandvar Lhs design Setrandvar Monte Carlo samplin Figure 1 4 Demonstrations provided with the NISP toolbox Chapter 2 Installation In this section we present the installation process for the toolbox We present the steps which are required to have a running version of the toolbox and presents the several checks which can be performed before using the toolbox 2 1 Introduction There are two possible ways of installing the NISP toolbox in Scilab e use the ATOMS system and get a binary version of the toolbox e build the toolbox from the sources The next two sections present these two ways of using the toolbox Before getting into the installation process let us present some details of the the internal components of the toolbox The following list is an overview of the content of the directories tbanisp demos demonstration scripts tbanisp doc the documentation tbanisp doc usermanual the TEXsources of this manual tbznisp etc startup and shutdo
36. e is a unique decompo sition Y fo K icht KS fig tum sho pl Woy eoe Tp 7 70 i 1 2 p 1 lt i lt j lt p where fo is a constant and the function of the decomposition satisfy the egualities 1 d gege Tris Les Zi dEi 0 7 71 0 for any k 1 2 8 and any indices 414 C 1 p The equalities 7 71 mean that the integral of each function with respect to one of its variables is zero The immediate consequence of this is that the decomposition functions are orthogonal i e 1 d hd EM DEE i CC E yas j d SSC dip 0 7 72 0 if in i Chyn Js This because if the two set of indices 41 4s and ji js this means that there is at least one index k which appears in one index and not in the other By the equality 7 71 this implies that if we integrate with respect to xy then the integral is zero Since the integral in 7 72 is for all the variables since implies that all the integral is zero We are now going to explicitely compute the decomposition functions fo fi fij etc by integration the decomposition using the orthogonality to simplify the results If we integrate the equation 7 70 with respect to all the variables we get Id e fo 7 73 If we integrate the equation 7 70 with respect to all the variables except 7 we get fx daa fo fA xa 7 74 where z is the vector x without its i th component i e Tri an ER ERN ip 7 75 If we integrate the equation 7 70 with
37. eate two random variables gt vui randvar_new Uniforme vul 3 vu2 randvar_new Uniforme vu2 4 Assume now that we have lost the token associated with the variable vu2 We can easily simulate this situation by using the clear which destroys a variable from Scilab s memory space gt clear vu2 randvar getvalue vu2 l error 4 Undefined variable vu2 It is now impossible to generate values from the variable vu2 Moreover it may be difficult to know exactly what went wrong and what exact variable is lost At any time we can use the randvar tokens function in order to get the list of current variables Deleting these variables allows to clean the memory space properly without memory loss randvar tokens ans 3 4 randvar destroy 3 ans 3 randvar destroy 4 ans 4 randvar tokens ans L 4 3 Examples In this section we present to examples of use of the randvar class The first example presents the simulation of a Normal random variable and the generation of 1000 random variables The second example presents the transformation of a Uniform outcome into a LogUniform outcome 21 4 3 1 A sample session We present a sample Scilab session where the randvar class is used to generate samples from the Normale law In the following Scilab session we create a Normale random variable and compute samples from this law The nisp_initseed function is used to initial
38. expansion are computed based on given numerical experiments which creates the association between the inputs and the outputs Once computed the expansion can be used as a regular function The mean standard deviation or quantile can also be directly retrieved The tool allows to get the following results e mean e variance e quantile e correlation etc Moreover we can generate the C source code which computes the output of the polynomial chaos expansion This C source code is stand alone that is it is independent of both the NISP library and Scilab It can be used as a meta model The figure 6 1 presents the most commonly used methods available in the polychaos class More methods are presented in figure 6 2 The inline help contains the detailed calling sequence for each function and will not be repeated here More than 50 methods are available and most of them will not be presented here More informations about the Oriented Object system used in this toolbox can be found in the section 4 2 2 6 2 Examples In this section we present to examples of use of the polychaos class 6 2 1 Product of two random variables In this section we present the polynomial expansion of the product of two random variables We analyse the Scilab script and present the methods which are available to perform the sensi Al Constructors pc polychaos new file pc polychaos new srv ny pc polychaos_ new pc nopt varopt
39. first notice that the integral of the sin function in the interval 7 m is zero since this function is symetric Hence E sin x dx 0 We are going to prove that D sin x dx T Indeed if we integrate the sin z function by part we get ina Labor Lote J TT TT 0 f cos x dx T On the other hand the equality cos x sin z 1 implies f sin z dz J EN yd T T an cos x dx T We now combine 7 102 and 7 104 and get cos 2 dx 2r gd cos z da The previous equality implies that S cos x dx 2m which leads to f cos x dx 7 Finally the previous equality combined with 7 102 immediately leads to 7 100 We are going to prove that Uu IT H sin a de e Indeed if we integrate the sin x function by part we get 8 sin r dr cos x sin z n cos 2 3 sin z cos x dx T T 0 3 cos a sin z da T 68 7 99 7 100 7 101 7 102 7 103 7 104 7 105 7 106 7 107 7 108 7 109 7 110 On the other hand the equality cos x sin z 1 implies n sin x dx sin x sin 1 dx 7 111 T L 1 cos x sin x dx 7 112 T sin x dx I cos a sin z da 7 113 T T We plug the equality 7 100 into the previous equation and get f i sin z da f cos a sin z da 7 114 TT 71 We combine 7 110 and 7 110 and get 3 H i cos z sin x d
40. he expectation of Y is E Y f X4 XoF Xi Xa dx1dx2 7 37 60 where iz 2 is the joint probability distribution function of the variables X and X Since X and X are independent variables we have F 21 2 FX E Xo 7 38 where Fj is the probability distribution function of X and F gt is the probability distribution function of X5 Then we have E Y i X X3 FACX1 F5 X5 dz1dz2 tk XF An lm D WEEN E X E X3 Therefore The variance of Y is The expectation E Y is EV fS IRE E XIF Xde H i XEM E Now we have which leads to Therefore E Y V X E X1 V X5 E X2 oi ui ea 15 Finally we get V Y ei ui e3 12 uua 61 7 39 7 40 7 41 7 49 7 43 7 44 7 45 7 46 7 47 We can expand the previous equality and get V Y 0703 0143 1702 pis uiuo The last two terms of the previous equality can be simplified so that we get V Y 0103 0113 1503 The sensitivity indices can be computed from the definitions s V YIX a _ V E Y X2 Es ange ANA VI We have E Y X1 E X2 X u2X1 Similarily F Y X2 ka Hence V m2X3 ELSE V X2 Sy VY S We get i pV Xy Le ECT V Y Di H V Y Finally the first order sensitivity indices are 8 H207 V Y g pio V Y Since 0202 gt 0 we have 2 2 2 2 2 2 22 22 1201 mo V Y 0703
41. ize the seed for the uniform random variable generator Then we use the randvar_new function to create a new random variable from the Normale law with mean 1 and standard deviation 0 5 The main loop allows to compute 1000 samples from this law based on calls to the randvar_getvalue function Once the samples are computed we use the Scilab function mean to check that the mean is close to 1 which is the expected value of the Normale law when the number of samples is infinite Finally we use the randvar_destroy function to destroy our random variable Once done we plot the empirical distribution function of this sample with 50 classes nisp_initseed 0 mu 1 0 Sigma 0 5 rv randvar new Normale mu sigma nbshots 1000 values zeros nbshots for i 1 nbshots values i randvar_getvalue rv end mymean mean values mysigma st deviation values myvariance variance values mprintf Meanj is j fi expected 2f in mymean mu mprintf Standard deviationy isy yhfyCexpectedy 4f n mysigma sigma mprintf Variance isy 1 4fu expected 0 4f n myvariance sigma 2 randvar_destroy rv histplot 50 values xtitle Histogram of X X P x The previous script produces the following output Mean is 0 988194 expected 1 000000 Standard deviation is 0 505186 expected 0 500000 Variance is 0 255213 expected 0 250000 The previous script also produces the figure 4 4 4 3 2
42. library Indeed the random number generator is based on the compiler so that its quality cannot be guaranted The Uniforme law is associated with the parameters a b R with a lt b It produces real values uniform in the interval a b To compute the uniform random number X in the interval a b a uniform random number in the interval 0 1 is generated and then scaled with X a b a X 4 6 Let us now analyse how the uniform random number X 0 1 is computed The uniform random generator is based on the C function rand which returns an integer n in the interval 18 0 RAND_MAX The value of the RAND M AX variable is defined in the file stdlib h and is compiler dependent For example with the Visual Studio C 2008 compiler the value is RAND_MAX 2 1 32767 4 7 A uniform value X in the range 0 1 is computed from ze N where N RAND MAX and n 0 RAN D M AX X 4 8 4 2 Methods In this section we give an overview of the methods which are available in the randvar class 4 2 1 Overview The figure 4 3 presents the methods available in the randvar class The inline help contains the detailed calling sequence for each function and will not be repeated here Constructors rv randvar_new type options Methods value randvar_getvalue rv options randvar_getlog rv Destructor randvar_destroy rv Static methods rvlist randvar_tokens nbrv
43. mean sampling 1 variance sampling 1 Check sampling of random variable 2 min sampling 2 max sampling 2 Plot scf histplot 40 sampling 1 x linspace 1 3 1000 p dnorm x 1 0 5 plot x p r 32 xtitle Empirical histogram of X1 X P X legend Empirical Exact scf histplot 40 sampling 2 x linspace 0 2 1000 p exppdf x 5 plot x p r xtitle Empirical histogram of X2 X P X legend Empirical Exact Clean up setrandvar_destroy srv randvar destroy rv1 randvar destroy rv2 The previous script produces the figures 5 6 and 5 7 Empirical histogram of X1 Empirical Exact Figure 5 6 Monte Carlo Sampling Histogram and exact distribution functions for the first variable 5 2 3 A LHS design In this section we present the creation of a Latin Hypercube Sampling In our example the DOE is based on two random variables the first being Normal with mean 1 0 and standard deviation 0 5 and the second being Uniform in the interval 2 3 We begin by defining two random variables with the randvar new function vul randvar_new Normale 1 0 0 5 vu2 randvar_new Uniforme 2 0 3 0 Then we create a collection of random variable
44. mi Sensitivity indices Fonction Ishigami Indice de sensibilit CCI totale CU premier ordre Figure 6 5 Ishigami function Sensitivity indices 51 Chapter 7 A tutorial introduction to sensitivity analysis In this chapter we extremely briefly present the theory which is used in the library This section is a tutorial introduction to the NISP module 7 1 Sensitivity analysis In this section we present the sensitivity analysis and emphasize the difference between global and local analysis Consider the model Y f X 7 1 where X Dx C R is the input and Y Dy CR is the output of the model The mapping f is presented in figure 7 1 Figure 7 1 Global analysis The assume that the input X is a random variable so that the output variable Y is also a random variable We are interested in measuring the sensitivity of the output depending on the uncertainty of the input More precisely we are interested in knowing e the input variables X which generate the most variability in the output Y e the input variables X which are not significant e a sub space of the input variables where the variability is maximum 92 e if input variables interacts Consider the mapping presented in figure 7 1 The f mapping transforms the domain Dy into the domain Dy If f is sufficiently smooth small perturbations of X generate small perturbations of Y The local sensitivity analysis focuses o
45. mos directory In the following Scilab session we use the atomsInstall function to download and install the binary version of the toolbox corresponding to the current operating system gt atomsInstall NISP ans INISP 2 1 allusers D Programs SC3623 1 contrib NISP 2 1 I The allusers option of the atomsInstall function can be used to install the toolbox for all the users of this computer We finally load the toolbox with the atomsLoad function gt atomsLoad NISP Start NISP Toolbox Load gateways Load help Load demos ans INISP 2 1 D Programs SC3623 1 contrib NISP 2 1 Now that the toolbox is loaded it will be automatically loaded at the next Scilab session 2 3 Installing the toolbox from the sources In this section we present the steps which are required in order to install the toolbox from the sources In order to install the toolbox from the sources a compiler is required to be installed on the machine This toolbox can be used with Scilab v5 1 and Scilab v5 2 We suppose that the archive has been unpacked in the tbxnisp directory The following is a short list of the steps which are required to setup the toolbox 1 build the toolbox run the tbrnisp builder sce script to create the binaries of the library create the binaries for the gateway generate the documentation 2 load the toolbox run the tbrnisp load sce script to load all commands and setup the documentation 11 3 setup the sta
46. n parameters srvu setrandvar_new setrandvar_addrandvar srvu vul setrandvar_addrandvar srvu vu2 5 Create the Design Of Experiments degre 2 setrandvar buildsample srvx Quadrature degre Setrandvar buildsample srvu srvx 6 Create the polynomial ny 1 pc polychaos new srvx ny np setrandvar getsize srvx mprintf Number of experiments 4d n np polychaos setsizetarget pc np T Perform the DOE inputdata setrandvar getsample srvu outputdata Exemple inputdata 63 polychaos settarget pc outputdata 8 Compute the coefficients of the polynomial expansion polychaos setdegree pc degre polychaos computeexp pc srvx Integration 9 Get the sensitivity indices average polychaos getmean pc Ey expectation mui mu2 var polychaos getvariance pc Vy expectation mu2 2 sigmai 2 mui 2 sigma2 2 sigmai 2 sigma2 2 mprintf Meanyuuu f y expectation y f 1n average Ey_expectation mprintf Varianceyyuu sf yu expectation f n var Vy_expectation mprintf First order sensitivity index n Si polychaos getindexfirst pc 1 S1 expectation mu272 sigmai172 Vy_expectation mprintf QuuuVariable X1 Y 4f Cexpectationy 4f n S1 S1_expectation re abs Si S81 expectation S1 expectation mprintf GQuuuuuuuRelativeyErrory 4f n re S2 polychaos_getindexfirst pc 2 S2_ex
47. n the behaviour of the mapping in the neighbourhood of a particular point X Dx toward a particular point Y Dy The global sensisitivity analysis models the propagation of uncertainties transforming the whole set Dx into the set Dy In the following we assume that there is only one output variable so that Y R There are two types of analysis that we are going to perform that is uncertainty analysis and sensitivity analysis In uncertainty analysis we assume that fx is the probability density function of the variable X and we are searching for the probability density function fy of the variable Y and by its cumulated density function Fy y P Y lt y This problem is difficult in the general case and this is why we often are looking for an estimate of the expectation of Y as defined by E Y y fy gydi 7 2 and an estimate of its variance Var Y f y E Y fy y dy 7 3 We might also be interested in the computation of the probability over a threshold In sensitivity analysis we focus on the relative importance of the input variable X on the uncertainty of Y This way we can order the input variables so that we can reduce the variability of the most important input variables in order to finally reduce the variability of Y More details on this topic can be found in the papers of Homma and Saltelli 1 or in the work of Sobol 4 The Thesis by Jacques 2 presents an overview of sensitivity analysis 7 2 St
48. nd Uniform variable is in the interval 1 0 2 5 This collection is stored in the variable srvu vul randvar_new Normale 1 0 0 5 vu2 randvar_new Uniforme 1 0 2 5 srvu setrandvar_new setrandvar_addrandvar srvu vul setrandvar_addrandvar srvu vu2 42 Methods output polychaos_gettarget pc np polychaos_getsizetarget pc polychaos getsample pc k ovar polychaos getquantile pc k polychaos getsample pc polychaos getquantile pc alpha polychaos getoutput pc polychaos getmultind pc polychaos_getlog pc polychaos getinvquantile pc threshold polychaos getindextotal pc polychaos getindexfirst pc ny polychaos_getdimoutput pc nx polychaos_getdiminput pc p polychaos getdimexp pc no polychaos_getdegree pc polychaos getcovariance pc polychaos getcorrelation pc polychaos_getanova pc polychaos generatecode pc filename funname polychaos computeoutput pc polychaos computeexp pc srv method polychaos computeexp pc pc2 invalue varopt polychaos buildsample pc type np order Figure 6 2 More methods from the polychaos class 43 The first design of experiment is build on the stochastic set srvx and based on a Quadrature type of DOE Then this DOE is transformed into a DOE for the uncertain collection of parameters srvu degre 2 setrandvar_buildsample
49. om variables is the sum of the expectations be the variables independent or not we have E Y E sin X E asin X4 E bX2 sin X E sin X aE sin X The expectation of the variable X3 sin X is E X3 sin X1 i T fa r3 sin z1 g1 2 3 21 22 z3 dzydzodzg ET set tr FIS x3 sin 21 g1 21 92 22 g3 x3 dz1dx2dx3 E X3 E sin X 70 bE X2 sin X4 7 136 7 137 7 138 7 139 7 140 Hence the expectation of Y is E Y E sin Xi aE sin X1 bE X3 E sin X 7 141 We now combine the equations 7 122 7 127 and 7 135 and get 1 1 E Y Dias den 0 7 142 a e 14 2 7 143 7 11 3 Variance The variance of the output Y is V Y E Y EY 7 144 E sin X asim X3 bX4 sin Xi EM 7 145 E sin X4 a sin X2 b X3 sin X1 7 146 2 sin X a sin X5 2sin X bX34 7 147 2a sin X2 bX3 sin X1 E Y 7 148 E sin X1 a E sin X b E X E sin X1 7 149 2a E sin X E sin X2 20E sin X4 E X1 7 150 2abE sin X3 E X3 E sin X1 E Y y 7 151 By the equality 7 122 the expectation of sin X is zero in the interval 7 7 Therefore the terms associated with E sin X1 can be simplified in the previous equality This leads to V Y E sin Xi a E sin X2 amp E X E sin X1 7 152 2bE X3 E sin X1 E Y 7 153 We now compute the terms appearing in the previous equalit
50. orm distribution function We store the tokens in the variables vu1 vu2 and vu3 These variables are regular Scilab double precision floating point numbers Each value is a token which represents a random variable stored in the toolbox memory space gt vul randvar_new Uniforme vul 0 vu2 randvar new Uniforme vu2 1 gt vu3 randvar_new Uniforme vu3 2 There are now 3 objects in current use as indicated by the following statements The tokenmatrix is a row matrix containing regular double precision floating point numbers nb randvar size nb 3 gt tokenmatrix randvar_tokens tokenmatrix 0 1 2 We assume that we have now made our job with the random variables so that it is time to destroy the random variables We call the randvar_destroy functions which destroys the variables gt randvar_destroy vul randvar destroy vu2 randvar destroy vu3 We can finally check that there are no random variables left in the memory space 20 gt nb randvar_size nb 0 gt tokenmatrix randvar tokens tokenmatrix L Scilab is a wonderful tool to experiment algorithms and make simulations It happens some times that we are managing many variables at the same time and it may happen that at some point we are lost The static methods provides tools to be able to recover from such a situation without closing our Scilab session In the following session we cr
51. oupempty cpp Compilation of sci polychaos getgroupinter cpp Compilation of sci polychaos getinvquantile cpp Compilation of sci polychaos buildsample cpp Compilation of sci polychaos getoutput cpp Compilation of sci polychaos getquantile cpp Compilation of sci polychaos getquantwilks cpp Compilation of sci polychaos getsample cpp Compilation of sci polychaos setgroupaddvar cpp Compilation of sci polychaos computeoutput cpp Compilation of sci polychaos setinput cpp Compilation of sci polychaos propagateinput cpp Compilation of sci_polychaos_getanova cpp Compilation of sci polychaos setanova cpp Compilation of sci polychaos getanovaord cpp Compilation of sci polychaos getanovaordco cpp Compilation of sci polychaos realisation cpp Compilation of sci polychaos save cpp Compilation of sci polychaos generatecode cpp Building shared library be patient Generate a cleaner file Generating loader gateway sce Building help Building the master document C tbxnisp help en_US Building the manual file javaHelp C tbxnisp help en_US Please wait building Generating loader sce in this can take a while The following script presents the messages which are generated when the loader of the toolbox 14 is launched The loader script performs the following steps e load the gateway and the NISP library e load the help e load the demo exec C tbxnisp loader sce Start NISP Toolbox Load gateways Load help Load demos
52. pe X1 et X2 0 557674 The function polychaos getanova prints the functionnal decomposition of the normalized variance polychaos getanova pc The previous script produces the following output 00 0 313953 0 442325 1 55229e 009 8 08643e 031 0 243721 7 26213e 031 1 6007e 007 k o HO HO KA Fk OO OO F Fa ki L k ki We can compute the density function associated with the output variable of the function In order to compute it we use the polychaos_buildsample function and create a Latin Hypercube Sampling with 10000 experiments The polychaos_getsample function allows to quiery the polynomial and get the outputs We plot it with the histplot Scilab graphic function which produces the figure 6 4 polychaos buildsample pc Lhs 10000 0 sample output polychaos getsample pc scf histplot 50 sample output xtitle Ishigami Histogram We can plot a bar graph of the sensitivity indices as presented in figure 6 5 for i 1 nx indexfirst i polychaos getindexfirst pc i indextotal i polychaos getindextotal pc i end scf 49 Fonction Ishigami Histogramme normalis 0 12 0 10 0 08 0 06 0 04 0 02 ll 10 5 0 0 00 mJ 10 15 20 Figure 6 4 Ishigami function Histogram of the variable on a LHS design with 10000 experiments 90 bar indextotal 0 2 blue bar indexfirst 0 15 yellow legend Total First jorder pos 1 xtitle Ishiga
53. pectation mui 2 sigma2 2 Vy expectation mprintf uuuuVariable X2u 14 fu expectationy 4f n S2 S2_expectation re abs S2 S2_expectation S2_expectation mprintf QuuuuuuuRelativeyErrory f n re mprintf Total sensitivity index n STi polychaos_getindextotal pc 1 mprintf j j ji Variable X1 7 4 f Nn ST1 ST2 polychaos_getindextotal pc 2 mprintf j j j Variable X2 7 4 f Nn 8T2 Clean up polychaos destroy pc randvar destroy vul randvar destroy vu2 randvar destroy vx1 randvar destroy vx2 setrandvar_destroy srvu setrandvar_destroy srvx The previous script produces the following output Mean 5 250000 expectation 5 250000 Variance 18 687500 expectation 18 687500 First order sensitivity index Variable X1 0 163880 expectation 0 163880 Relative Error 0 000000 Variable X2 0 752508 expectation Relative Error 0 000000 Total sensitivity index Variable X1 0 247492 Variable X2 0 836120 0 752508 We see that the polynomial chaos performs an exact computation 64 7 7 Sobol decomposition Sobol 4 introduced the sensitivity index based on V E Y X by decomposing the function f as a sum of function with an increasing number of parameters We consider the function f Y edd dish 7 69 where x x1 p 0 1 P If f can be integrated in 0 1 then ther
54. plot for X4 25 Figure 7 5 Scatter plot for an affine model Variable X4 99 The first order sensivity indice of Y to the variable X is defined by s V E Y X5 vay 7 30 fori 1 2 p The sensitivity indice measures the part of the variance which is caused by the uncertainty in X We can compute the sensitivity indice when the function f is linear Assume that the output Y depends linearily on the input X Y B M GX 7 31 i 1 2 p where 5 ER for 0 1 2 p Then E Y X Bo y E X GiXi 7 32 since the expection of a sum is the sum of expectations Then V E Y A V 8 Xi 7 33 since the variance of a constant term is zero Therefore the sensitivity index of Y to the variable X is BV Xi S Yy 7 35 for i 1 2 p Hence if we make the assumption that a model is affine then the empirical linear correlation coefficient can be used to estimate the sensitivity indices 7 6 The effect of the interactions In this example we consider a non linear non additive model made of the product of two inde pendent random variables The goal of this example is to show that in some cases we have to consider the interations between the variables Consider the function Y X Xo 7 36 where X and X3 are two independent normal random variables with mean 44 and uz and variances 2 2 oi and 03 Let us compute the expectation of the random variable Y T
55. printf First order sensitivity index n mprintf Guu Variable X1 4f n polychaos_getindexfirst pc 1 mprintf Quu Variable X2 4f n polychaos_getindexfirst pc 2 mprintf uuuuVariablewX3u u4fAn polychaos_getindexfirst pc 3 mprintf Total sensitivity index n mprintf uuuuVariable X10 14f n polychaos_getindextotal pc 1 mprintf uuuuVariable X20 14f n polychaos_getindextotal pc 2 mprintf uuuuVariablewX3u u4fAn polychaos_getindextotal pc 3 The previous script produces the following output Mean 3 500000 Variance 13 842473 First order sensitivity index Variable X1 0 313953 Variable X2 0 442325 Variable X3 0 000000 48 Total sensitivity index Variable X1 0 557675 Variable X2 0 442326 Variable X3 0 243721 We now focus on the variance generated by the variables 1 and 3 We set the group to the empty group with the polychaos_setgroupempty function and add variables with the polychaos_setgroupaddvar function groupe 1 3 polychaos_setgroupempty pc polychaos_setgroupaddvar pc groupe 1 polychaos_setgroupaddvar pc groupe 2 mprintf Fraction of the variance of a group of variables n mprintf uuuuGroupeuXiswetuX2u 4fAn polychaos_getgroupind pc The previous script produces the following output Fraction of the variance of a group of variables Grou
56. respect to all the variables except and 7 we get 1 f die fo fixi Tm Ta i ae 7 76 0 65 If we integrate the equation 7 70 with respect to all the variables except and j and k we get 1 Pr ld j k fo fi zi AE fx Gk 7 77 Tog Fik Ti Y es Tx dl Tx 7 78 The previous computations allows to get the decomposition functions f x dx 7 79 1 fe fo f f w dz 7 80 0 fine fo file w iss 7 81 agi n NE Ek fijlTi Pal Lp Tx 7 82 1 hie drv 7 83 0 until the last term Jig phone m films 7 84 1245 K i peng ip 1 EH E TT 7 85 1 lt i1 lt lt ip_1 lt p The last term is obviously so that the equality 7 70 is satisfied We have considered a function where the variables are in 0 1 P In fact when we consider the more general model Y f X X where the random variables X are independent and uniform in 0 1 the decomposition 7 70 is still valid 7 8 Decomposition of the expectation We can consider the decomposition 7 70 in terms of expectation and variance If we compute the expectation of Y by the expression 7 70 we get E Y fo 7 86 which is an obvious consequence of the zero integral property 7 71 We can compute the first order decomposition functions f by computing the conditional expectation with respect to X Indeed since the conditional expectation with respect to X is E Y X f ro ds 7 87 for
57. rmations available in the randvar class 23 value randvar getvalue rv rv2 value2 returns a random value from the distribution function of the random variable rv by transformation of value2 from the distribution function of random variable rv2 In the following session we transform a uniform random variable sample into a LogUniform variable sample We begin to create a random variable rv from a LogUniform law and parameters a 10 b 20 Then we create a second random variable rv2 from a Uniforme law and parameters a 2 b 3 The main loop is based on the transformation of a sample computed from rv2 into a sample from rv The mean allows to check that the transformed samples have an mean value which corresponds to the random variable rv nisp initseed O a 10 0 b 20 0 rv randvar_new LogUniforme a b rv2 randvar new Uniforme 2 3 nbshots 1000 valuesLou zeros nbshots for i 1 nbshots valuesUni i randvar_getvalue rv2 valuesLou i randvar getvalue rv rv2 valuesUni i end computed mean valuesLou mu b a log b log a expected mu mprintf Expectation 4 5f expected 4 5f n computed expected scf histplot 50 valuesUni xtitle Empirical histogram Uniform variable X P X scf histplot 50 valuesLou xtitle Empirical histogram Log Uniform variable X P X randvar_destroy rv randvar_destroy rv2 The pre
58. rtup configuration file of your Scilab system so that the toolbox is known at startup see below for details 4 run the unit tests run the tbrnisp runtests sce script to perform all unit tests and check that the toolbox is OK 5 run the demos run the tbrnisp rundemos sce script to run all demonstration scripts and get a quick interactive overview of its features The following script presents the messages which are generated when the builder of the toolbox is launched The builder script performs the following steps e compile the NISP C library e compile the C gateway library the glue between the library and Scilab e generate the Java help files from the xml files e generate the loader script exec C tbxnisp builder sce Building sources Generate a loader file Generate a Makefile Running the Makefile Compilation of utils cpp Compilation of blasi_d cpp Compilation of dcdflib cpp Compilation of faure cpp Compilation of halton cpp Compilation of linpack_d cpp Compilation of niederreiter cpp Compilation of reversehalton cpp Compilation of sobol cpp Building shared library be patient Generate a cleaner file Generate a loader file Generate a Makefile Running the Makefile Compilation of nisp gc cpp Compilation of nisp gva cpp Compilation of nisp ind cpp Compilation of nisp index cpp Compilation of nisp inv cpp Compilation of nisp math cpp Compilation of nisp msg cpp Compilation of nisp_conf cpp Compilation
59. rv It is easy to change the type of sampling by modifying the second argument of the setrandvar_buildsa function This way we can create the Petras Quadrature and Sobol sampling presented in figures 5 14 5 15 and 5 16 3T Monte Carlo Sampling Variable 2 Uniforme 1 0 1 0 1 0 0 8 0 6 0 4 0 2 0 0 0 2 0 4 0 6 0 8 1 0 Variable 1 Uniforme 1 0 1 0 Figure 5 11 Monte Carlo Sampling Two uniform variables in the interval 1 1 Variable 1 Uniforme 1 0 1 0 0 84 0 77 0 6 4 0 5 044 0 3 0 2 0 1 0 0 Figure 5 12 Latin Hypercube Sampling First uniform variable in 1 1 38 Variable 2 Uniforme 1 0 1 0 0 8 0 7 7 0 6 4 0 5 1 0 0 8 0 6 0 4 0 2 0 0 0 2 0 4 0 6 0 8 1 0 Figure 5 13 Latin Hypercube Sampling Second uniform variable in 1 1 Petras Sampling LR SE SER EJE SER SEE BE eS ee SSL 0 7 7 IN E do D e o o pe H Ho
60. s for Non Intrusive Spectral Projection This work has been realized in the context of the OPUS project http opus project fr Open Source Platform for Uncertainty treatments in Simulation funded by ANR the french Agence Nationale pour la Recherche http www agence nationale recherche fr The toolbox is released under the Lesser General Public Licence LGPL as all components of the OPUS project 1 2 The NISP library The NISP library is based on a set of 3 C classes so that it provides an object oriented framework for uncertainty analysis The Scilab toolbox provides a pseudo object oriented interface to this library so that the two approaches are consistent The NISP library is release under the LGPL licence The NISP library provides three tools which are detailed below e The randvar class allows to manage random variables specified by their distribution law and their parameters Once a random variable is created one can generate random numbers from the associated law e The setrandvar class allows to manage a collection of random variables This collection is associated with a sampling method such as MonteCarlo Sobol Quadrature etc It is possible to build the sample and to get it back so that the experiments can be performed e The polychaos class allows to manage a polynomial representation of the simulated model One such object must be associated with a set of experiments which have been performed
61. s with the setrandvar_new function which creates here an empty collection of random variables Then we add the two random variables to the collection 33 Empirical histogram of X2 Empirical Exact Figure 5 7 Monte Carlo Sampling Histogram and exact distribution functions for the second variable srv setrandvar new setrandvar addrandvar srv vul setrandvar_addrandvar srv vu2 We can now build the DOE so that it is a LHS sampling with 1000 experiments setrandvar buildsample srv Lhs 1000 At this point the DOE is stored in the memory space of the NISP library but we do not have a direct access to it We now call the setrandvar_getsample function and store that DOE into the sampling matrix sampling setrandvar getsample srv The sampling matrix has 1000 rows corresponding to each experiment and 2 columns cor responding to each input random variable The following script allows to plot the sampling which is is presented in figure 5 8 my handle scf clf my handle reset plot sampling 1 sampling 2 my_handle children children children line_mode off my_handle children children children mark_mode on my_handle children children children mark_size 2 my_handle children title text Latin Hypercube Sampling
62. section present several examples 5 1 Introduction The setrandvar class allows to manage a collection of random variables and to build a Design Of Experiments DOE Several types of DOE are provided e Monte Carlo e Latin Hypercube Sampling e Smolyak Once a DOE is created we can retrieve the information experiment by experiment or the whole matrix of experiments This last feature allows to benefit from the fact that Scilab can natively manage matrices so that we do not have to perform loops to manage the complete DOE Hence good performances can be observed even if the language still is interpreted The figure 5 1 presents the methods available in the setrandvar class A complete description of the input and output arguments of each function is available in the inline help and will not be repeated here More informations about the Oriented Object system used in this toolbox can be found in the section 4 2 2 5 2 Examples In this section we present examples of use of the setrandvar class In the first example we present a Scilab session where we create a Latin Hypercube Sampling In the second part we present various types of DOE which can be generated with this class 5 2 1 A Monte Carlo design with 2 variables In the following example we build a Monte Carlo design of experiments with 2 input random variables The first variable is associated with a Normal distribution function and the second 27 Constructors srv
63. sizetarget function ny 13 pc polychaos new srvx ny np setrandvar getsize srvx polychaos setsizetarget pc np 44 In the next step we perform the simulations prescribed by the DOE We perform this loop in the Scilab language and make a loop over the index k which represents the index of the current experiment while np is the total number of experiments to perform For each loop we get the input from the uncertain collection srvu with the setrandvar_getsample function pass it to the Exemple function get back the output which is then transferred to the polynomial pc by the polychaos_settarget function This is slow for k 1 np inputdata setrandvar_getsample srvu k outputdata Exemple inputdata mprintf Experiment d inputy 4f 4f uoutputu u f n k inputdata 1 inputdata 2 outputdata polychaos settarget pc k outputdata end The previous script produces the following output Experiment 1 input 0 133975 1 169052 output 0 156623 Experiment 2 input 0 133975 1 750000 output 0 234456 Experiment 3 input 0 133975 2 330948 output 0 312288 Experiment 4 input 1 000000 1 169052 output 1 169052 Experiment 5 input 1 000000 1 750000 output 1 750000 Experiment 6 input 1 000000 2 330948 output 2 330948 Experiment 7 input 1 866025 1 169052 output 2 181482 Experiment 8 input 1 866025 1 750000 output 3 265544 Experiment 9 input 1 866025
64. u setrandvar_new setrandvar_addrandvar srvu rvul setrandvar_addrandvar srvu rvu2 setrandvar_addrandvar srvu rvu3 setrandvar_addrandvar srvu rvu4 Create a sampling by a Latin Hypercube Sampling with size 5000 nbshots 5000 setrandvar_buildsample srvu Lhs nbshots sampling setrandvar_getsample srvu Perform the experiments y Exemple sampling Scatter plots y depending on X_i for k 1 4 scf plot y sampling k rx xistr X string k xtitle Scatteryplotyfory xistr xistr Y end Compute the sample linear correlation coefficients rhol lincorrcoef sampling 1 y D SRC1 rho172 SRCiexpected 1 30 mprintf SRC_1 5f expected 5f n SRC1 SRClexpected 77 rho2 lincorrcoef sampling 2 y 96 SRC2 rho272 SRC2expected 4 30 mprintf SRC_2 5f expected 5f n SRC2 SRC2expected 1 rho3 lincorrcoef sampling 3 y SRC3 rho372 SRC3expected 9 30 mprintf SRC_3 5f expected 5f n SRC3 SRC3expected rho4 lincorrcoef sampling 4 y SRC4 rho472 SRC4expected 16 30 mprintf SRC_4 5f expected 5f n SRC4 SRC4expected 77 SUM SRC1 SRC2 SRC3 SRC4 SUMexpected 1 mprintf SUM 5f expected 4 5f n SUM SUMexpected 77 Clean up randvar_destroy rvul randvar_destroy rvu2 randvar_destroy rvu3 randvar destroy rvu4 setrandvar_destroy srvu The previous s
65. ue cpp Sci setrandvar new cpp Sci setrandvar tokens cpp Sci setrandvar size cpp Sci setrandvar destroy cpp Sci setrandvar freememory cpp Sci setrandvar addrandvar cpp Sci setrandvar getlog cpp sci_setrandvar_getdimension cpp sci_setrandvar_getsize cpp sci_setrandvar_getsample cpp sci_setrandvar_setsample cpp Sci setrandvar save cpp Sci setrandvar buildsample cpp Sci polychaos new cpp Sci polychaos destroy cpp Sci polychaos tokens cpp Sci polychaos size cpp Sci polychaos setdegree cpp Sci polychaos getdegree cpp Sci polychaos freememory cpp 13 Compilation of sci polychaos getdimoutput cpp Compilation of sci polychaos setdimoutput cpp Compilation of sci polychaos getsizetarget cpp Compilation of sci polychaos setsizetarget cpp Compilation of sci polychaos freememtarget cpp Compilation of sci polychaos settarget cpp Compilation of sci polychaos gettarget cpp Compilation of sci polychaos getdiminput cpp Compilation of sci polychaos getdimexp cpp Compilation of sci polychaos getlog cpp Compilation of sci polychaos computeexp cpp Compilation of sci polychaos getmean cpp Compilation of sci polychaos getvariance cpp Compilation of sci polychaos getcovariance cpp Compilation of sci polychaos getcorrelation cpp Compilation of sci polychaos getindexfirst cpp Compilation of sci polychaos getindextotal cpp Compilation of sci polychaos getmultind cpp Compilation of sci polychaos getgroupind cpp Compilation of sci polychaos setgr
66. vious script produces the following output Expectation 14 63075 expected 14 42695 The previous script also produces the figures 4 6 and 4 7 The transformation depends on the mother random variable rv1 and on the daughter ran dom variable rv Specific transformations are provided for all many combinations of the two distribution functions These transformations will be analysed in the next sections 24 Empirical histogram Uniform variable x E amp 08 0 6 7 0 47 0 0 2 0 2 1 2 2 2 3 2 4 2 5 2 6 2 7 2 8 2 9 3 0 Figure 4 6 The histogram of a Uniform random variable with 1000 samples 25 Empirical histogram Log Uniform variable 0 18 0 12 7 Figure 4 7 The histogram of a Log Uniform random variable with 1000 samples 26 Chapter 5 The setrandvar class In this chapter we presen the setrandvar class The first section gives a brief outline of the features of this class and the second
67. w scripts for the toolbox tbrnisp help inline help pages tbanisp macros Scilab macros files sci tbrnisp sci gateway the sources of the gateway tbrnisp src the sources of the NISP library tbinisp tests tests tbxnisp tests nonreg tests tests after some bug has been identified tbxnisp tests unit tests unit tests The current version is based on the NISP Library v2 1 9 2 2 Installing the toolbox from ATOMS The ATOMS component is the Scilab tool which allows to search download install and load toolboxes ATOMS comes with Scilab v5 2 The Scilab NISP toolbox has been packaged and is provided mainly by the ATOMS component The toolbox is provided in binary form depending on the user s operating system The Scilab NISP toolbox is available for the following platforms e Windows 32 bits e Linux 32 bits 64 bits e Mac OS X The ATOMS component allows to use a toolbox based on compiled source code without having a compiler installed in the system Installing the Scilab NISP toolbox from ATOMS requires the following steps e atomsList prints the list of current toolboxes e atomsShow prints informations about a toolbox e atomsInstall installs a toolbox on the system e atomsLoad loads a toolbox Once installed and loaded the toolbox will be available on the system from session to session so that there is no need to load the toolbox again it will be available right from the start of the session In
68. x with the setrandvar_getsample function The sampling matrix has np rows and 2 columns one for each input variable nisp_initseed 0 rvui randvar new Normale 1 3 rvu2 randvar new Uniforme 2 3 srvu setrandvar_new setrandvar_addrandvar srvu rvul setrandvar_addrandvar srvu rvu2 i np 5000 setrandvar_buildsample srvu MonteCarlo np sampling setrandvar getsample srvu Check sampling of random variable 1 mean sampling 1 Expectation 1 Check sampling of random variable 2 mean sampling 2 Expectation 2 5 77 scf histplot 50 sampling 1 xtitle Empirical jhistogram of X1 scf histplot 50 sampling 2 xtitle Empirical histogram of X2 Clean up setrandvar_destroy srvu randvar destroy rvu1 randvar destroy rvu2 The previous script produces the following output gt mean sampling 1 Expectation 1 ans 1 0064346 mean sampling 2 Expectation 2 5 ans 2 5030984 The prevous script also produces the figures 5 2 and 5 3 We may now want to add the exact distribution to these histograms and compare The Normal distribution function is not provided by Scilab but is provided by the Stixbox module Indeed the dnorm function of the Stixbox module computes the Normal probability distribution function 20 Empirical histogram of X1
69. y Actually we do not have much to compute since the equalities 7 127 and 7 135 are already available Indeed the equality 7 127 immediately leads to E sn X E sin X2 7 154 o E 7 155 ND What remains to compute is E sin X3 and E X3 By definition the expectation of the random variable X4 is E X3 zd mn 7 156 a lod 7 157 B EEN 7 158 _ Gm 9 7 159 n 8 7 160 71 On the other hand the expectation of the random variable sin X3 is 2 1 4 E sin X3 sin z2 dx 7 161 T Jo The equality 7 108 then implies that IT E sin X2 a dea 7 162 3 a 7 163 We now plug the equalities 7 155 7 135 7 160 and 7 163 into 7 153 and get VY 2 Pon Dent S a 7 164 1 3a br bri a 5 is 5 SE 7 165 l a T T p T z 7 166 7 11 4 Sobol decomposition In this section we perform the Sobol decompotion of the function f as presented in the section rA A By the equation 7 80 we have AGE he OLOR 7 167 TODO 7 11 5 The Sobol method for sensitivity analysis TODO 7 11 6 The Ishigami function by the Sobol method In this section we compute the sensitivity indices of the Ishigami function by the Sobol method We consider three random variables uniform in 7 7 We use Monte Carlo experiments to compute the sensitivity indices The following script allows to perform the analysis function y ishigami x a 7 b 0 1 si sin x 1
70. y hypothesis the variance of each variable is V Xi l V X2 4 7 22 V X3 9 VCX4 16 7 23 Since the variables are independent the variance of the output Y is V Y V X V X3 V X3 V X4 30 7 24 The standardized regression coefficient is BV X SRC 1 29 V Y 30 KS for 1 2 3 4 More specifically we have 1 4 e 7 26 SRC on SRO 7 26 9 16 SRC3 SRC4 Got 3 30 4 30 We have the following inequalities SRC gt S RCs gt SRC gt S RC 7 28 This means that the variable which causes the most variance in the output is X4 while the variable which causes the least variance in the output is X4 The script below performs the analysis with the NISP module The sampling is based on a Latin Hypercube Sampling design with 5000 points 99 function y Exemple x y x 1 20066 x 3 x 4 endfunction function r lincorrcoef x y Returns the linear correlation coefficient of x and y The variables are expected to be column matrices with the same size x x y y mx mean x my mean y sx sqrt sum x mx 2 sy sqrt sum y my 72 r x mx y my sx sy endfunction Initialisation de la graine aleatoire nisp_initseed 0 Create the random variables rvui randvar new Normale 0 1 rvu2 randvar new Normale 0 2 rvu3 randvar_new Normale 0 3 rvu4 randvar_new Normale 0 4 srv

NISP Toolbox Manual - Forge

Contents

Download Pdf Manuals

Related Search

Related Contents