Visual explanation of first order Sobol' indices

sobol_visually.py

r"""Visual explanation of Sobol' indices.

Sobol' indices are metrics to express sensitivity of the output from
perturbations comming from input parameters. First order indices write

.. math:: S_{x_i} = \frac{\mathbb{V}_i(Y)}{\mathbb{V}[Y]} =
          \frac{\mathbb{\mathbb{V}}[\mathbb{E}(Y|x_i)]}{\mathbb{V}[Y]}

The following is using the Ishigami function

.. math:: F = \sin(x_1)+7\sin(x_2)^2+0.1x_3^4\sin(x_1), x\in [-\pi, \pi]^3.

This function exhibits strong nonlinearity and nonmonotonicity. It is
particularly interesting because the first order indice of :math:`S_{x_3} = 0`
whereas it's total order is :math:`S_{T_{x_3}} = 0.244`. Note that on second
order indices, :math:`S_{x_1,x_3} = 0.244`. It means that almost 25% of the
observed variance on the output is due to the correlations between :math:`x_3`
and :math:`x_1`, although :math:`x_3` by itself has no impact on the output.

Starting from a scatter plots of the output with respect to each parameter.
By conditioning the output value by given values of the parameter,
the conditional output mean is computed (dots). It corresponds to the term
:math:`\mathbb{E}(Y|x_i)]`. Taking the variance of this term gives the
numerator of the Sobol' indices. Looking at :math:`x_3`, the variance of the
mean is close to zero leading to :math:`S_{x_3} = 0`. But we can further
observe that the variance of the output is not constant along the parameter
values of :math:`x_3`. This heteroscedasticity is explained by higher order
interactions. Moreover, an heteroscedasticity is also noticeable on
:math:`x_1` leading to an interaction between :math:`x_3` and :math:`x_1`.
On :math:`x_2`, the variance seems to be constant and thus null interaction
with this parameter can be supposed. This case is fairly simple to analyse
visually---although it is only a qualitative analysis. Nevertheless, when the
number of input parameters increases such analysis becomes unrealistic as it
would be difficult to conclude on high order terms.

---------------------------

MIT License

Copyright (c) 2018 Pamphile Tupui ROY

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
"""
import numpy as np
from batman.functions import Ishigami
from batman.space import Space
import matplotlib.pyplot as plt

p_labels = ['$x_1$', '$x_2$', '$x_3$']
sample = Space([[-np.pi, -np.pi, -np.pi], [np.pi, np.pi, np.pi]])
sample.sampling(1000)

n_dim = sample.shape[1]

func = Ishigami()
output = func(sample).flatten()

var_y = np.var(output)
mini = np.min(output)
maxi = np.max(output)
n_bins = 10
bins = np.linspace(-np.pi, np.pi, num=n_bins, endpoint=False)
dx = bins[1] - bins[0]

# Sobol explanation
fig, ax = plt.subplots(1, n_dim)
for i in range(n_dim):
    xi = sample[:, i]

    ax[i].scatter(xi, output, marker='+')
    ax[i].set_xlabel(p_labels[i])

    for bin_ in bins:

        idx = np.where((bin_ <= xi) & (xi <= bin_ + dx))
        xi_ = xi[idx]
        y_ = output[idx]

        ave_y_ = np.mean(y_)

        ax[i].plot([bin_ + dx / 2] * 2, [mini, maxi], c='k')
        ax[i].scatter(bin_ + dx / 2, ave_y_, c='r')

ax[0].set_ylabel('Y')

plt.tight_layout()
plt.show()