sensitivity_introduction_interactive.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- dom:TITLE: A practical introduction to sensitivity analysis -->\n",
    "\n",
    "# A practical introduction to sensitivity analysis\n",
    "\n",
    "<!-- dom:AUTHOR: Leif Rune\n",
    "Hellevik at Department of Structural Engineering, NTNU -->\n",
    "<!-- Author: -->\n",
    "**Leif Rune Hellevik**, Department of Structural Engineering, NTNU\n",
    "\n",
    "Date: **Sep\n",
    "29, 2018**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "1"
    }
   },
   "outputs": [],
   "source": [
    "# import modules\n",
    "import ipywidgets as widgets\n",
    "from IPython.display import display\n",
    "import numpy as np\n",
    "import chaospy as cp\n",
    "import os, sys, inspect\n",
    "# Use this if you want to include modules from a subfolder\n",
    "cmd_subfolder = os.path.realpath(os.path.abspath(os.path.join(os.path.split(inspect.getfile( inspect.currentframe() ))[0],\"python_source\")))\n",
    "if cmd_subfolder not in sys.path:\n",
    "     sys.path.insert(0, cmd_subfolder)\n",
    "\n",
    "from monte_carlo import generate_sample_matrices_mc, calculate_sensitivity_indices_mc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction <div id=\"sec:introduction\"></div>\n",
    "\n",
    "This practical introduction to\n",
    "sensitivity analysis is based on the\n",
    "presentation and examples available from\n",
    "[saltelli_global_2008](#saltelli_global_2008). To give\n",
    "the reader an even\n",
    "better hands on experience of the topic, we have\n",
    "integrated the computations in\n",
    "a python notebook format.\n",
    "\n",
    "Many sensitivity analyses reported in the literature\n",
    "are based on\n",
    "derivatives at set point or point of interest. Indeed such\n",
    "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\n",
    "be thought of as the mathematical definition of the\n",
    "sensitivity of $Y_i$ versus\n",
    "$X_j$.\n",
    "\n",
    "However, what is important to keep in mind is that local derivatives\n",
    "are\n",
    "only informative at the set point in the parameter space at which\n",
    "they are\n",
    "computed, and do not provide information for the rest of the\n",
    "parameter space.\n",
    "Naturally, such a linearisation will matter little\n",
    "for linear models, but care must be\n",
    "taken for general, nonlinear models.  This is especially important\n",
    "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",
    "* Identify critical regions in the input parameter space\n",
    "\n",
    "## Local approaches based on derivatives\n",
    "\n",
    "Many\n",
    "sensitivity analyses found in the scientific literature are based\n",
    "on\n",
    "derivatives.  This fact has naturally a rational basis as the\n",
    "partial derivative\n",
    "$\\partial y/\\partial Z_i$ of a model prediction $y$\n",
    "with respect to an input\n",
    "$Z_i$, can be understood as the mathematical\n",
    "representation of the sensitivity\n",
    "of $y$ with respect to $Z_i$.\n",
    "\n",
    "Although a local, partial derivative approach\n",
    "is computationally\n",
    "inexpensive, it has in general limited usage for nonlinear\n",
    "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\n",
    "extrapolated to provide information on the sensitivity in other\n",
    "regions of the\n",
    "parameter space in the case of a linear model.\n",
    "\n",
    "To illustrate the fraction of\n",
    "the parameter space one may\n",
    "explore with the local partial\n",
    "derivative approach (also called the\n",
    "one factor at the time (OAT) approach), we\n",
    "provide a code snippet which\n",
    "calculates the ratio of a\n",
    "[hypersphere](https://en.wikipedia.org/wiki/N-sphere#Recurrences) to a\n",
    "hypercube."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "2"
    }
   },
   "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": [
    "<!-- dom:FIGURE: [figs/hypersphere.png, width=400 frac=0.7] Ratio of hypersphere\n",
    "volume to hypercube volume.  -->\n",
    "<!-- begin figure -->\n",
    "\n",
    "<p>Ratio of hypersphere\n",
    "volume to hypercube volume.</p>\n",
    "<img src=\"figs/hypersphere.png\" width=400>\n",
    "\n",
    "<!--\n",
    "end figure -->\n",
    "\n",
    "Based on the brief motivation above, we will present of\n",
    "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 a local\n",
    "derivative approach at the center of the parameter space.\n",
    "\n",
    "To introduce the\n",
    "methods of sensitivity analysis, we shall\n",
    "start from derivatives and illustrate\n",
    "them on a very simple linear\n",
    "model.\n",
    "\n",
    "# A simple linear model\n",
    "\n",
    "As an simple\n",
    "linear model example consider:\n",
    "\n",
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"eq:linear_model\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "Y = \\sum_{i=1}^{r} \\Omega_i \\, Z_i\n",
    "\\label{eq:linear_model}\n",
    "\\tag{1}\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "where the input factors are $\\mathbf{X} = (\\Omega_1, \\Omega_2, \\ldots,\n",
    "\\Omega_r,\n",
    "Z_1, Z_2, \\ldots, Z_r)$. For simplicity we assume that the\n",
    "model output $Y$ of\n",
    "([1](#eq:linear_model)) is a scalar and\n",
    "that the $\\Omega s$ are fixed\n",
    "coefficients or weights.\n",
    "\n",
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto1\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    " \\Omega_1=\\Omega_2=\\ldots=\\text{constant}\n",
    "\\label{_auto1}\n",
    "\\tag{2}\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "Consequently, the true factors of ([1](#eq:linear_model)) are just\n",
    "$(Z_1, Z_2,\n",
    "\\ldots, Z_r)$. The individual variables\n",
    "$Z_i$ are taken as normally distributed\n",
    "with mean zero\n",
    "\n",
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"eq:NZi\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation} Z_i \\sim N(0, \\sigma_{Z_i}), \\qquad i=1,2, \\ldots, r\n",
    "\\label{eq:NZi} \\tag{3} \\end{equation}\n",
    "$$\n",
    "\n",
    "As the predicted value $Y$ of the model in ([1](#eq:linear_model)) is\n",
    "linear\n",
    "combination of normally distributed factors, it is easy to\n",
    "verify (see exercices\n",
    "in [saltelli_global_2008](https://onlinelibrary.wiley.com/doi/book/10.1002/9780470725184) that $Y$ also will\n",
    "be\n",
    "normally distributed with:\n",
    "\n",
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"eq:analytic_mean_std\"></div>\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\bar{y} = \\sum_{i=1}^{r} \\Omega_i \\; \\bar{z}_i, \\qquad\n",
    "\\sigma_Y = \\sqrt{\\sum_{i=1}^{r} \\Omega_i^2 \\, \\sigma_{Z_i}^2}\n",
    "\\label{eq:analytic_mean_std} \\tag{4}\t\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "Furthermore, we order the factors from the most certain to the less\n",
    "certain,\n",
    "i.e.:\n",
    "\n",
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto2\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    " \\sigma_{Z_1} <  \\sigma_{Z_2} <  \\ldots  <  \\sigma_{Z_r}\n",
    "\\label{_auto2} \\tag{5}\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "# Scatterplots versus derivatives\n",
    "\n",
    "We have implemented the simple linear model\n",
    "in ([1](#eq:linear_model)) in python as:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "3"
    }
   },
   "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\n",
    "numpy-array of size $r\\times 2$, with one row per factor,\n",
    "where the first column\n",
    "holds the mean whereas the second column holds\n",
    "the standard deviation. The\n",
    "weights $\\Omega_{1\\ldots r}$ are stored in a numpy-vector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "4"
    }
   },
   "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 weights\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$\n",
    "samples from the distributions of each factor and an input sample is thus\n",
    "produced:\n",
    "\n",
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"eq:mc_sample\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\mathbf{Z} = \\left [\n",
    "\\begin{array}{cccc}\n",
    "Z_{1,1} & Z_{1,2}  &\n",
    "\\ldots & Z_{1,N} \\\\ \n",
    "Z_{2,1} & Z_{2,2}  & \\ldots & Z_{2,N}\\\\ \n",
    "\\vdots & \\vdots &\n",
    "\\vdots & \\vdots \\\\ \n",
    "Z_{r,1} & Z_{r,2}  & \\ldots & Z_{r,N}\n",
    "\\end{array} \n",
    "\\right ]\n",
    "\\label{eq:mc_sample} \\tag{6}\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "We may the compute a value of $Y$ from ([1](#eq:linear_model)) for each\n",
    "column\n",
    "in ([6](#eq:mc_sample)) to produce a solution vector\n",
    "$\\mathbf{Y}$. Having\n",
    "sampled $N$ values from each input factor we may\n",
    "produce $r$ scatter plots, by\n",
    "projecting in turn the $N$ values of\n",
    "$\\mathbf{Y}$ against the $N$ values of each\n",
    "of the $r$ input factors.\n",
    "\n",
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"eq:mc_solution\"></div>\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{7}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "5"
    }
   },
   "outputs": [],
   "source": [
    "# Generate distributions for each element in z and sample\n",
    "Ns = 50\n",
    "\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())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the assumption of independent factors $Z_i$ allows us to sample\n",
    "each\n",
    "$Z_i$ independently from its own marginal distribution. We store\n",
    "all the samples\n",
    "for all the factors $Z_i$ in the the numpy array\n",
    "`Z[i,:]`, where $i$ corresponds\n",
    "to $Z_i$ as:\n",
    "\n",
    "pdf.append(cp.Normal(z[0],z[1]))\n",
    "Z[i,:]=pdf[i].sample(N)\n",
    "\n",
    "From the scatterplots generated by the python code above we\n",
    "intuitively get the\n",
    "impression that $Y$ is more sensitive to $Z_4$\n",
    "than to $Z_3$, and that $Y$ is\n",
    "more sensitive to $Z_3$ than to $Z_3$,\n",
    "and that we may order the factors my\n",
    "influence on $Y$ as:\n",
    "\n",
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"eq:scatter_plot_rank\"></div>\n",
    "$$\n",
    "\\begin{equation}\n",
    "Z_4 > Z_3 > Z_2 > Z_1 \n",
    "\\label{eq:scatter_plot_rank} \\tag{8}\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "Our intuitive notion of influence is based on that there is more shape\n",
    "(or\n",
    "better pattern) in the plot for $Z_4$ than for $Z_3$ and likewise.\n",
    "\n",
    "For our\n",
    "simple linear model in ([1](#eq:linear_model)) we are in the\n",
    "fortunate situation\n",
    "that we may compute the local derivatives analyticaly:\n",
    "\n",
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"eq:Sp\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "S_{Z_i}^{p} = \\frac{\\partial Y}{\\partial Z_i} = \\Omega_i\n",
    "\\label{eq:Sp} \\tag{9}\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "In our code example we set all the $\\Omega_i=2$ for $i=1,\\ldots,4$,\n",
    "and\n",
    "according to the local sensitivity measure $S_{Z_i}^{p}$ in\n",
    "([9](#eq:Sp)) all\n",
    "the input factors $Z_i$s are equally important and\n",
    "independent of the variation\n",
    "of each factor. This measure is clearly\n",
    "at odds with the ranking of influence\n",
    "based on the scatterplots in\n",
    "([8](#eq:scatter_plot_rank)) and is an indication\n",
    "of the usefulness of\n",
    "scatterplots in sensitivy analysis. However, the\n",
    "bidimensional\n",
    "scatterplots may in some cases be deceiving and lead to type II\n",
    "errors (i.e. failure to identify influential parameters) [Saltelli et al. 2004](https://www.wiley.com/en-us/Sensitivity+Analysis+in+Practice%3A+A+Guide+to+Assessing+Scientific+Models-p-9780470870938)\n",
    "\n",
    "Most sensitivity measures aim to preserve the rich information\n",
    "provided\n",
    "by the scatterplots in a condensed format. The challenge is\n",
    "how to rank the\n",
    "factors rapidly and automatically without having to\n",
    "inspect many scatterplots in\n",
    "situations with many input\n",
    "factors. Another challenge with scatterplots is that\n",
    "sensitivities for\n",
    "sets cannot be visualized, while luckily compact sensitivity\n",
    "measures may be\n",
    "defined in such cases."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "15"
    }
   },
   "outputs": [],
   "source": [
    "# start widgets for scatter plots\n",
    "from ipywidgets import FloatSlider, VBox\n",
    "import ipywidgets\n",
    "\n",
    "# Create array of plots\n",
    "f, ax = plt.subplots(int(Nrv/2), int(Nrv/2),sharey=True)\n",
    "\n",
    "weight_sliders=[]\n",
    "for k in range(Nrv):\n",
    "    wtitle = 'w' + str(k)\n",
    "    weight_sliders.append(FloatSlider(description=wtitle, min=1, max=5, value=2))\n",
    "\n",
    "\n",
    "def update_computations_and_plot(w0,w1,w2,w3):\n",
    "    w=np.array([w0,w1,w2,w3])\n",
    "    \n",
    "    Y = linear_model(w, Z.transpose())\n",
    "    \n",
    "    for k in range(Nrv):\n",
    "        Ztitle = 'Z' + str(k)\n",
    "        if (k<Nrv/2):\n",
    "            ax[0,k%2].clear()\n",
    "            ax[0,k%2].scatter(Z[k,:],Y[:])\n",
    "        else:\n",
    "            ax[1,k%2].clear()\n",
    "            ax[1,k%2].scatter(Z[k,:],Y[:])\n",
    "\n",
    "\n",
    "ipywidgets.interactive(update_computations_and_plot,w0=weight_sliders[0],\n",
    "                       w1=weight_sliders[1],w2=weight_sliders[2],w3=weight_sliders[3])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Normalized derivatives\n",
    "\n",
    "A simple way to improve the derivative sensitivity\n",
    "measure $S_{Z_i}^{p}$ in\n",
    "([9](#eq:Sp)) is to scale the input-output variables\n",
    "with their standard deviations:\n",
    "\n",
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"eq:Ss\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "S_{Z_i}^{\\sigma} = \\frac{\\partial Y/\\sigma_Y}{\\partial\n",
    "Z_i/\\sigma_{Z_i}} = \\frac{\\sigma_{Z_i}}{\\sigma_{Y}} \\; \\frac{\\partial\n",
    "Y}{\\partial Z_i}\n",
    "\\label{eq:Ss} \\tag{10}\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "In case of our simple linear model ([1](#eq:linear_model)) we get from\n",
    "([10](#eq:Ss)):\n",
    "\n",
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"eq:Ss_simple\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\left (S_{Z_i}^{\\sigma} \\right)^2 = \\left(\n",
    "\\frac{\\sigma_{Z_i}}{\\sigma_{Y}}\\right)^2 \\; \\left (\\frac{\\partial Y}{\\partial\n",
    "Z_i}\\right)^2 = \\left( \\frac{\\sigma_{Z_i}\\, \\Omega_i}{\\sigma_{Y}}\\right)^2 \\;\n",
    "\\qquad \\textsf{which may be rearranged to:} \\qquad \\sigma_y^2 \\,\n",
    "(S_{Z_i}^{\\sigma})^2 = \\left ( \\Omega_{i} \\sigma_{Y} \\right )^2\n",
    "\\label{eq:Ss_simple} \\tag{11}\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "Based on the linearity of our model we previously found\n",
    "([4](#eq:analytic_mean_std)) which also yields:\n",
    "\n",
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"eq:Ss_model_ded\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    " \\sigma_Y^2 = \\sum_{i=1}^{r} \\left(\\Omega_i^2 \\,\n",
    "\\sigma_{Z_i}\\right)^2\n",
    "\\label{eq:Ss_model_ded} \\tag{12}\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "As both ([12](#eq:Ss_model_ded)) and ([11](#eq:Ss_simple)) must hold\n",
    "simultaneously we get\n",
    "\n",
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"eq:Ss1\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\left (S_{Z_i}^{\\sigma} \\right)^2=1 \n",
    "\\label{eq:Ss1} \\tag{13}\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "The normalized derivative measure of sensitivity in ([10](#eq:Ss)) is\n",
    "more\n",
    "convincing than ([9](#eq:Sp)): first, as it involves both the\n",
    "weights $\\Omega_i$\n",
    "and the factors $Z_i$ in ([1](#eq:linear_model));\n",
    "second as the measures are\n",
    "properly scaled and sums up to one,\n",
    "which allows for an easy interpretation\n",
    "of the output sensitivity with\n",
    "respect to each of the input factors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "7"
    }
   },
   "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",
    "    import pandas as pd\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 subsequent model\n",
    "evaluations, we may estimate the standard deviation\n",
    "of $\\mathbf{Y}$ and compute\n",
    "the relative error with respect to the\n",
    "theoretical value. You may change the\n",
    "number of sample above,\n",
    "i.e. $N$, and see how $N$ influence the estimates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "8"
    }
   },
   "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\n",
    "number of times we have evaluated our simple linear model\n",
    "([1](#eq:linear_model)). The evaluation of the model is normally the\n",
    "most\n",
    "computationally expensive part of the analysis, and for that\n",
    "reasons `Ns` is\n",
    "referred to as the `cost` of the analysis.\n",
    "\n",
    "# Conditional variances\n",
    "\n",
    "As noted\n",
    "previously, the importance of a factor $Z_i$ is manifested\n",
    "the existence of a\n",
    "`shape` or `pattern` in the scatter plot of model outputs\n",
    "$Y$ against $Z_i$. Conversely, a uniform cloud of\n",
    "output points $Y$ as a function of\n",
    "$Z_i$ is a symptom - not a proof -\n",
    "indicating that $Z_i$ is a\n",
    "noninfluential factor. In this section we seek to\n",
    "demonstrate that\n",
    "conditional variances is a useful means to quantify the\n",
    "`shape` or\n",
    "`pattern` in the outputs.\n",
    "\n",
    "The shape in the outputs $Y$ for a given\n",
    "$Z_i$, may be seen in the\n",
    "scatterplot as of $Y$ versus $Z_i$. In particular, we\n",
    "may cut the\n",
    "$Z_i$-axis into slices and assess how the distribution of the\n",
    "outputs\n",
    "$Y$ changes from slice to slice. This is illustrated in the code\n",
    "snippet\n",
    "below, where the slices are identified with vertical dashed\n",
    "lines at equidistant\n",
    "locations on each $Z_i$-axis, $i=1, \\ldots,4$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "9"
    }
   },
   "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\n",
    "$Z_i$ fixed while averaging over all output values of $Y$ due\n",
    "to all-but $Z_i$,\n",
    "which corresponds to the conditional expected value:\n",
    "\n",
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto3\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "E_{Z_{\\sim i}} (Y\\;|\\;Z_i) \n",
    "\\label{_auto3} \\tag{14}\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "For convenience we let $Z_{\\sim i}$ denote `all-but` $Z_i$. Naturally,\n",
    "a measure\n",
    "of how much $E_{Z_{\\sim i}} (Y\\;|\\;Z_i)$ varies in the range\n",
    "of $Z_i$ is given\n",
    "by the conditional variance:\n",
    "\n",
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto4\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\text{V}_{Z_i}(E_{Z_{\\sim i}} (Y\\;|\\;Z_i))\n",
    "\\label{_auto4}\n",
    "\\tag{15}\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "Further, the variance the output $Y$ may be decomposed into:\n",
    "\n",
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"eq:VarDecomp\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\text{V}(Y) = E_{Z_i} ( V_{Z_{\\sim i}} (Y \\; | Z_{i})) +\n",
    "\\text{V}_{Z_i}(E_{Z_{\\sim i}} (Y\\;|\\;Z_i))\n",
    "\\label{eq:VarDecomp} \\tag{16}\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "A large $\\text{V}_{Z_i}(E_{Z_{\\sim i}} (Y\\;|\\;Z_i))$ will imply that\n",
    "$Z_i$ is an\n",
    "important factor and is therefore coined the first-order\n",
    "effect of $Z_i$ on $Y$,\n",
    "and its fraction of the total variation of $Y$ is expressed by $S_i$, `the\n",
    "first-order sensitivity index` of $Z_i$ on $Y$:\n",
    "\n",
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto5\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "S_i = \\frac{\\text{V}_{Z_i}(E_{Z_{\\sim i}}\n",
    "(Y\\;|\\;Z_i))}{\\text{V}(Y)}\n",
    "\\label{_auto5} \\tag{17}\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "By ([16](#eq:VarDecomp)), $S_i$ is number always in the range $[0,1]$,\n",
    "and a\n",
    "high value implies an important factor."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "10"
    }
   },
   "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\n",
    "sensitivity indices may be\n",
    "computed with two approaches; the Monte Carlo method\n",
    "and the\n",
    "polynomial chaos expansion method.\n",
    "\n",
    "### Monte Carlo\n",
    "\n",
    "Below some code\n",
    "snippets are provided to illustrate how we may compute\n",
    "the Sobol indices with\n",
    "the MCM. For the interested reader we have also\n",
    "written a seperate and more\n",
    "detailed notebook [A brief introduction to\n",
    "UQ and SA with the Monte Carlo\n",
    "method](monte_carlo.ipynb)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "11"
    }
   },
   "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": 12,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "12"
    }
   },
   "outputs": [],
   "source": [
    "# Monte Carlo\n",
    "# get joint distributions\n",
    "    from sensitivity_examples_nonlinear import generate_distributions\n",
    "    \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\n",
    "to illustrate how we may compute\n",
    "the Soboil indices with the polynomial chaos\n",
    "expansions using `chaospy`. A more in-depth view on `chaospy` and its usage\n",
    "can be found in the notebook [A practical introduction to polynomial\n",
    "chaos with the chaospy package](introduction_gpc.ipynb)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "13"
    }
   },
   "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": 14,
   "metadata": {
    "attributes": {
     "classes": [],
     "id": "",
     "n": "14"
    }
   },
   "outputs": [],
   "source": [
    "# Polychaos convergence\n",
    "    from numpy import linalg as LA\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')"
   ]
  }
 ],
 "metadata": {},
 "nbformat": 4,
 "nbformat_minor": 2
}