# *Sobol'* sensitivity indices

In this example we are going to quantify the correlation between the input variables and the output variable of a model thanks to the *Sobol'* indices.

The *Sobol'* indices allow to evaluate the importance of a single variable or a specific set of variables. 

*Sobol'* indices range is $\left[0; 1\right]$ ; the more the indice value is close to 1 the more the variable is important toward the output variance of the function. The *Sobol'* indices can be computed at different orders.

The first order indices evaluate the importance of one variable at a time ($d$ indices, with $d$ the input dimension of the model). The $d$ total indices give the relative importance of every variables except the variable $x_i$, for every variable.

Here the *Sobol'* indices are estimated on an analytical function: *Ishigami*. It writes

$$ F(\mathbf{x}) = \sin(x_1)+7\sin(x_2)^2+0.1x_3^4\sin(x_1), \quad \mathbf{x}\in [-\pi, \pi]^3 $$
Analytical values of *Sobol'* indices for this function are available:

\begin{align}
    S_{1, 2, 3} &= [0.314, 0.442, 0.], \\
    S_{T_{1, 2, 3}} &= [0.558, 0.442, 0.244].
\end{align}

This function is interesting because it exhibits an interaction between $x_1$ and $x_3$ although $x_3$ by itself do not play a role at first order.

## The Ishigami function


```python
import openturns as ot
import numpy as np
```


```python
# Create the model and input distribution
formula = ['sin(X1)+7*sin(X2)^2+0.1*X3^4*sin(X1)']
input_names = ['X1', 'X2', 'X3']
dimension = 3
corners = [[-np.pi] * dimension, [np.pi] * dimension]
model = ot.SymbolicFunction(input_names, formula)
distribution = ot.ComposedDistribution([ot.Uniform(corners[0][i], corners[1][i])
                                        for i in range(dimension)])
true_sobol = [[0.314, 0.442, 0.], [0.558, 0.442, 0.244]]  # [first, total]
```

## Without Surrogate


```python
# Create X/Y data
ot.RandomGenerator.SetSeed(0)
size = 10000
sample = ot.SobolIndicesExperiment(distribution, size, True).generate()
output = model(sample)
```


```python
# Compute Sobol' indices using the Saltelli estimator
## first order indices
sensitivityAnalysis = ot.SaltelliSensitivityAlgorithm(sample, output, size)
first_order = sensitivityAnalysis.getFirstOrderIndices()
print(f"First order Sobol' indices: {first_order}\n"
      f"Relative error: {abs(np.array(first_order) - true_sobol[0]) * 100}%")  # maximum is 1

## total order indices
total_order = sensitivityAnalysis.getTotalOrderIndices()
print(f"Total order Sobol' indices: {total_order}\n"
      f"Relative error: {abs(np.array(total_order) - true_sobol[1]) * 100}%")  # maximum is 1
```

    First order Sobol' indices: [0.302751,0.460825,0.00669407]
    Relative error: [1.12489872 1.88252502 0.66940742]%
    Total order Sobol' indices: [0.57499,0.427147,0.256687]
    Relative error: [1.69897217 1.48529071 1.26868602]%


## With Surrogate


```python
# Create a Gaussian process surrogate model
from otsurrogate import SurrogateModel
surrogate = SurrogateModel('kriging', corners, input_names)
```


```python
# Generation of data to fit the surrogate model
sequence_type = ot.LowDiscrepancyExperiment(ot.SobolSequence(), distribution, 256)
learning_sample = sequence_type.generate()
learning_output = model(learning_sample)
```


```python
surrogate.fit(learning_sample, learning_output)
```


```python
# Compute sensitivity indices
## first order indices
output, _ = surrogate(sample)  # Do not return the information about the variance
sensitivityAnalysis = ot.SaltelliSensitivityAlgorithm(sample, output, size)
first_order = sensitivityAnalysis.getFirstOrderIndices()
print(f"First order Sobol' indices: {first_order}\n"
      f"Relative error: {abs(np.array(first_order) - true_sobol[0]) * 100}%")  # maximum is 1

## total order indices
total_order = sensitivityAnalysis.getTotalOrderIndices()
print(f"Total order Sobol' indices: {total_order}\n"
      f"Relative error: {abs(np.array(total_order) - true_sobol[1]) * 100}%")  # maximum is 1
```

    First order Sobol' indices: [0.312765,0.472573,0.0070339]
    Relative error: [0.1235497  3.05730132 0.70338973]%
    Total order Sobol' indices: [0.562337,0.447125,0.237793]
    Relative error: [0.43372069 0.51250457 0.62068279]%


## Using otsklearn


```python
surrogate = SurrogateModel('otsklearn.Kriging()', corners, input_names)
surrogate.fit(learning_sample, learning_output)
```


```python
# Compute sensitivity indices
## first order indices
#output = surrogate.evaluate(sample)
output, _ = surrogate(sample)  # Do not return the information about the variance
sensitivityAnalysis = ot.SaltelliSensitivityAlgorithm(sample, output, size)
first_order = sensitivityAnalysis.getFirstOrderIndices()
print(f"First order Sobol' indices: {first_order}\n"
      f"Relative error: {abs(np.array(first_order) - true_sobol[0]) * 100}%")  # maximum is 1

## total order indices
total_order = sensitivityAnalysis.getTotalOrderIndices()
print(f"Total order Sobol' indices: {total_order}\n"
      f"Relative error: {abs(np.array(total_order) - true_sobol[1]) * 100}%")  # maximum is 1
```

    First order Sobol' indices: [0.303545,0.459592,0.00691174]
    Relative error: [1.04549598 1.75916887 0.69117411]%
    Total order Sobol' indices: [0.577597,0.426664,0.257912]
    Relative error: [1.95973662 1.53359806 1.3912165 ]%



```python

```

7d059de69ee1490032e15076773fe8ad9ba1e43a
ipynb
Jupyter Notebook
doc/examples/sensitivity_sobol.ipynb
mbaudin47/otsurrogate
2c2f043f456b8e8a752a38980c6c1515c0c90fd5
MIT
doc/examples/sensitivity_sobol.ipynb
mbaudin47/otsurrogate
2c2f043f456b8e8a752a38980c6c1515c0c90fd5
MIT
2019-02-18T14:57:50.000Z
2019-03-13T16:18:56.000Z
doc/examples/sensitivity_sobol.ipynb
tupui/otsurrogate
2c2f043f456b8e8a752a38980c6c1515c0c90fd5
MIT
2022-03-29T08:21:58.000Z
2022-03-29T08:21:58.000Z
Qwen/Qwen-72B
1. YES
2. YES