{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "# A practical introduction to sensitivity analysis\n", "\n", " \n", "**Leif Rune Hellevik**, Department of Structural Engineering, NTNU\n", "\n", "Date: **Jul 13, 2018**" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# ipython magic\n", "%matplotlib notebook\n", "%load_ext autoreload\n", "%autoreload 2" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# plot configuration\n", "import matplotlib\n", "import matplotlib.pyplot as plt\n", "plt.style.use(\"ggplot\")\n", "# import seaborn as sns # sets another style\n", "matplotlib.rcParams['lines.linewidth'] = 3\n", "fig_width, fig_height = (7.0,5.0)\n", "matplotlib.rcParams['figure.figsize'] = (fig_width, fig_height)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# import modules\n", "import numpy as np\n", "from numpy import linalg as LA\n", "import chaospy as cp\n", "from sensitivity_examples_nonlinear import generate_distributions\n", "from monte_carlo import generate_sample_matrices_mc\n", "from monte_carlo import calculate_sensitivity_indices_mc\n", "import pandas as pd\n", "from operator import index" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Introduction\n", "
\n", "\n", "This practical introduction to sensitivity analysis is based on the\n", "presentation and examples found in [[saltelli_global_2008]](#saltelli_global_2008). To give\n", "the reader an even better hands on experience of the topic, we have\n", "integrated the computations in a python notebook format.\n", "\n", "Many sensitivity analyses reported in the literature are based on\n", "derivatives at set point or point of interest. Indeed such approaches\n", "are based on the fact that the derivative of $\\partial Y_i/\\partial\n", "X_j$ of quantity of interest $Y_i$ as a function of an input variable\n", "$X_j$ can be thought of as the mathematical definition of the\n", "sensitivity of $Y_i$ versus $X_j$.\n", "\n", "However, what is important to keep in mind is that local derivatives\n", "are only informative at the set point in the parameter space at which\n", "they are computed, and do not provide information for the rest of the\n", "parameter space. Naturally, such a linearisation will matter little\n", "for linear models, but for general, nonlinear models, care must be\n", "taken. In particular this is important in situations when the input\n", "parameters are uncertain.\n", "\n", "# Local versus global sensitivity analysis\n", "\n", "Motivation and useful purposes of sensitivity analysis\n", "\n", " * Parameter prioritization of parameters of high sensitivity (importance)\n", "\n", " * Parameter fixation of parameters of low sensitivity (importance)\n", "\n", " * Reveal surprising relations/properties of the model\n", "\n", " * Indentify critical regions in the input parameter space\n", "\n", "## Local approaches based on derivatives\n", "\n", "Many sensitivity analyses found in the scientific literature are based\n", "on derivatives. This fact has naturally a rational basis as the\n", "partial derivative $\\partial y/\\partial Z_i$ of a model predicion $y$\n", "with respect to an input $Z_i$, can be understood as the mathematical\n", "representation of the sensitivity of $y$ with respect to $Z_i$.\n", "\n", "Even though a local, partial derivative approach is computationallyXS\n", "inexpensive it has in general limited usage for nonlinear models. The\n", "derivatives are linearizations of the model sensitivities around the\n", "point in the parameter space at which they are evaluated, and may only\n", "be extrapolated to provide information on the sensitivity in other\n", "regions of the parameter space in the case of a linear model.\n", "\n", "As a first motivation you may consider the ratio areas of circle with\n", "diameter $D$ to the corresponding square with sides $D$, as an optimistic approximation of the portion of the parameter space you might represent with a one factor at the time (OAT) approach:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "\\frac{1}{4} \\; \\frac{\\pi \\, D^2}{D^2} \\approx 3/4 \n", "\\label{_auto1} \\tag{1}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "

Area ratio of circle to a square.

\n", "\n", "\n", "\n", "\n", "\n", "\n", "You may compute the corresonding ratio for a sphere and a cube\n", "yourself. For higher dimensions we provide a code snippet below which\n", "calculates the ratio of a [hypersphere](https://en.wikipedia.org/wiki/N-sphere#Recurrences) to a\n", "hypercube. This is you motivate why you should avoid \"one factor at\n", "the time\" analyses for sensitivity analysis in higher dimensions with multiple parameters." ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# See https://en.wikipedia.org/wiki/N-sphere#Recurrences\n", "\n", "%matplotlib nbagg\n", "import ipywidgets as widgets\n", "import matplotlib.pyplot as plt\n", "import interactive_pie\n", "from interactive_pie import update_pie\n", "\n", "plt.suptitle('Why you should avoid \"one factor at the time\" in higher dimensions.')\n", "w_dim = widgets.IntSlider(min=1, max=20, value=2, description='Dimensions')\n", "widgets.interactive(update_pie, Ndim=w_dim)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Based on the brief motivation above we will present of methods based\n", "on the exploration of the input parameter space by judiciously\n", "selecting samples in that space. Such approaches result in more robust and\n", "informative sensitivity measures, than what would be the result from a local\n", "derivative app\t\troach at the center of the parameter space.\n", "\n", "To introduce the methods of sensitivity analysis, we shall\n", "start from derivatives and illustrate them on a very simple linear\n", "model.\n", "\n", "\n", "# A simple linear model\n", "\n", "As an simple linear model example consider:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "Y = \\sum_{i=1}^{r} \\Omega_i \\, Z_i\n", "\\label{eq:linear_model} \\tag{2}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "where the input factors are $\\mathbf{X} = (\\Omega_1, \\Omega_2, \\ldots,\n", "\\Omega_r, Z_1, Z_2, \\ldots, Z_r)$. For simplicity we assume that the\n", "model output $Y$ of ([2](#eq:linear_model)) is a single variable and\n", "that the $\\Omega s$ are fixed coefficients or weights." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", " \\Omega_1=\\Omega_2=\\ldots=\\text{constant}\n", "\\label{_auto2} \\tag{3}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Consequently, the true factors of ([2](#eq:linear_model)) are just\n", "$(Z_1, Z_2, \\ldots, Z_r)$. The individual variables\n", "$Z_i$ are taken to normally distributed with mean zero" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} Z_i \\sim N(0, \\sigma_{Z_i}), \\qquad i=1,2, \\ldots, r\n", "\\label{eq:NZi} \\tag{4} \\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As the predicted value $Y$ of the model in ([2](#eq:linear_model)) is\n", "linear combination of normally distributed factors, it is easy to\n", "verify (see exercices in [[saltelli_global_2008]](#saltelli_global_2008)) that $Y$ also will\n", "be normally distributed with:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "\\bar{y} = \\sum_{i=1}^{r} \\Omega_i \\; \\bar{z}_i, \\qquad \\sigma_Y = \\sqrt{\\sum_{i=1}^{r} \\Omega_i^2 \\, \\sigma_{Z_i}^2}\n", "\\label{eq:analytic_mean_std} \\tag{5}\t\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Furthermore, we order the factors from the most certain to the less\n", "certain, i.e.:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", " \\sigma_{Z_1} < \\sigma_{Z_2} < \\ldots < \\sigma_{Z_r}\n", "\\label{_auto3} \\tag{6}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Scatterplots versus derivatives\n", "\n", "We have implemented the simple linear model in ([2](#eq:linear_model)) in python as:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# start the linear model\n", "def linear_model(w, z):\n", " return np.sum(w*z, axis=1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To hold the mean and the standard deviation of all the input factors\n", "we use a numpy-array of size $r\\times 2$, with one row per factor,\n", "where the first column holds the mean whereas the second column holds\n", "the standard deviation. The weights $\\Omega_{1\\ldots r}$ are stored in a numpy-vector." ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [], "source": [ " # Set mean (column 0) and standard deviations (column 1) for each factor z. Nrv=nr. rows\n", " Nrv = 4 # number of random variables \n", " zm = np.array([[0., i] for i in range(1, Nrv + 1)])\n", "\n", " # The above \"list comprehension\" is equivalent to the next for lines\n", " # zm = np.zeros((Nrv, 2))\n", " # zm[0, 1] = 1\n", " # zm[1, 1] = 2\n", " # zm[2, 1] = 3\n", " # zm[3, 1] = 4\n", "\n", " # Set the weight\n", " c = 2\n", " w = np.ones(Nrv) * c" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We may now perform a Monte Carlo experiment on our model by generating $N$ samples from the distributions of each factor and an input sample is thus produced:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "\\mathbf{Z} = \\left [\n", "\\begin{array}{cccc}\n", "Z_{1,1} & Z_{1,2} & \\ldots & Z_{1,N} \\\\ \n", "Z_{2,1} & Z_{2,2} & \\ldots & Z_{2,N}\\\\ \n", "\\vdots & \\vdots & \\vdots & \\vdots \\\\ \n", "Z_{r,1} & Z_{r,2} & \\ldots & Z_{r,N}\n", "\\end{array} \n", "\\right ]\n", "\\label{eq:mc_sample} \\tag{7}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We may the compute a value of $Y$ from ([2](#eq:linear_model)) for each\n", "column in ([7](#eq:mc_sample)) to produce a solution vector\n", "$\\mathbf{Y}$. Having sampled $N$ values from each input factor we may\n", "produce $r$ scatter plots, by projecting in turn the $N$ values of\n", "$\\mathbf{Y}$ against the $N$ values of each of the $r$ input factors." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "\\mathbf{Y} = \\left [\n", "\\begin{array}{c}\n", "y_1 \\\\ \n", "y_2 \\\\ \n", "\\vdots \\\\ \n", "y_N\n", "\\end{array}\n", "\\right ]\n", "\\label{eq:mc_solution} \\tag{8}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [], "source": [ " # Generate distributions for each element in z and sample\n", " Ns = 500\n", " # jpdf = generate_distributions(zm)\n", " \n", " pdfs = []\n", "\n", " for i, z in enumerate(zm):\n", " pdfs.append(cp.Normal(z[0], z[1]))\n", "\n", " jpdf = cp.J(*pdfs)\n", "\n", " # generate Z\n", " Z = jpdf.sample(Ns)\n", " # evaluate the model\n", " Y = linear_model(w, Z.transpose())\n", " print(np.var(Y))\n", "\n", " # Scatter plots of data for visual inspection of sensitivity\n", " fig=plt.figure()\n", " for k in range(Nrv):\n", " plt.subplot(2, 2, k + 1)\n", " plt.plot(Z[k, :], Y[:], '.')\n", " xlbl = 'Z' + str(k)\n", " plt.xlabel(xlbl)\n", " \n", " fig.tight_layout() # adjust subplot(s) to the figure area." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that the assumption of independent factors $Z_i$ allows us to sample\n", "each $Z_i$ independently from its own marginal distribution. We store\n", "all the samples for all the factors $Z_i$ in the the numpy array\n", "`Z[i,:]`, where $i$ corresponds to $Z_i$ as:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " pdf.append(cp.Normal(z[0],z[1]))\n", " Z[i,:]=pdf[i].sample(N)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "From the scatterplots generated by the python code above we\n", "intuitively get the impression that $Y$ is more sensitive to $Z_4$\n", "than to $Z_3$, and that $Y$ is more sensitive to $Z_3$ than to $Z_3$,\n", "and that we may order the factors my influence on $Y$ as:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "Z_4 > Z_3 > Z_2 > Z_1 \n", "\\label{eq:scatter_plot_rank} \\tag{9}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Our intuitive notion of influence is based on that there is more shape\n", "(or better pattern) in the plot for $Z_4$ than for $Z_3$ and likewise.\n", "\n", "For our simple linear model in ([2](#eq:linear_model)) we are in the\n", "fortunate situation that we may compute the local derivatives analyticaly:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "S_{Z_i}^{p} = \\frac{\\partial Y}{\\partial Z_i} = \\Omega_i\n", "\\label{eq:Sp} \\tag{10}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In our code example we set all the $\\Omega_i=2$ for $i=1,\\ldots,4$,\n", "and according to the local sensitivity meansure $S_{Z_i}^{p}$ in\n", "([10](#eq:Sp)) all the input factors $Z_i$s are equally important and\n", "independent of the variation of each factor. This measure is clearly\n", "at odds with the ranking of influence based on the scatterplots in\n", "([9](#eq:scatter_plot_rank)) and is an indication of the usefullness of\n", "scatterplots in sensititivy analysis. However, the bidimensional\n", "scatterplots may in some cases be deceiving and lead to type II\n", "errors (i.e. failure to identify influential parameters). ref to Saltelli 2004...\n", "\n", "Most sensitivity measures aim to preserve the rich information\n", "provided by the scatterplots in a condensed format. The challenge is\n", "how to rank the factors rapidly and automatically without having to\n", "inspect many scatterplots in situations with many input\n", "factors. Another challenge with scatterplots is that sensitivities for\n", "sets cannot be visualized, while luckily compact sensitivity measures may be\n", "defined in such cases.\n", "\n", "# Normalized derivatives\n", "\n", "A simple way to improve the derivative sensitivity measure $S_{Z_i}^{p}$ in\n", "([10](#eq:Sp)) is to scale the input-output variables with their standard deviations:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "S_{Z_i}^{\\sigma} = \\frac{\\partial Y/\\sigma_Y}{\\partial Z_i/\\sigma_{Z_i}} = \\frac{\\sigma_{Z_i}}{\\sigma_{Y}} \\; \\frac{\\partial Y}{\\partial Z_i}\n", "\\label{eq:Ss} \\tag{11}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In case of our simple linear model ([2](#eq:linear_model)) we get from\n", "([11](#eq:Ss)):" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "\\left (S_{Z_i}^{\\sigma} \\right)^2 = \\left( \\frac{\\sigma_{Z_i}}{\\sigma_{Y}}\\right)^2 \\; \\left (\\frac{\\partial Y}{\\partial Z_i}\\right)^2 = \\left( \\frac{\\sigma_{Z_i}\\, \\Omega_i}{\\sigma_{Y}}\\right)^2 \\; \\qquad \\textsf{which may be rearranged to:} \\qquad \\sigma_y^2 \\, (S_{Z_i}^{\\sigma})^2 = \\left ( \\Omega_{i} \\sigma_{Y} \\right )^2\n", "\\label{eq:Ss_simple} \\tag{12}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Based on the linearity of our model we previously found ([5](#eq:analytic_mean_std)) which also yields:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", " \\sigma_Y^2 = \\sum_{i=1}^{r} \\left(\\Omega_i^2 \\, \\sigma_{Z_i}\\right)^2\n", "\\label{eq:Ss_model_ded} \\tag{13}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As both ([13](#eq:Ss_model_ded)) and ([12](#eq:Ss_simple)) must hold simultaneously we get" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "\\left (S_{Z_i}^{\\sigma} \\right)^2=1 \n", "\\label{eq:Ss1} \\tag{14}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The normalized derivative measure of sensitivity in ([11](#eq:Ss)) is\n", "more convincing than ([10](#eq:Sp)): first, as it involves both the\n", "weights $\\Omega_i$ and the factors $Z_i$ in ([2](#eq:linear_model));\n", "second as the measures are properly scaled and summarizes to one,\n", "which allows for an easy interpretation of the output sensitivity with\n", "respect to each of the input factors." ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [], "source": [ " # Theoretical sensitivity indices\n", " std_y = np.sqrt(np.sum((w * zm[:, 1])**2))\n", " s = w * zm[:,1]/std_y\n", " \n", " print(\"\\nTheoretical sensitivity indices\\n\")\n", " row_labels= ['S_'+str(idx) for idx in range(1,Nrv+1)]\n", " print(pd.DataFrame(s**2, columns=['S analytic'],index=row_labels).round(3))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Based on samples of the random input variables and \n", "subsequent model evaluations, we may estimate the standard deviation\n", "of $\\mathbf{Y}$ and compute the relative error with respect to the\n", "theoretical value. You may change the number of sample above,\n", "i.e. $N$, and see how $N$ influence the estimates." ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [], "source": [ " # Expectation and variance from sampled values\n", " \n", " print(\"Expectation and std from sampled values\\n\")\n", " print('std(Y)={:2.3f} and relative error={:2.3f}'.format(np.std(Y, 0), (np.std(Y, 0) - std_y) / std_y))\n", " print('mean(Y)={:2.3f} and E(Y)={:2.3}'.format(np.mean(Y, 0), np.sum(zm[:,0]*w)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that `Ns` is the size of our Monte Carlo experiment, corresponding\n", "to the number of times we have evaluated our simple linear model\n", "([2](#eq:linear_model)). The evaluation of the model is normally the\n", "most computationally expensive part of the analysis, and for that\n", "reasons `Ns` is referred to as the `cost` of the analysis.\n", "\n", "# Conditional variances\n", "\n", "As noted previously, the importance of a factor $Z_i$ is manifested\n", "the existence of a `shape` or `pattern` in the model outputs\n", "$Y$. Conversely, a uniform cloud of output points $Y$ as a function of\n", "$Z_i$ is a symptom, albeit not a proof, indicating that $Z_i$ is a\n", "noninfluential factor. In this section we seek to demonstrate that\n", "conditional variances is a usefull means to quantify the `shape` or\n", "`pattern` in the outputs.\n", "\n", "The shape in the outputs $Y$ for a given $Z_i$, may be seen in the\n", "scatterplot as of $Y$ versus $Z_i$. In particular, we may cut the\n", "$Z_i$-axis into slices and assess how the distribution of the outputs\n", "$Y$ changes from slice to slice. This is illustrated in the code\n", "snippet below, where the slices are identified with vertical dashed\n", "lines at equidistant locations on each $Z_i$-axis, $i=1, \\ldots,4$." ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false }, "outputs": [], "source": [ " # # Scatter plots of data, z-slices, and linear model\n", " fig=plt.figure()\n", "\n", " Ndz = 10 # Number of slices of the Z-axes\n", "\n", " Zslice = np.zeros((Nrv, Ndz)) # array for mean-values in the slices\n", " ZBndry = np.zeros((Nrv, Ndz + 1)) # array for boundaries of the slices\n", " dz = np.zeros(Nrv)\n", "\n", " for k in range(Nrv):\n", " plt.subplot(2, 2, k + 1)\n", "\n", " zmin = np.min(Z[k, :])\n", " zmax = np.max(Z[k, :]) # each Z[k,:] may have different extremas\n", " dz[k] = (zmax - zmin) / Ndz\n", "\n", " ZBndry[k, :] = np.linspace(zmin, zmax, Ndz + 1) # slice Zk into Ndz slices\n", " Zslice[k, :] = np.linspace(zmin + dz[k] / 2., zmax - dz[k] / 2., Ndz) # Midpoint in the slice\n", "\n", " # Plot the the vertical slices with axvline\n", " for i in range(Ndz):\n", " plt.axvline(ZBndry[k, i], np.amin(Y), np.amax(Y), linestyle='--', color='.75')\n", "\n", " # Plot the data\n", " plt.plot(Z[k, :], Y[:], '.')\n", " xlbl = 'Z' + str(k)\n", " plt.xlabel(xlbl)\n", " plt.ylabel('Y')\n", "\n", " Ymodel = w[k] * Zslice[k, :] # Produce the straight line\n", "\n", " plt.plot(Zslice[k, :], Ymodel)\n", "\n", " ymin = np.amin(Y); ymax = np.amax(Y)\n", " plt.ylim([ymin, ymax])\n", " \n", " fig.tight_layout() # adjust subplot(s) to the figure area." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note, that average value of $Y$ in a very thin slice, corresponds to\n", "keeping $Z_i$ fixed while averaging over all output values of $Y$ due\n", "to all-but $Z_i$, which corresponds to the conditional expected value:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "E_{Z_{\\sim i}} (Y\\;|\\;Z_i) \n", "\\label{_auto4} \\tag{15}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For convenience we let $Z_{\\sim i}$ denote `all-but` $Z_i$. Naturally,\n", "a measure of how much $E_{Z_{\\sim i}} (Y\\;|\\;Z_i)$ varies in the range\n", "of $Z_i$ is given by the conditional variance:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "\\text{V}_{Z_i}(E_{Z_{\\sim i}} (Y\\;|\\;Z_i))\n", "\\label{_auto5} \\tag{16}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Further, the variance the output $Y$ may be decomposed into:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "\\text{V}(Y) = E_{Z_i} ( V_{Z_{\\sim i}} (Y \\; | Z_{i})) + \\text{V}_{Z_i}(E_{Z_{\\sim i}} (Y\\;|\\;Z_i))\n", "\\label{eq:VarDecomp} \\tag{17}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A large $\\text{V}_{Z_i}(E_{Z_{\\sim i}} (Y\\;|\\;Z_i))$ will imply that\n", "$Z_i$ is an important factor and is therefore coined the first-order\n", "effect of $Z_i$ on $Y$, and its fraction of the total variation of $Y$ is expressed by $S_i$, `the first-order sensitivity index` of $Z_i$ on $Y$:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "S_i = \\frac{\\text{V}_{Z_i}(E_{Z_{\\sim i}} (Y\\;|\\;Z_i))}{\\text{V}(Y)}\n", "\\label{_auto6} \\tag{18}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "By ([17](#eq:VarDecomp)), $S_i$ is number always in the range $[0,1]$,\n", "and a high value implies an important factor." ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false }, "outputs": [], "source": [ " # # Scatter plots of averaged y-values per slice, with averaged data\n", "\n", " Zsorted = np.zeros_like(Z)\n", " Ysorted = np.zeros_like(Z)\n", " YsliceMean = np.zeros((Nrv, Ndz))\n", "\n", " fig=plt.figure()\n", " for k in range(Nrv):\n", " plt.subplot(2, 2, k + 1)\n", "\n", " # sort values for Zk, \n", " sidx = np.argsort(Z[k, :]) #sidx holds the indexes for the sorted values of Zk\n", " Zsorted[k, :] = Z[k, sidx].copy()\n", " Ysorted[k, :] = Y[sidx].copy() # Ysorted is Y for the sorted Zk\n", "\n", " for i in range(Ndz):\n", " plt.axvline(ZBndry[k, i], np.amin(Y), np.amax(Y), linestyle='--', color='.75')\n", "\n", " # find indexes of z-values in the current slice\n", " zidx_range = np.logical_and(Zsorted[k, :] >= ZBndry[k, i], Zsorted[k, :] < ZBndry[k, i + 1])\n", "\n", " if np.any(zidx_range): # check if range has elements\n", " YsliceMean[k, i] = np.mean(Ysorted[k, zidx_range])\n", " else: # set value to None if noe elements in z-slice\n", " YsliceMean[k, i] = None\n", "\n", " plt.plot(Zslice[k, :], YsliceMean[k, :], '.')\n", " \n", " \n", "\n", " # # Plot linear model\n", " Nmodel = 3\n", " zmin = np.min(Zslice[k, :])\n", " zmax = np.max(Zslice[k, :])\n", "\n", " zvals = np.linspace(zmin, zmax, Nmodel)\n", " #linear_model\n", " Ymodel = w[k] * zvals\n", " plt.plot(zvals, Ymodel)\n", "\n", " xlbl = 'Z' + str(k)\n", " plt.xlabel(xlbl)\n", "\n", " plt.ylim(ymin, ymax)\n", " \n", " fig.tight_layout() # adjust subplot(s) to the figure area.\n", " \n", " SpoorMan=[np.nanvar(YsliceMean[k,:],axis=0)/np.var(Y) for k in range(4)] \n", " print(SpoorMan)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# How to compute the sensitivity indices\n", "\n", "Below we will demostrate how the Sobol sensitivity indices may be\n", "computed with two approaches; the Monte Carlo method and the\n", "polynomial chaos expansion method.\n", "\n", "### Monte Carlo\n", "\n", "Below some code snippets are provided to illustrate how we may compute\n", "the Soboil indices with the MCM. For the interested reader we have also\n", "writen a seperate and more detailed notebook [A brief introduction to\n", "UQ and SA with the Monte Carlo method](monte_carlo.ipynb)." ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# calculate sens indices of non additive model\n", "def mc_sensitivity_linear(Ns, jpdf, w, sample_method='R'):\n", "\n", " Nrv = len(jpdf)\n", "\n", " # 1. Generate sample matrices\n", " A, B, C = generate_sample_matrices_mc(Ns, Nrv, jpdf, sample_method)\n", "\n", " # 2. Evaluate the model\n", " Y_A, Y_B, Y_C = evaluate_linear_model(A, B, C, w)\n", "\n", " # 3. Approximate the sensitivity indices\n", " S, ST = calculate_sensitivity_indices_mc(Y_A, Y_B, Y_C)\n", "\n", " return A, B, C, Y_A, Y_B, Y_C, S, ST\n", "# end calculate sens indices of non additive model\n", "\n", "\n", "# model evaluation\n", "def evaluate_linear_model(A, B, C, w):\n", "\n", " number_of_parameters = A.shape[1]\n", " number_of_sampless = A.shape[0]\n", " # 1. evaluate sample matrices A\n", " Y_A = linear_model(w, A)\n", "\n", " # 2. evaluate sample matrices B\n", " Y_B = linear_model(w, B)\n", "\n", " # 3. evaluate sample matrices C\n", " Y_C = np.empty((number_of_sampless, number_of_parameters))\n", " for i in range(number_of_parameters):\n", " z = C[i, :, :]\n", " Y_C[:, i] = linear_model(w, z)\n", "\n", " return Y_A, Y_B, Y_C" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false }, "outputs": [], "source": [ " # Monte Carlo\n", " # get joint distributions\n", " jpdf = generate_distributions(zm)\n", "\n", " Ns_mc = 1000000\n", " # calculate sensitivity indices\n", " A_s, B_s, C_s, f_A, f_B, f_C, S_mc, ST_mc = mc_sensitivity_linear(Ns_mc, jpdf, w)\n", "\n", " Sensitivities=np.column_stack((S_mc,s**2))\n", " row_labels= ['S_'+str(idx) for idx in range(1,Nrv+1)]\n", " print(\"First Order Indices\")\n", " print(pd.DataFrame(Sensitivities,columns=['Smc','Sa'],index=row_labels).round(3))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Polynomial chaos expansion\n", "\n", "As for the MCM some code snippets are provided to illustrate how we may compute\n", "the Soboil indices with the polynomial chaos expansions using `chaospy`. A more in dept treatment of `chaospy` and its usage is provided in the separate notebook [A practical introduction to polynomial chaos with the chaospy package](introduction_gpc.ipynb)." ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false }, "outputs": [], "source": [ " # Polychaos computations\n", " Ns_pc = 80\n", " samples_pc = jpdf.sample(Ns_pc)\n", " polynomial_order = 4\n", " poly = cp.orth_ttr(polynomial_order, jpdf)\n", " Y_pc = linear_model(w, samples_pc.T)\n", " approx = cp.fit_regression(poly, samples_pc, Y_pc, rule=\"T\")\n", "\n", " exp_pc = cp.E(approx, jpdf)\n", " std_pc = cp.Std(approx, jpdf)\n", " print(\"Statistics polynomial chaos\\n\")\n", " print('\\n E(Y) | std(Y) \\n')\n", " print('pc : {:2.5f} | {:2.5f}'.format(float(exp_pc), std_pc))\n", " \n", " \n", " S_pc = cp.Sens_m(approx, jpdf)\n", "\n", " Sensitivities=np.column_stack((S_mc,S_pc, s**2))\n", " print(\"\\nFirst Order Indices\")\n", " print(pd.DataFrame(Sensitivities,columns=['Smc','Spc','Sa'],index=row_labels).round(3))\n", "\n", "# print(\"\\nRelative errors\")\n", "# rel_errors=np.column_stack(((S_mc - s**2)/s**2,(S_pc - s**2)/s**2))\n", "# print(pd.DataFrame(rel_errors,columns=['Error Smc','Error Spc'],index=row_labels).round(3))" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false }, "outputs": [], "source": [ " # Polychaos convergence\n", " Npc_list = np.logspace(1, 3, 10).astype(int)\n", " error = []\n", "\n", " for i, Npc in enumerate(Npc_list):\n", " Zpc = jpdf.sample(Npc)\n", " Ypc = linear_model(w, Zpc.T)\n", " Npol = 4\n", " poly = cp.orth_chol(Npol, jpdf)\n", " approx = cp.fit_regression(poly, Zpc, Ypc, rule=\"T\")\n", " s_pc = cp.Sens_m(approx, jpdf)\n", " error.append(LA.norm((s_pc - s**2)/s**2))\n", "\n", " plt.figure()\n", " plt.semilogy(Npc_list, error)\n", " _=plt.xlabel('Nr Z')\n", " _=plt.ylabel('L2-norm of error in Sobol indices')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# References\n", "\n", " 1.
**A. Saltelli**. \n", " *Global Sensitivity Analysis : the Primer*,\n", " John Wiley,\n", " 2008." ] } ], "metadata": {}, "nbformat": 4, "nbformat_minor": 2 }