What is Sensitivity analysis?

Sensitivity analysis is a method for quantitatively evaluating which input uncertainty is derived from the uncertainty of the output of a mathematical model or system. https://en.wikipedia.org/wiki/Sensitivity_analysis

When the output is $ Y $ and the input is $ \ bf X $, the mathematical model is expressed as follows. Uncertainty in an input $ X_i $ also causes uncertainty in $ Y $, but quantitatively how much the uncertainty in $ Y $ depends on which $ X_i $ evaluate.

Y = f({\bf X})

It is used for the following purposes.

--Test how robust the results obtained are when there is uncertainty in the model (for example, assuming uncertain quantities that cannot be observed directly). --Understand the relationship between model inputs and outputs --Improve the certainty of the calculation result by identifying the input with high sensitivity and focusing on the input to increase the certainty. --Model simplification by fixing insensitive inputs --Used to determine which inputs to focus on when calibrating the parameters of a model with many parameters

Sensitivity analysis method

There are several methods for sensitivity analysis. Use properly according to the calculation cost, purpose, and characteristics of the mathematical model.

One-at-a-time (OAT) is the simplest method. As the name suggests, it is a method to check how much the result changes by moving only one parameter by methods such as partial differential coefficient and linear regression. It is very simple and easy to understand, and the calculation cost is low, but since the parameters other than the moving parameters are fixed to the baseline values, the input space cannot be fully searched.

Also, if you try to evaluate with the partial derivative, you can only evaluate the value at that particular baseline locally. Unable to evaluate interactions between input variables or when $ y $ changes non-monotonically with respect to $ x_i $

global sensitivity analysis (Sobol method)

There is a variation-based sensitivity analysis as a method for performing a global sensitivity analysis that solves these problems. It is often called the Sobol method. https://en.wikipedia.org/wiki/Variance-based_sensitivity_analysis

In this method, inputs and outputs are considered as random variables. Break down the output variance into contributions from one input (or multiple inputs). For example, it can be said that 70% of the output variance comes from the fluctuation of $ X_1 $, 20% comes from $ X_2 $, and the remaining 10% comes from the interaction between $ X_1 $ and $ X_2 $.

This method can handle (i) global and (ii) non-linear reactions, and (iii) can measure the interaction of input parameters (effects that appear when multiple parameters are changed at the same time). There is. On the other hand, since it is processed statistically, it is necessary to sufficiently sample the input space, and in general, a large number of calculations with various inputs are required, and the calculation cost is heavy.

What kind of amount should be calculated specifically? Now let the model be a $ d $ dimensional variable $ \ bf X = \ {X_1, X_2, \ dots, X_d } $. It is assumed that all variables are normalized as $ X_i \ in [0,1] $ and are uniformly distributed. The output $ Y $ can generally be written as:

Y=f_{0}+\sum _{i=1}^{d}f_{i}(X_{i})+\sum _{i<j}^{d}f_{ij}(X_{i},X_{j})+\cdots +f_{1,2,\dots ,d}(X_{1},X_{2},\dots ,X_{d})

Where $ f_i $ is a function of $ X_i $ and $ f_ {ij} $ is a function of $ X_i $, $ X_j $.

Calculating the variance of $ Y $

\operatorname {Var}(Y)=\sum _{{i=1}}^{d}V_{i}+\sum _{{i<j}}^{{d}}V_{{ij}}+\cdots +V_{{12\dots d}}

Can be written as the sum of the terms whose variance depends on each $ i $. Here, $ V_i $ is defined by the following formula.

V_{{i}}=\operatorname {Var}_{{X_{i}}}\left(E_{{{\textbf  {X}}_{{\sim i}}}}(Y\mid X_{{i}})\right),

This formula first fixes $ X_i $ to a specific value and integrates it against variables other than $ i $ to calculate the expected value of $ Y $ for a specific $ X_i $, and that amount of $ X_i $. It means to calculate the variance when the is changed.

The quadratic term $ V_ {ij} $ is also defined as follows.

V_{ij} = \operatorname{Var}_{X_{ij}} \left( E_{\textbf{X}_{\sim ij}} \left( Y \mid X_i, X_j\right)\right) - \operatorname{V}_{i} - \operatorname{V}_{j}

There are mainly two types of indexes for evaluating the sensitivity of the parameter $ i $: first-order index $ S_i $ and total-effect index $ S_ {Ti} $. Each is defined as follows.

S_{i}={\frac  {V_{i}}{\operatorname {Var}(Y)}}

S_{{Ti}}={\frac  {E_{{{\textbf  {X}}_{{\sim i}}}}\left(\operatorname {Var}_{{X_{i}}}(Y\mid {\mathbf  {X}}_{{\sim i}})\right)}{\operatorname {Var}(Y)}}=1-{\frac  {\operatorname {Var}_{{{\textbf  {X}}_{{\sim i}}}}\left(E_{{X_{i}}}(Y\mid {\mathbf  {X}}_{{\sim i}})\right)}{\operatorname {Var}(Y)}}

$ S_i $ is a method to evaluate the variance of $ Y $ when $ i $ "only" is changed. When $ i $ is changed alone, it fluctuates by about $ S_i $ on average. On the other hand, the Total-effect index is the contribution of $ i $ when all parameters are changed, including the contribution of interaction with other parameters. If $ S_ {Ti} $ is small, its value will not affect the result even if it is fixed, so it can be used for model simplification.

Monte Carlo sampling of the parameter space is performed to calculate these indicators. Random sampling may be used, but the Quasi-Monte Carlo method is often used as a device for more efficient calculation. In fact, the libraries described below do so. Usually, sampling of several thousand x (number of parameters) is required.

How to use SALib

SALib is one of the libraries that perform global sensitivity analysis. Here, how to use it will be explained. Several methods are implemented in this library, including the Sobol method, but only the Sobol method will be explained here.

Reference: https://www.youtube.com/watch?v=gkR_lz5OptU&feature=youtu.be

The installation method is pip install SALib Requires numpy, scipy, matplotlib as prerequisites.

SALib is performed in the following 4 steps.

Definition of model input
Sample generation
Execution of the mathematical model for which sensitivity analysis is desired
Calculate the sensitivity analysis indicators ($ S_i $ and $ S_ {Ti} $)

The simplest sample code is as follows. Here, a test function called ʻIshigami` is used as a mathematical model.

from SALib.sample import saltelli
from SALib.analyze import sobol
from SALib.test_functions import Ishigami
import numpy as np

# Define the model inputs
problem = {
    'num_vars': 3,
    'names': ['x1', 'x2', 'x3'],
    'bounds': [[-3.14159265359, 3.14159265359],
               [-3.14159265359, 3.14159265359],
               [-3.14159265359, 3.14159265359]]
}

# Generate samples
param_values = saltelli.sample(problem, 1000)

# Run model (example)
Y = Ishigami.evaluate(param_values)

# Perform analysis
Si = sobol.analyze(problem, Y, print_to_console=True)

The output looks like this:

Parameter S1 S1_conf ST ST_conf
x1 0.307975 0.068782 0.560137 0.076688
x2 0.447767 0.057872 0.438722 0.042571
x3 -0.004255 0.052080 0.242845 0.025605

Parameter_1 Parameter_2 S2 S2_conf
x1 x2 0.012205 0.082296
x1 x3 0.251526 0.094662
x2 x3 -0.009954 0.070276

The output looks like this: The plot of the first-order index, second-order index, and total-order index is as follows for easy understanding.

$ S_1 $ and $ S_2 $ have large values, but $ S_3 $ is almost zero. This is because the ʻIshigami` function used here is defined as follows, and $ S_3 $ of $ x_3 $ becomes 0, but there is an interaction with $ x_1 $, so $ S_ {13} $ Will be a positive value.

f(x) = \sin(x_1) + a \sin^2(x_2) + b x_3^4 \sin(x_1)

When you execute the saltelli.sample method, $ N (2D + 2) $ sample points are generated, and you need to find $ Y $ for each. If $ N $ is too small, the error bars will be large, so you need to increase $ N $ until the error bars are small enough. Generally around $ N = 1000 $.

Click here for jupyter notebook including full code

https://gist.github.com/yohm/cbca9052200e377bfdc9dc2bc83b848d

Batch execution using OACIS

I ran a simple function in python here, but in general you'll probably want to use a heavier mathematical model. For example, if one calculation takes several minutes, it becomes unmanageable on one PC, and you want to execute it on a cluster machine. Here, I will introduce the procedure for batch execution using oacis. Since oacis has a Python API, you can use it to submit jobs and refer to the results.

Here, as an example, the script that calculates the Ishigami function used earlier is registered as a simulator. After that, if you use your own mathematical model, please read it as appropriate.

First, I prepared a script to calculate the Ishigami function as follows. It takes parameters as arguments so that it can be executed from oacis, and writes the output to a JSON file called _output.json.

`ishigami.py`


import math,json

with open('_input.json') as f:
    x = json.load(f)
x1,x2,x3 = x['x1'],x['x2'],x['x3']
y = math.sin(x1) + 7.0*(math.sin(x2))**2 + 0.1*(x2**4)*math.sin(x1)
with open('_output.json', 'w') as f:
    json.dump({"y": y}, f)

Name: Ishigami
Definition of parameters: [(x1,Float), (x2,Float), (x3,Float)]
command: "python ~/path/to/ishigami.py"

Next, prepare the following script.

from SALib.sample import saltelli
from SALib.analyze import sobol
from SALib.test_functions import Ishigami
import numpy as np

# import oacis module
import os,sys
sys.path.append( os.environ['OACIS_ROOT'] )
import oacis

# Define the model inputs
problem = {
    'num_vars': 3,
    'names': ['x1', 'x2', 'x3'],
    'bounds': [[-3.14159265359, 3.14159265359],
               [-3.14159265359, 3.14159265359],
               [-3.14159265359, 3.14159265359]]
}
# define problems
param_values = saltelli.sample(problem, 1000)
param_values

# find ishigami simulator and create ParameterSets and Runs
sim = oacis.Simulator.find_by_name('ishigami')
host = oacis.Host.find_by_name('localhost')
ps_array = []
for param in param_values:
    ps = sim.find_or_create_parameter_set( {"x1":param[0],"x2":param[1],"x3":param[2]} )
    ps.find_or_create_runs_upto( 1, submitted_to=host )
    ps_array.append(ps)

This will make a PS. After that, oacis will submit the job without permission, so wait until it is completed. (Here, we use a mathematical model that completes the calculation in an instant, but since overhead submits each one as a job, this process also takes about half a day. However, a mathematical model that takes several minutes or more. It is very useful if you want to run in parallel on a cluster etc., because such overhead can be ignored)

After completion, if you get the result as follows, the process will be exactly the same.

# output values
Y = [ ps.average_result("y")[0] for ps in ps_array ]
Y = np.array(Y)
Y

Si = sobol.analyze(problem, Y)
Si

Click here for the jupyter notebook of all the codes used this time

https://gist.github.com/yohm/e22be8bbdff338eb90849c3f1507769a

[PYTHON] Global sensitivity analysis using sensitivity analysis libraries salib and oacis

What is Sensitivity analysis?

Sensitivity analysis method

How to use SALib

Batch execution using OACIS

ishigami.py

`ishigami.py`