Compute sensitivity indices by brute force.
s = nisp_bruteforcesa ( func , nx ) s = nisp_bruteforcesa ( func , nx , randgen ) s = nisp_bruteforcesa ( func , nx , randgen , n ) [ s , nbevalf ] = nisp_bruteforcesa ( ... )
a function or a list, the name of the function to be evaluated.
a 1-by-1 matrix of floating point integers, the number of inputs of the function.
a function or a list, the random number generator. (default = uniform random variables)
a 1-by-1 matrix of floating point integers (default n=10000), the number of Monte-Carlo experiments, for each sensitivity index
a nx-by-1 matrix of doubles, the first order sensitivity indices
a nx-by-1 matrix of doubles, the actual number of function evaluations.
The algorithm uses the Sobol method to compute the first order sensitivity indices.
This function does not make any hypothesis about the dependency structure of the input random variables, i.e. it can also work if the input random variables are dependent.
Any optional input argument equal to the empty matrix will be set to its default value.
The function should have header y = func ( x ) where x is a m-by-nx matrix of doubles, where m is the number of experiments to perform, nx is the number of input random variables, and y is a m-by-1 matrix of doubles.
It might happen that the function requires additionnal arguments to be evaluated. In this case, we can use the following feature. The function func can also be a list, with header y = f ( x , a1 , a2 , ... ). In this case the func variable should hold the list (f,a1,a2,...) and the input arguments a1, a2, will be automatically be appended at the end of the calling sequence of rg.
The random number generator must have header x = randgen ( nx , n , step , i ) where
There are two phases in the algorithm. During step=1, a basic Monte-Carlo sampling must be generated. During step=2, the values associated with variable #i must be replicated, so that there are only sqrt(n) different values in x(:,i). This allows to compute E(Y|Xi), with sqrt(n) different values of Xi.
In the case where all the input variables are independent, it is straightforward to generate independent samples. In the case where the input variables are dependent, the creation of the sampling must be done according to the dependency structure.
It might happen that the random number generator requires additionnal arguments to be evaluated. In this case, we can use the following feature. The function randgen can also be a list, with header x = rg ( n , step , i , a1 , a2 , ... ). In this case the randgen variable should hold the list (rg,a1,a2,...) and the input arguments a1, a2, will be automatically be appended at the end of the calling sequence of rg.
TODO : test it on a 4 input variables test case.
TODO : compute the total sensitivity indices.
// Compute the first order sensitivity indices of the ishigami function. // Three random variables uniform in [-pi,pi]. function y=ishigami(x) a=7. b=0.1 s1=sin(x(:,1)) s2=sin(x(:,2)) x34 = x(:,3).^4 y(:,1) = s1 + a.*s2.^2 + b.*x34.*s1 endfunction function x=myrandgen(nx, n, step, i) if ( step == 1 ) then x = grand(n,nx,"unf",-%pi,%pi) else sn = sqrt(n) xi = grand(sn,1,"unf",-%pi,%pi) xr = grand(n,nx-1,"unf",-%pi,%pi) x = [xr(:,1:i-1),xi.*.ones(sn,1),xr(:,i:nx-1)] end endfunction a=7.; b=0.1; exact = nisp_ishigamisa(a,b); n = 500; nx = 3; // See the sampling in step 1 step = 1; x = myrandgen ( nx, 4, step, []) // See the sampling in step 2. // See how the variable #3 is replicated. step = 2; x = myrandgen ( nx, 16, step, 3) [ s , nbevalf ] = nisp_bruteforcesa ( ishigami , nx , myrandgen , n ); mprintf("S1 : %f (expected = %f)\n", s(1), exact.S1); mprintf("S2 : %f (expected = %f)\n", s(2), exact.S2); mprintf("S3 : %f (expected = %f)\n", s(3), exact.S3); | ![]() | ![]() |