Skip to content

Commit

Permalink
Improved documentation Sobol indices
Browse files Browse the repository at this point in the history
+ Fomatting
+ Detailed descriptions
  • Loading branch information
Prateek Bhustali committed May 9, 2022
1 parent 26d27db commit d842ec6
Show file tree
Hide file tree
Showing 8 changed files with 132 additions and 100 deletions.
23 changes: 17 additions & 6 deletions docs/code/sensitivity/sobol/README.rst
Original file line number Diff line number Diff line change
Expand Up @@ -7,18 +7,29 @@ Single output models
======================
We demonstrate the computation of the Sobol indices for models with a single output using the following examples:

1. Ishigami function
2. Exponential function
3. Sobol function with parameters a := [0., 0.5, 3., 9., 99., 99.] : Example from [2] page 11
1. **Additive function**

This is a beginner-friendly example for introducing Sobol indices. The function is a linear combination of two inputs which produces a scalar output.

2. **Ishigami function**

The Ishigami function is a non-linear, non-monotonic function that is commonly used to benchmark uncertainty and senstivity analysis methods.

3. **Sobol function**

The Sobol function is non-linear function that is commonly used to benchmark uncertainty
and senstivity analysis methods. Unlike the Ishigami function which has 3 input
variables, the Sobol function can have any number of input variables (see [2]_).

Multiple output models
========================

We demonstrate the computation of the Sobol indices for models with multiple outputs using the following example:

1. Mechanical oscillator ODE (numerical model): Example from [1] page 19
1. **Mechanical oscillator ODE**

The Sobol indices are computed for a mechanical oscillator governed by a second-order differential equation [1]_. The model outputs the displacement of the oscillator for a given time period. Here the sensitivity of the model parameters are computed at each point in time (see [1]_).

[1] Gamboa F, Janon A, Klein T, Lagnoux A, others. Sensitivity analysis for multidimensional and functional outputs. Electronic journal of statistics 2014; 8(1): 575-603.
.. [1] Gamboa F, Janon A, Klein T, Lagnoux A, others. Sensitivity analysis for multidimensional and functional outputs. Electronic journal of statistics 2014; 8(1): 575-603.
[2] Saltelli, A. (2002). Making best use of model evaluations to compute indices.
.. [2] Saltelli, A. (2002). Making best use of model evaluations to compute indices.
21 changes: 21 additions & 0 deletions docs/code/sensitivity/sobol/local_additive.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,21 @@
"""
Auxiliary file
==============================================
"""

import numpy as np


def evaluate(X, params) -> np.array:
r"""A linear function that is used to demonstrate sensitivity indices.
.. math::
f(x) = a \cdot x_1 + b \cdot x_2
"""
a, b = params

Y = a * X[:, 0] + b * X[:, 1]

return Y
20 changes: 0 additions & 20 deletions docs/code/sensitivity/sobol/local_exponential.py

This file was deleted.

16 changes: 13 additions & 3 deletions docs/code/sensitivity/sobol/plot_mechanical_oscillator_ODE.py
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,10 @@
.. math::
m \sim \mathcal{U}(10, 12), c \sim \mathcal{U}(0.4, 0.8), k \sim \mathcal{U}(70, 90), \ell \sim \mathcal{U}(-1, -0.25).
Here, we compute the Sobol indices for each point in time and are called
pointwise-in-time Sobol indices. These indices describe the sensitivity of the model
parameters at each point in time.
"""

# %%
Expand All @@ -28,7 +32,9 @@
from UQpy.distributions.collection.JointIndependent import JointIndependent
from UQpy.sensitivity.sobol import Sobol

# %%
# %% [markdown]
# **Define the model and input distributions**

# Create Model object
model = PythonModel(
model_script="local_mechanical_oscillator_ODE.py",
Expand All @@ -46,13 +52,17 @@
L = Uniform(-1, (-0.25 - -1))
dist_object = JointIndependent([M, C, K, L])

# %%
# %% [markdown]
# **Compute Sobol indices**

# %% [markdown]
SA = Sobol(runmodel_obj, dist_object)

computed_indices = SA.run(n_samples=500)

# %%
# Plot the Sobol indices
# **Plot the Sobol indices**

t_0 = 0
t_f = 40
dt = 0.05
Expand Down
Original file line number Diff line number Diff line change
@@ -1,10 +1,10 @@
"""
Exponential function
Additive function
==============================================
.. math::
f(x) := \exp(x_1 + 2x_2), \quad x_1, x_2 \sim \mathcal{N}(0, 1)
f(x) = a \cdot X_1 + b \cdot X_2, \quad X_1, X_2 \sim \mathcal{N}(0, 1), \quad a,b \in \mathbb{R}
"""

Expand All @@ -15,16 +15,21 @@
from UQpy.distributions.collection.JointIndependent import JointIndependent
from UQpy.sensitivity.sobol import Sobol

# %%
# %% [markdown]
# **Define the model and input distributions**

# Create Model object
a, b = 1, 2

model = PythonModel(
model_script="local_exponential.py",
model_script="local_additive.py",
model_object_name="evaluate",
var_names=[
"X_1",
"X_2",
],
delete_files=True,
params=[a, b],
)

runmodel_obj = RunModel(model=model)
Expand All @@ -33,28 +38,21 @@
dist_object = JointIndependent([Normal(0, 1)] * 2)

# %% [markdown]
# Compute Sobol indices
# **Compute Sobol indices**

# %%
# %% [markdown]
SA = Sobol(runmodel_obj, dist_object)

# Compute Sobol indices using the pick and freeze algorithm
computed_indices = SA.run(
n_samples=100_000, num_bootstrap_samples=1_000, confidence_level=0.95
)
computed_indices = SA.run(n_samples=50_000)

# %% [markdown]
# Expected first order Sobol indices (computed analytically):
# **First order Sobol indices**
#
# Expected first order Sobol indices:
#
# X1: 0.0118
# :math:`\mathrm{S}_1 = \frac{a^2 \cdot \mathbb{V}[X_1]}{a^2 \cdot \mathbb{V}[X_1] + b^2 \cdot \mathbb{V}[X_2]} = \frac{1^2 \cdot 1}{1^2 \cdot 1 + 2^2 \cdot 1} = 0.2`
#
# X2: 0.3738
# :math:`\mathrm{S}_2 = \frac{b^2 \cdot \mathbb{V}[X_2]}{a^2 \cdot \mathbb{V}[X_1] + b^2 \cdot \mathbb{V}[X_2]} = \frac{2^2 \cdot 1}{1^2 \cdot 1 + 2^2 \cdot 1} = 0.8`

# %%
computed_indices["sobol_i"]

# %% [markdown]
# Confidence intervals for first order Sobol indices

# %%
computed_indices["CI_sobol_i"]
78 changes: 45 additions & 33 deletions docs/code/sensitivity/sobol/plot_sobol_func.py
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,10 @@
Sobol function
==============================================
The Sobol function is non-linear function that is commonly used to benchmark uncertainty
and senstivity analysis methods. Unlike the Ishigami function which has 3 input
variables, the Sobol function can have any number of input variables.
.. math::
g(x_1, x_2, \ldots, x_D) := \prod_{i=1}^{D} \frac{|4x_i - 2| + a_i}{1 + a_i},
Expand All @@ -23,7 +27,9 @@
from UQpy.distributions.collection.JointIndependent import JointIndependent
from UQpy.sensitivity.sobol import Sobol

# %%
# %% [markdown]
# **Define the model and input distributions**

# Create Model object
num_vars = 6
a_vals = np.array([0.0, 0.5, 3.0, 9.0, 99.0, 99.0])
Expand All @@ -42,83 +48,89 @@
dist_object = JointIndependent([Uniform(0, 1)] * num_vars)

# %% [markdown]
# #### Compute Sobol indices
# **Compute Sobol indices**

# %%
# %% [markdown]
SA = Sobol(runmodel_obj, dist_object)

# Compute Sobol indices using the pick and freeze algorithm
computed_indices = SA.run(n_samples=50_000, estimate_second_order=True)

# %% [markdown]
# First order Sobol indices
# **First order Sobol indices**
#
# Expected first order Sobol indices:
#
# $S_1$ = 5.86781190e-01
# :math:`S_1` = 5.86781190e-01
#
# $S_2$ = 2.60791640e-01
# :math:`S_2` = 2.60791640e-01
#
# $S_3$ = 3.66738244e-02
# :math:`S_3` = 3.66738244e-02
#
# $S_4$ = 5.86781190e-03
# :math:`S_4` = 5.86781190e-03
#
# $S_5$ = 5.86781190e-05
# :math:`S_5` = 5.86781190e-05
#
# $S_6$ = 5.86781190e-05
# :math:`S_6` = 5.86781190e-05

# %%
computed_indices["sobol_i"]

# %% [markdown]
# Total order Sobol indices
# **Total order Sobol indices**
#
# $S_{T_1}$ = 6.90085892e-01
# Expected total order Sobol indices:
#
# $S_{T_2}$ = 3.56173364e-01
# :math:`S_{T_1}` = 6.90085892e-01
#
# $S_{T_3}$ = 5.63335422e-02
# :math:`S_{T_2}` = 3.56173364e-01
#
# $S_{T_4}$ = 9.17057664e-03
# :math:`S_{T_3}` = 5.63335422e-02
#
# $S_{T_5}$ = 9.20083854e-05
# :math:`S_{T_4}` = 9.17057664e-03
#
# $S_{T_6}$ = 9.20083854e-05
# :math:`S_{T_5}` = 9.20083854e-05
#
# :math:`S_{T_6}` = 9.20083854e-05
#

# %%
computed_indices["sobol_total_i"]

# %% [markdown]
# Second-order Sobol indices
# **Second order Sobol indices**
#
# Expected second order Sobol indices:
#
# $S_{12}$ = 0.0869305
# :math:`S_{T_{12}}` = 0.0869305
#
# $S_{13}$ = 0.0122246
# :math:`S_{T_{13}}` = 0.0122246
#
# $S_{14}$ = 0.00195594
# :math:`S_{T_{14}}` = 0.00195594
#
# $S_{15}$ = 0.00001956
# :math:`S_{T_{15}}` = 0.00001956
#
# $S_{16}$ = 0.00001956
# :math:`S_{T_{16}}` = 0.00001956
#
# $S_{23}$ = 0.00543316
# :math:`S_{T_{23}}` = 0.00543316
#
# $S_{24}$ = 0.00086931
# :math:`S_{T_{24}}` = 0.00086931
#
# $S_{25}$ = 0.00000869
# :math:`S_{T_{25}}` = 0.00000869
#
# $S_{26}$ = 0.00000869
# :math:`S_{T_{26}}` = 0.00000869
#
# $S_{34}$ = 0.00012225
# :math:`S_{T_{34}}` = 0.00012225
#
# $S_{35}$ = 0.00000122
# :math:`S_{T_{35}}` = 0.00000122
#
# $S_{36}$ = 0.00000122
# :math:`S_{T_{36}}` = 0.00000122
#
# $S_{45}$ = 0.00000020
# :math:`S_{T_{45}}` = 0.00000020
#
# $S_{46}$ = 0.00000020
# :math:`S_{T_{46}}` = 0.00000020
#
# $S_{56}$ = 2.0e-9
# :math:`S_{T_{56}}` = 2.0e-9

# %%
computed_indices["sobol_ij"]
Loading

0 comments on commit d842ec6

Please sign in to comment.