Sensitivity analysis

This notebook is an element of the risk-engineering.org courseware. It can be distributed under the terms of the Creative Commons Attribution-ShareAlike licence.

Author: Eric Marsden [email protected]


This notebook contains an introduction to use of Python, SciPy, SymPy and the SALib library for sensitivity analysis. Consult the accompanying lecture slides for details of the applications of sensitivity analysis and some intuition and theory of the technique.

In [17]:
import numpy
import scipy.stats
import pandas
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
%matplotlib inline
%config InlineBackend.figure_formats=['svg']

The Rosenbrock function is a classic in uncertainty analysis and sensitivity analysis.

In [18]:
def rosenbrock(x1, x2):
    return 100 * (x2 - x1**2)**2 + (1 - x1)**2

Since we are lucky enough to be working in a small number of dimensions, let's plot the function over the domain $[-2, 2]^2$ to get a feel for its shape.

In [19]:
N = 100
fig = plt.figure()
ax = fig.gca(projection='3d')
x = numpy.linspace(-2, 2, N)
y = numpy.linspace(-2, 2, N)
X, Y = numpy.meshgrid(x, y)
ax.plot_surface(X, Y, rosenbrock(X, Y), alpha=0.8)
ax.set_axis_bgcolor('white')
plt.xlabel(r"$x_1$")
plt.ylabel(r"$x_2$");
/opt/anaconda/anaconda3/lib/python3.5/site-packages/ipykernel/__main__.py:8: MatplotlibDeprecationWarning: The set_axis_bgcolor function was deprecated in version 2.0. Use set_facecolor instead.

Local sensitivity analysis

We can undertake a local sensitivity analysis by using the SymPy package to do symbolic differentiation of the Rosenbrock function, with respect to the two input parameters. You may need to install the SymPy package.

You may be interested in the minireference.com tutorial on SymPy.

In [20]:
import sympy
from sympy.interactive import printing
printing.init_printing(use_latex='mathjax')
In [21]:
x1 = sympy.Symbol('x1')
x2 = sympy.Symbol('x2')
rosen = 100 * (x2 - x1**2)**2 + (1 - x1)**2
d1 = sympy.diff(rosen, x1)
d1
Out[21]:
$$- 400 x_{1} \left(- x_{1}^{2} + x_{2}\right) + 2 x_{1} - 2$$
In [22]:
d2 = sympy.diff(rosen, x2)
d2
Out[22]:
$$- 200 x_{1}^{2} + 200 x_{2}$$

The rosenbrock function looks pretty flat around $(0, 0)$; let's check the local sensitivity in that location. First check $\frac{∂f}{∂x_1}(0, 0)$, then $\frac{∂f}{∂x_2}(0, 0)$

In [23]:
d1.subs({x1: 0, x2: 0})
Out[23]:
$$-2$$
In [24]:
d2.subs({x1: 0, x2: 0})
Out[24]:
$$0$$

The function looks much steeper (higher local sensitivity) around $(-2, -2)$; let's check that numerically.

In [25]:
d1.subs({x1: -2, x2: -2})
Out[25]:
$$-4806$$
In [26]:
d2.subs({x1: -2, x2: -2})
Out[26]:
$$-1200$$

At $(-2, 2)$ the sensitivity should be somewhere in between these two points.

In [27]:
d1.subs({x1: -2, x2: 2})
Out[27]:
$$-1606$$
In [28]:
d2.subs({x1: -2, x2: 2})
Out[28]:
$$-400$$

We can use SciPy's optimization functionality to find the minimum of the Rosenbrock function on the domain $[-2, 2]^2$, then check that (as we expect) the local sensitivity at the minimum is zero. The last argument [2, 2] to the function scipy.optimize.fmin is the starting point of the optimization search.

In [29]:
import scipy
scipy.optimize.fmin(lambda x: rosenbrock(x[0], x[1]), [2, 2])
Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 62
         Function evaluations: 119
Out[29]:
array([ 0.99998292,  0.99996512])

The subs function in SymPy does variable substitution; it allows you to evaluate an expression with given values for the variables (x1 and x2 in this case).

In [30]:
d1.subs({x1: 1, x2: 1})
Out[30]:
$$0$$
In [31]:
d2.subs({x1: 1, x2: 1})
Out[31]:
$$0$$

Global sensitivity analysis

We can use the SALib library (available for download from https://github.com/jdherman/SALib) to undertake a global sensitivity analysis, using the Sobol' method. You may need to install this library, with

pip install SALib

In [32]:
# this will fail if you haven't installed the SALib library properly
from SALib.sample import saltelli
from SALib.analyze import sobol
In [33]:
N = 1000
problem = {
    'num_vars': 2, 
    'names': ['x1', 'x2'], 
    'bounds': [[-2, 2], [-2, 2]]
}
sample = saltelli.sample(problem, N)
Y = numpy.empty([sample.shape[0]])
for i in range(len(Y)):
    x = sample[i]
    Y[i] = rosenbrock(x[0], x[1])
sobol.analyze(problem, Y)
Out[33]:
{'S1': array([ 0.49911118,  0.30018188]),
 'S1_conf': array([ 0.09114489,  0.05701945]),
 'S2': array([[        nan,  0.21626586],
        [        nan,         nan]]),
 'S2_conf': array([[        nan,  0.19439404],
        [        nan,         nan]]),
 'ST': array([ 0.707303  ,  0.51033275]),
 'ST_conf': array([ 0.08995183,  0.08555191])}

S1 contains the first-order sensitivity indices, which tell us how much $x_1$ and $x_2$ each contribute to the overall output variability of the rosenbrock function over the domain $[-2, 2]^2$.

Interpretation: we note that $x_1$ (whose sensitivity index is around 0.5) contributes to roughly half of total output uncertainty, and is roughly two times more influential (or sensitive) over this domain than $x_2$. ST contains the total indices, which include the interaction effects with other variables. The total sensitivity of $x_1$ (around 0.7) indicates that a significant amount (around 20%) of our total output uncertainty is due to the interaction of $x_1$ with other input variables. Since there are only two input variables, we know that this interaction effect must be with $x_2$.

Some sensitivity analysis methods are also able to provide second and third order sensitivity indices. A second order index $s_{i,j}$ tells you the level of interaction effects between $x_i$ and $x_j$. A third order index $s_{i,j,k}$ tells you the level of interaction between three parameters $x_i$, $x_j$ and $x_k$.

Exercise

The Ishigami function is a well known test function for uncertainty analysis and sensitivity analysis (it is highly non-linear). Its definition is given below.

In [34]:
def ishigami(x1, x2, x3):
    return numpy.sin(x1) + 7*numpy.sin(x2)**2 + 0.1 * x3**4 * numpy.sin(x1)

Task: undertake a global sensitivity analysis of the Ishigami function over the domain $[-\pi, \pi]^3$ and estimate the first-order and total sensitivity indices.