Sensitivity Analysis

Quantifying the sensitivity of a model output $f$ to the model parameters $z$ can be an important component of any modeling exercise. This section demonstrates how to use popular local and global sensitivity analysis.

Sobol Indices

Any function $f$ with finite variance parameterized by a set of independent variables $z$ with $\rho (z)=\prod _{j=1}^d\rho (z_j)$ and support $\mathrm{\Gamma }=\underset{j=1}{\overset{d}{⨂}}{\mathrm{\Gamma }}_j$ can be decomposed into a finite sum, referred to as the ANOVA decomposition,

$$f(z)={\hat{f}}_0+\sum _{i=1}^d{\hat{f}}_i(z_i)+\sum _{i,j=1}^d{\hat{f}}_{i,j}(z_i,z_j)+\cdots +{\hat{f}}_{1,\dots ,d}(z_1,\dots ,z_d)$$

or more compactly

$$f(z)=\sum _{\mathit{u}\subseteq \mathcal{D}}{\hat{f}}_{\mathit{u}}(z_{\mathit{u}})$$

where ${\hat{f}}_{\mathit{u}}$ quantifies the dependence of the function $f$ on the variable dimensions $i\in \mathit{u}$ and $\mathit{u}=(u_1,\dots ,u_s)\subseteq \mathcal{D}=\{1,\dots ,d\}$.

The functions ${\hat{f}}_{\mathit{u}}$ can be obtained by integration, specifically

$${\hat{f}}_{\mathit{u}}(z_{\mathit{u}})={\int }_{{\mathrm{\Gamma }}_{\mathcal{D}\setminus \mathit{u}}}f(z)\text{d}{\rho }_{\mathcal{D}\setminus \mathit{u}}(z)-\sum _{\mathit{v}\subset \mathit{u}}{\hat{f}}_{\mathit{v}}(z_{\mathit{v}}),$$

where $\text{d}{\rho }_{\mathcal{D}\setminus \mathit{u}}(z)=\prod _{j\notin \mathit{u}}\text{d}{\rho }_j(z)$ and ${\mathrm{\Gamma }}_{\mathcal{D}\setminus \mathit{u}}=\underset{j\notin \mathit{u}}{⨂}{\mathrm{\Gamma }}_j$.

The first-order terms ${\hat{f}}_{\mathit{u}}(z_i)$, $\Vert \mathit{u}{\Vert }_0=1$ represent the effect of a single variable acting independently of all others. Similarly, the second-order terms $\Vert \mathit{u}{\Vert }_0=2$ represent the contributions of two variables acting together, and so on.

The terms of the ANOVA expansion are orthogonal, i.e. the weighted $L^2$ inner product $({\hat{f}}_{\mathit{u}},{\hat{f}}_{\mathit{v}}{)}_{L_{\rho }^2}=0$, for $\mathit{u}\ne \mathit{v}$. This orthogonality facilitates the following decomposition of the variance of the function $f$

$$\mathbb{V}\left[f\right]=\sum _{\mathit{u}\subseteq \mathcal{D}}\mathbb{V}\left[{\hat{f}}_{\mathit{u}}\right],\mathbb{V}\left[{\hat{f}}_{\mathit{u}}\right]={\int }_{{\mathrm{\Gamma }}_{\mathit{u}}}{f}_{\mathit{u}}^2\text{d}{\rho }_{\mathit{u}},$$

where $\text{d}{\rho }_{\mathit{u}}(z)=\prod _{j\in \mathit{u}}\text{d}{\rho }_j(z)$.

The quantities $\mathbb{V}\left[{\hat{f}}_{\mathit{u}}\right]/\mathbb{V}\left[f\right]$ are referred to as Sobol indices [SMCS2001] and are frequently used to estimate the sensitivity of $f$ to single, or combinations of input parameters. Note that this is a global sensitivity, reflecting a variance attribution over the range of the input parameters, as opposed to the local sensitivity reflected by a derivative. Two popular measures of sensitivity are the main effect and total effect indices given respectively by

$$S_i=\frac{\mathbb{V}\left[{\hat{f}}_{{\mathit{e}}_i}\right]}{\mathbb{V}\left[f\right]},S_i^T=\frac{\sum _{\mathit{u}\in \mathcal{J}}\mathbb{V}\left[{\hat{f}}_{\mathit{u}}\right]}{\mathbb{V}\left[f\right]}$$

where ${\mathit{e}}_i$ is the unit vector, with only one non-zero entry located at the $i$-th element, and $\mathcal{J}=\{\mathit{u}:i\in \mathit{u}\}$.

Sobol indices can be computed different ways. In the following we will use polynomial chaos expansions, as in [SRESS2008].

import matplotlib.pyplot as plt
from pyapprox.benchmarks import setup_benchmark
from pyapprox.surrogates import approximate
from pyapprox import analysis
benchmark = setup_benchmark("ishigami", a=7, b=0.1)

num_samples = 1000
train_samples = benchmark.variable.rvs(num_samples)
train_vals = benchmark.fun(train_samples)

approx_res = approximate(
    train_samples, train_vals, 'polynomial_chaos',
    {'basis_type': 'hyperbolic_cross', 'variable': benchmark.variable,
     'options': {'max_degree': 8}})
pce = approx_res.approx

res = analysis.gpc_sobol_sensitivities(pce, benchmark.variable)

Now lets compare the estimated values with the exact value

print(res.main_effects[:, 0])
print(benchmark.main_effects[:, 0])

[3.13751427e-01 4.42645977e-01 6.47617333e-08]
[0.31390519 0.44241114 0.        ]

We can visualize the sensitivity indices using the following

fig, axs = plt.subplots(1, 3, figsize=(3*8, 6))
analysis.plot_main_effects(benchmark.main_effects, axs[0])
analysis.plot_total_effects(benchmark.total_effects, axs[1])
analysis.plot_interaction_values(benchmark.sobol_indices,
                                 benchmark.sobol_interaction_indices, axs[2])
axs[0].set_title(r'$\mathrm{Main\;Effects}$')
axs[1].set_title(r'$\mathrm{Total\;Effects}$')
axs[2].set_title(r'$\mathrm{Sobol\;Indices}$')
plt.show()

References

[SMCS2001]

I.M. Sobol. Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates. Mathematics and Computers in Simulation, 55(3): 271-280, 2001.

[SRESS2008]

B. Sudret. Global sensitivity analysis using polynomial chaos expansions. Reliability Engineering & System Safety, 93(7): 964-979, 2008.

#..
# .. [MT1991]  `M.D. Morris. Factorial Sampling Plans for Preliminary Computational Experiments, Technometrics, 33:2, 161-174, 1991 <https://doi.org/10.1080/00401706.1991.10484804>`_