# Sobol Indices¶

It is required to first build a POD object based on the Kriging metamodel or on the polynomial chaos in order to compute the Sobol indices. It also can be used only if the input parameters dimension is greater than 2 (without counting the defect).

```
In [1]:
```

```
# import relevant module
import openturns as ot
import otpod
# enable display figure in notebook
try:
%matplotlib inline
except:
pass
```

## Generate data¶

```
In [2]:
```

```
inputSample = ot.NumericalSample(
[[4.59626812e+00, 7.46143339e-02, 1.02231538e+00, 8.60042277e+01],
[4.14315790e+00, 4.20801346e-02, 1.05874908e+00, 2.65757364e+01],
[4.76735111e+00, 3.72414824e-02, 1.05730385e+00, 5.76058433e+01],
[4.82811977e+00, 2.49997658e-02, 1.06954641e+00, 2.54461380e+01],
[4.48961094e+00, 3.74562922e-02, 1.04943946e+00, 6.19483646e+00],
[5.05605334e+00, 4.87599783e-02, 1.06520409e+00, 3.39024904e+00],
[5.69679328e+00, 7.74915877e-02, 1.04099514e+00, 6.50990466e+01],
[5.10193991e+00, 4.35520544e-02, 1.02502536e+00, 5.51492592e+01],
[4.04791970e+00, 2.38565932e-02, 1.01906882e+00, 2.07875350e+01],
[4.66238956e+00, 5.49901237e-02, 1.02427200e+00, 1.45661275e+01],
[4.86634219e+00, 6.04693570e-02, 1.08199374e+00, 1.05104730e+00],
[4.13519347e+00, 4.45225831e-02, 1.01900124e+00, 5.10117047e+01],
[4.92541940e+00, 7.87692335e-02, 9.91868726e-01, 8.32302238e+01],
[4.70722074e+00, 6.51799251e-02, 1.10608515e+00, 3.30181002e+01],
[4.29040932e+00, 1.75426222e-02, 9.75678838e-01, 2.28186756e+01],
[4.89291400e+00, 2.34997929e-02, 1.07669835e+00, 5.38926138e+01],
[4.44653744e+00, 7.63175936e-02, 1.06979154e+00, 5.19109415e+01],
[3.99977452e+00, 5.80430585e-02, 1.01850716e+00, 7.61988190e+01],
[3.95491570e+00, 1.09302814e-02, 1.03687664e+00, 6.09981789e+01],
[5.16424368e+00, 2.69026464e-02, 1.06673711e+00, 2.88708887e+01],
[5.30491620e+00, 4.53802273e-02, 1.06254792e+00, 3.03856837e+01],
[4.92809155e+00, 1.20616369e-02, 1.00700410e+00, 7.02512744e+00],
[4.68373805e+00, 6.26028935e-02, 1.05152117e+00, 4.81271603e+01],
[5.32381954e+00, 4.33013582e-02, 9.90522007e-01, 6.56015973e+01],
[4.35455857e+00, 1.23814619e-02, 1.01810539e+00, 1.10769534e+01]])
signals = ot.NumericalSample(
[[ 37.305445], [ 35.466919], [ 43.187991], [ 45.305165], [ 40.121222], [ 44.609524],
[ 45.14552 ], [ 44.80595 ], [ 35.414039], [ 39.851778], [ 42.046049], [ 34.73469 ],
[ 39.339349], [ 40.384559], [ 38.718623], [ 46.189709], [ 36.155737], [ 31.768369],
[ 35.384313], [ 47.914584], [ 46.758537], [ 46.564428], [ 39.698493], [ 45.636588],
[ 40.643948]])
```

```
In [3]:
```

```
# signal detection threshold
detection = 38.
```

## Build POD with Kriging model¶

### Running the Kriging based POD¶

```
In [4]:
```

```
krigingPOD = otpod.KrigingPOD(inputSample, signals, detection)
# we can change all simulation size parameters as we are not interested in having an accurate POD curve
krigingPOD.setSamplingSize(200)
krigingPOD.setSimulationSize(50)
krigingPOD.run()
```

```
Start optimizing covariance model parameters...
Kriging optimizer completed
kriging validation Q2 (>0.9): 1.0000
Computing POD per defect: [==================================================] 100.00% Done
```

### Show POD graphs¶

```
In [5]:
```

```
fig, ax = krigingPOD.drawPOD(probabilityLevel=0.9, confidenceLevel=0.95,
name='figure/PODKriging.png')
# The figure is saved in PODPolyChaos.png
fig.show()
```

```
/home/dumas/anaconda2/lib/python2.7/site-packages/matplotlib/figure.py:397: UserWarning: matplotlib is currently using a non-GUI backend, so cannot show the figure
"matplotlib is currently using a non-GUI backend, "
```

## Build POD with polynomial chaos model¶

### Running the chaos based POD¶

```
In [6]:
```

```
chaosPOD = otpod.PolynomialChaosPOD(inputSample, signals, detection)
# we can change all simulation size parameters as we are not interested in having an accurate POD curve
chaosPOD.setSamplingSize(200)
chaosPOD.setSimulationSize(50)
chaosPOD.run()
```

```
Start build polynomial chaos model...
Polynomial chaos model completed
Polynomial chaos validation R2 (>0.8) : 0.9999
Polynomial chaos validation Q2 (>0.8) : 0.9987
Computing POD per defect: [==================================================] 100.00% Done
```

### Show POD graphs¶

```
In [7]:
```

```
fig, ax = chaosPOD.drawPOD(probabilityLevel=0.9, confidenceLevel=0.95,
name='figure/PODChaos.png')
# The figure is saved in PODPolyChaos.png
fig.show()
```

## Run the sensitivity analysis¶

The sensitivity analysis can only be performed with POD computed with a kriging metamodel or a polynomial chaos.

The Sobol indices are aggregated indices computed for the defect sizes defined in the POD study.

## Using the kriging model¶

```
In [8]:
```

```
# number of simulations
N = 1000
sobol = otpod.SobolIndices(krigingPOD, N)
sobol.run()
```

### Draw the figure with given labels¶

The default labels are \(Xi\) but the user can specify its own input
labels. Besides, the figure can be saved specifying the attribute
*name*.

```
In [9]:
```

```
label = ['X1', 'variable2', 'param3']
```

```
In [10]:
```

```
fig, ax = sobol.drawAggregatedIndices(label, name='figure/Sobol.png')
fig.show()
```

```
In [11]:
```

```
fig, ax = sobol.drawFirstOrderIndices(label, name='figure/FirstOrderSobol.png')
fig.show()
```

```
In [12]:
```

```
fig, ax = sobol.drawTotalOrderIndices(label, name='figure/TotalOrderSobol.png')
fig.show()
```

### Get the numerical results¶

The Sobol indices are given in the OpenTURNS object
SobolIndicesAlgorithm.
The method *getSensitivityResult* allows the get this object and then
get back all wanted results.

```
In [13]:
```

```
# get the OpenTURNS object
sobol_result = sobol.getSensitivityResult()
# get aggregated indices
print("Aggregated first order: {}".format(sobol_result.getAggregatedFirstOrderIndices()))
print("Aggregated total order: {}".format(sobol_result.getAggregatedTotalOrderIndices()))
# get the confidence interval
print('\nFirst order confidence interval:')
print(sobol_result.getFirstOrderIndicesInterval())
print('\nTotal order confidence interval:')
print(sobol_result.getTotalOrderIndicesInterval())
```

```
Aggregated first order: [0.772678,0.0317259,0.000271994]
Aggregated total order: [0.968083,0.226075,0.00066017]
First order confidence interval:
[0.719211, 0.829002]
[0.0158451, 0.0448602]
[-0.000689527, 0.00120284]
Total order confidence interval:
[0.901833, 1.02883]
[0.171566, 0.269238]
[-0.0438167, 0.0374588]
```

It is also possible to retreive the Sobol indices for one defect size among the list.

As example, we want the indices for the 4th defect size in the list. It
may **return an error** if the indices cannot be computed because no
variability exists. It is the case when the POD is equal to 0 or 1.

```
In [14]:
```

```
print("Defect sizes: {}".format(sobol.getDefectSizes()))
# get indices for a specific defect sizes
i = 3 # correspond with the 4th value.
print("\nIndices for defect size {:.3f}: {}".format(sobol.getDefectSizes()[i], sobol_result.getFirstOrderIndices(i)))
```

```
Defect sizes: [ 3.9549157 4.04659347 4.13827123 4.229949 4.32162677 4.41330454
4.5049823 4.59666007 4.68833784 4.78001561 4.87169337 4.96337114
5.05504891 5.14672668 5.23840444 5.33008221 5.42175998 5.51343775
5.60511551 5.69679328]
Indices for defect size 4.230: [0.727317,0.0140776,0.000158337]
```

### Change the defect sizes list¶

It is possible to modify the list of the defect sizes either to reduce the range or to compute the indices for a specific defect value. If only one defect size is provided, then the aggregated indices correspond to the indices.

```
In [15]:
```

```
sobol.setDefectSizes([4.5])
sobol.run()
```

```
In [16]:
```

```
# get the OpenTURNS object
sobol_result = sobol.getSensitivityResult()
# get aggregated indices
print("Aggregated first order: {}".format(sobol_result.getAggregatedFirstOrderIndices()))
print("Aggregated total order: {}".format(sobol_result.getAggregatedTotalOrderIndices()))
# get indices
print("First order: {}".format(sobol_result.getFirstOrderIndices()))
print("Total order: {}".format(sobol_result.getTotalOrderIndices()))
```

```
Aggregated first order: [0.851168,0.0515212,0.000266685]
Aggregated total order: [1.02607,0.159002,0.000775292]
First order: [0.851168,0.0515212,0.000266685]
Total order: [1.02607,0.159002,0.000775292]
```

### Change the method to compute the indices¶

OpenTURNS implements 4 methods : Saltelli, Martinez, Jansen and
Mauntz-Kucherenko. These methods can be chosen using the method
*setSensitivityMethod*.

```
In [17]:
```

```
sobol.setSensitivityMethod("Martinez")
sobol.run()
```

```
In [18]:
```

```
# get the OpenTURNS object
sobol_result = sobol.getSensitivityResult()
# get aggregated indices
print("Aggregated first order: {}".format(sobol_result.getAggregatedFirstOrderIndices()))
print("Aggregated total order: {}".format(sobol_result.getAggregatedTotalOrderIndices()))
```

```
Aggregated first order: [0.858261,0.0297892,-0.0040736]
Aggregated total order: [0.982881,0.153882,0.000308148]
```

## Case with polynomial chaos¶

With polynomial chaos, the POD is computed simulating several polynomial
chaos coefficients. Then it requires more times than with Kriging. The
number of simulations is initially set to 1000 but it can be changed
using the method *setSimulationSize*.

```
In [19]:
```

```
# number of simulations
N = 1000
sobol2 = otpod.SobolIndices(chaosPOD, N)
#sobol2.setSimulationSize(500)
sobol2.run()
```

```
In [20]:
```

```
# get the OpenTURNS object
sobol_result2 = sobol2.getSensitivityResult()
# get aggregated indices
print("Aggregated first order: {}".format(sobol_result2.getAggregatedFirstOrderIndices()))
print("Aggregated total order: {}".format(sobol_result2.getAggregatedTotalOrderIndices()))
# get the confidence interval
print('\nFirst order confidence interval:')
print(sobol_result2.getFirstOrderIndicesInterval())
print('\nTotal order confidence interval:')
print(sobol_result2.getTotalOrderIndicesInterval())
```

```
Aggregated first order: [0.733816,0.0303936,0.00202081]
Aggregated total order: [0.912059,0.221395,0.0131351]
First order confidence interval:
[0.674632, 0.779876]
[0.011269, 0.0475021]
[-0.00224056, 0.00590774]
Total order confidence interval:
[0.860978, 0.966444]
[0.176179, 0.263411]
[-0.0289697, 0.0499621]
```

```
In [21]:
```

```
fig, ax = sobol2.drawAggregatedIndices()
fig.show()
```