Skip to content

Commit

Permalink
Improved documentation for CramervonMises 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 21c0289 commit 86dcd1a
Show file tree
Hide file tree
Showing 4 changed files with 59 additions and 16 deletions.
10 changes: 10 additions & 0 deletions docs/code/sensitivity/cramer_von_mises/README.rst
Original file line number Diff line number Diff line change
@@ -1,3 +1,13 @@
Cramér-von Mises Sensitivity indices
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
These examples serve as a guide for using the Cramér-von Mises sensitivity module. They have been taken from various papers to enable validation of the implementation and have been referenced accordingly.

1. **Exponential function**

For the Exponential model, analytical Cramér-von Mises indices are available [1]_.

2. **Sobol function**

The Cramér-von Mises indices are computed using the Pick and Freeze approach [1]_. These model evaluations can be used to estimate the Sobol indices as well. We demonstrate this using the Sobol function.

.. [1] Gamboa, F., Klein, T., & Lagnoux, A. (2018). Sensitivity Analysis Based on Cramér-von Mises Distance. SIAM/ASA Journal on Uncertainty Quantification, 6(2), 522-548. doi:10.1137/15M1025621. (`Link <https://doi.org/10.1137/15M1025621>`_)
32 changes: 26 additions & 6 deletions docs/code/sensitivity/cramer_von_mises/plot_cvm_exponential.py
Original file line number Diff line number Diff line change
Expand Up @@ -3,9 +3,16 @@
Exponential function
==============================================
The exponential function was used in [1]_ to demonstrate the
Cramér-von Mises indices.
.. math::
f(x) := \exp(x_1 + 2x_2), \quad x_1, x_2 \sim \mathcal{N}(0, 1)
.. [1] Gamboa, F., Klein, T., & Lagnoux, A. (2018). Sensitivity Analysis Based on \
Cramér-von Mises Distance. SIAM/ASA Journal on Uncertainty Quantification, 6(2), \
522-548. doi:10.1137/15M1025621. (`Link <https://doi.org/10.1137/15M1025621>`_)
"""

# %%
Expand All @@ -15,7 +22,9 @@
from UQpy.distributions.collection.JointIndependent import JointIndependent
from UQpy.sensitivity.cramer_von_mises import CramervonMises as cvm

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

# Create Model object
model = PythonModel(
model_script="local_exponential.py",
Expand All @@ -30,29 +39,40 @@
dist_object = JointIndependent([Normal(0, 1)] * 2)

# %% [markdown]
# Compute Cramer-von Mises indices
# **Compute Cramér-von Mises indices**

# %%
# create cvm object
SA = cvm(runmodel_obj, dist_object)

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

# %% [markdown]
# Cramer-von Mises sensitivity analysis
# **Cramér-von Mises indices**
#
# Expected value of the sensitivity indices:
#
# $S^1_{CVM} = \frac{6}{\pi} \operatorname{arctan}(2) - 2 \approx 0.1145$
# :math:`S^1_{CVM} = \frac{6}{\pi} \operatorname{arctan}(2) - 2 \approx 0.1145`
#
# $S^2_{CVM} = \frac{6}{\pi} \operatorname{arctan}(\sqrt{19}) - 2 \approx 0.5693$
# :math:`S^2_{CVM} = \frac{6}{\pi} \operatorname{arctan}(\sqrt{19}) - 2 \approx 0.5693`

# %%
computed_indices["CVM_i"]

# %% [markdown]
# **Estimated first order Sobol indices**
#
# Expected first order Sobol indices:
#
# :math:`S_1` = 0.0118
#
# :math:`S_2` = 0.3738

# %%
computed_indices["sobol_i"]

# %% [markdown]
# **Estimated total order Sobol indices**

# %%
computed_indices["sobol_total_i"]
31 changes: 22 additions & 9 deletions docs/code/sensitivity/cramer_von_mises/plot_cvm_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,10 +27,12 @@
from UQpy.distributions.collection.JointIndependent import JointIndependent
from UQpy.sensitivity.cramer_von_mises import CramervonMises as cvm

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

# Create Model object
num_vars = 6
a_vals = np.array([1.0, 2.0, 3.0, 4.0, 5.0, 6.0])
a_vals = np.array([0.0, 0.5, 3.0, 9.0, 99.0, 99.0])

model = PythonModel(
model_script="local_sobol_func.py",
Expand All @@ -41,30 +47,37 @@
# Define distribution object
dist_object = JointIndependent([Uniform(0, 1)] * num_vars)

# %% [markdown]
# **Compute Cramér-von Mises indices**

# %%
SA = cvm(runmodel_obj, dist_object)

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

# %% [markdown]
# **Cramér-von Mises indices**

# %%
computed_indices["CVM_i"]

# %% [markdown]
# Sobol indices computed analytically
# **Estimated Sobol indices**
#
# $S_1$ = 0.46067666
# Expected first order Sobol indices:
#
# $S_2$ = 0.20474518
# :math:`S_1` = 5.86781190e-01
#
# $S_3$ = 0.11516917
# :math:`S_2` = 2.60791640e-01
#
# $S_4$ = 0.07370827
# :math:`S_3` = 3.66738244e-02
#
# $S_5$ = 0.0511863
# :math:`S_4` = 5.86781190e-03
#
# $S_6$ = 0.03760626
# :math:`S_5` = 5.86781190e-05
#
# :math:`S_6` = 5.86781190e-05

# %%
computed_indices["sobol_i"]
2 changes: 1 addition & 1 deletion docs/source/sensitivity/cramer_von_mises.rst
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
Cramér-von Mises indices
----------------------------------------

A sensitivity index based on the Cramér-von Mises distance. In contrast to variance based Sobol indices it takes into account the whole distribution of the model output and is therefore considered as a moment-free method [1]_. Furthermore the index can be naturally extended to multivariate model outputs (not implemented yet in UQPy).
A sensitivity index based on the Cramér-von Mises distance. In contrast to the variance based Sobol indices, it takes into account the whole distribution of the model output and is therefore considered as a moment-free method [1]_. Furthermore the index can be naturally extended to multivariate model outputs (not implemented yet in UQPy).

Consider a model :math:`Y=f(X): \mathbb{R}^d \rightarrow \mathbb{R}^k` with :math:`d` inputs :math:`X_{(1)}, X_{(2)}, \ldots, X_{(d)}` and :math:`k` outputs :math:`Y_{(1)}, Y_{(2)}, \ldots, Y_{(k)}`. We define the cumulative distribution function :math:`F(t)` of :math:`Y` as:

Expand Down

0 comments on commit 86dcd1a

Please sign in to comment.