Skip to content

Commit

Permalink
Added examples Cramér-von Mises sensitivity
Browse files Browse the repository at this point in the history
  • Loading branch information
Prateek Bhustali committed May 8, 2022
1 parent 27305fc commit 4d4a012
Show file tree
Hide file tree
Showing 5 changed files with 192 additions and 0 deletions.
20 changes: 20 additions & 0 deletions docs/code/sensitivity/cramer_von_mises/local_exponential.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,20 @@
"""
Auxiliary file
==============================================
"""

import numpy as np


def evaluate(X: np.array) -> np.array:
r"""A non-linear function that is used to demonstrate sensitivity index.
.. math::
f(x) = \exp(x_1 + 2*x_2)
"""

Y = np.exp(X[:, 0] + 2 * X[:, 1])

return Y
42 changes: 42 additions & 0 deletions docs/code/sensitivity/cramer_von_mises/local_sobol_func.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,42 @@
"""
Auxiliary file
==============================================
"""

import numpy as np
import copy


def evaluate(X, a_values):

dims = len(a_values)
g = 1

for i in range(dims):
g_i = (np.abs(4 * X[:, i] - 2) + a_values[i]) / (1 + a_values[i])
g *= g_i

return g


def sensitivities(a_values):

dims = len(a_values)

Total_order = np.zeros((dims, 1))

V_i = (3 * (1 + a_values) ** 2) ** (-1)

total_variance = np.prod(1 + V_i) - 1

First_order = V_i / total_variance

for i in range(dims):

rem_First_order = copy.deepcopy(V_i)
rem_First_order[i] = 0
Total_order[i] = V_i[i] * np.prod(rem_First_order + 1) / total_variance

return First_order.reshape(-1, 1), Total_order
58 changes: 58 additions & 0 deletions docs/code/sensitivity/cramer_von_mises/plot_cvm_exponential.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,58 @@
"""
Exponential function
==============================================
.. math::
f(x) := \exp(x_1 + 2x_2), \quad x_1, x_2 \sim \mathcal{N}(0, 1)
"""

# %%
from UQpy.run_model.RunModel import RunModel
from UQpy.run_model.model_execution.PythonModel import PythonModel
from UQpy.distributions import Normal
from UQpy.distributions.collection.JointIndependent import JointIndependent
from UQpy.sensitivity.cramer_von_mises import CramervonMises as cvm

# %%
# Create Model object
model = PythonModel(
model_script="local_exponential.py",
model_object_name="evaluate",
var_names=[r"$X_1$", "$X_2$"],
delete_files=True,
)

runmodel_obj = RunModel(model=model)

# Define distribution object
dist_object = JointIndependent([Normal(0, 1)] * 2)

# %% [markdown]
# Compute Cramer-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
#
# Expected value of the sensitivity indices:
#
# $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$

# %%
computed_indices["CVM_i"]

# %%
computed_indices["sobol_i"]

# %%
computed_indices["sobol_total_i"]
70 changes: 70 additions & 0 deletions docs/code/sensitivity/cramer_von_mises/plot_cvm_sobol_func.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,70 @@
r"""
Sobol function
==============================================
.. math::
g(x_1, x_2, \ldots, x_D) := \prod_{i=1}^{D} \frac{|4x_i - 2| + a_i}{1 + a_i},
where,
.. math::
x_i \sim \mathcal{U}(0, 1), \quad a_i \in \mathbb{R}.
"""

# %%
import numpy as np

from UQpy.run_model.RunModel import RunModel
from UQpy.run_model.model_execution.PythonModel import PythonModel
from UQpy.distributions import Uniform
from UQpy.distributions.collection.JointIndependent import JointIndependent
from UQpy.sensitivity.cramer_von_mises import CramervonMises as cvm

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

model = PythonModel(
model_script="local_sobol_func.py",
model_object_name="evaluate",
var_names=[r"$X_1$", "$X_2$"],
delete_files=True,
a_values=a_vals,
)

runmodel_obj = RunModel(model=model)

# Define distribution object
dist_object = JointIndependent([Uniform(0, 1)] * num_vars)

# %%
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)

# %%
computed_indices["CVM_i"]

# %% [markdown]
# Sobol indices computed analytically
#
# $S_1$ = 0.46067666
#
# $S_2$ = 0.20474518
#
# $S_3$ = 0.11516917
#
# $S_4$ = 0.07370827
#
# $S_5$ = 0.0511863
#
# $S_6$ = 0.03760626
#

# %%
computed_indices["sobol_i"]
2 changes: 2 additions & 0 deletions docs/source/conf.py
Original file line number Diff line number Diff line change
Expand Up @@ -92,6 +92,7 @@
"../code/transformations/nataf",
"../code/sensitivity/morris",
"../code/sensitivity/sobol",
"../code/sensitivity/cramer_von_mises",
"../code/stochastic_processes/bispectral",
"../code/stochastic_processes/karhunen_loeve",
"../code/stochastic_processes/spectral",
Expand Down Expand Up @@ -127,6 +128,7 @@
"auto_examples/transformations/nataf",
"auto_examples/sensitivity/morris",
"auto_examples/sensitivity/sobol",
"auto_examples/sensitivity/cramer_von_mises",
"auto_examples/stochastic_processes/bispectral",
"auto_examples/stochastic_processes/karhunen_loeve",
"auto_examples/stochastic_processes/spectral",
Expand Down

0 comments on commit 4d4a012

Please sign in to comment.