{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- dom:TITLE: Higher order sensitivity indices for interaction models -->\n",
    "# Higher order sensitivity indices for interaction models\n",
    "<!-- dom:AUTHOR: Vinzenz Gregor Eck at Expert Analytics -->\n",
    "<!-- Author: -->  \n",
    "**Vinzenz Gregor Eck**, Expert Analytics  \n",
    "<!-- dom:AUTHOR: Jacob Sturdy at Department of Structural Engineering, NTNU -->\n",
    "<!-- Author: --> **Jacob Sturdy**, Department of Structural Engineering, NTNU  \n",
    "<!-- dom:AUTHOR: Leif Rune Hellevik at Department of Structural Engineering, NTNU -->\n",
    "<!-- Author: --> **Leif Rune Hellevik**, Department of Structural Engineering, NTNU\n",
    "\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)\n",
    "\n",
    "# font = {'family' : 'sans-serif',\n",
    "#         'weight' : 'normal',\n",
    "#         'size'   : 18.0}\n",
    "# matplotlib.rc('font', **font)  # pass in the font dict as kwar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# import modules\n",
    "import numpy as np\n",
    "import chaospy as cp\n",
    "from monte_carlo import generate_sample_matrices_mc\n",
    "from monte_carlo import calculate_sensitivity_indices_mc\n",
    "from xlwt.Utils import col_by_name"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# A model with interaction\n",
    "<div id=\"sec:non_additive_model\"></div>\n",
    "\n",
    "For convenience and simplicity, we consider the same model as before:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"eq:non_add_model\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "Y = \\sum_{i=1}^{r} \\Omega_i \\, Z_i\n",
    "\\label{eq:non_add_model} \\tag{1}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The only difference is that we render ([1](#eq:non_add_model))\n",
    "non-additive by letting both $\\Omega_i$ and $Z_i$ be random\n",
    "variables. By doing so our model will have products of two distinct\n",
    "random variables, and consequently it will be a non-additive model, as it no longer fulfills the criterion:\n",
    "\n",
    "\n",
    "\n",
    "<hr/>\n",
    "**Definition: Additive model.**\n",
    "\n",
    "An [additive model](https://chemicalstatistician.wordpress.com/2014/03/07/applied-statistics-lesson-of-the-day-additive-models-vs-interaction-models-in-2-factor-experimental-designs/) is the arithmetic sum of its individual random variables.\n",
    "<hr/>\n",
    "\n",
    "\n",
    "\n",
    "In constrast, a `non-additive model`, is `not` simply a sum of its individual\n",
    "random variables. For our model ([1](#eq:non_add_model)) the effect of\n",
    "$Z_i$ on $Y$ depends on the value of $\\Omega_i$, and for this reason\n",
    "non-additive moldels are also often called `interaction models`.\n",
    "\n",
    "Note that $ Y = \\sum_{i=1}^{r} \\Omega_i \\, Z_i^2$, will be an additive\n",
    "model for constant $\\Omega_i$, but non-linear the individual $Z_i$.\n",
    "\n",
    "For our non-additive (interaction) model ([1](#eq:non_add_model)) we\n",
    "assume $\\Omega_i$ and $Z_i$ to be normaly distributed random\n",
    "variables:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto1\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "Z_i  \\sim N(0, \\sigma_{Z_i}), \\qquad i=1,2, \\ldots, r \n",
    "\\label{_auto1} \\tag{2}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto2\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation} \n",
    "\\Omega_i  \\sim N(\\mu_{\\Omega_i}, \\sigma_{\\Omega_i})\n",
    "\\label{_auto2} \\tag{3}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Which means that the distributions of $Z_i$ are unchanged, whereas the\n",
    "mean values $\\mu_{\\Omega_i}$ of the distributions of $\\Omega_i$ are\n",
    "chosen to be non-zero:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto3\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\mu_{\\Omega_i} = 0.5 i, \\qquad i=1,2, \\ldots, r \\qquad \\textsf{such that} \\qquad\n",
    "\\bar{\\Omega}_1 < \\bar{\\Omega}_3 < \\bar{\\Omega}_3 <\\bar{\\Omega}_4\n",
    "\\label{_auto3} \\tag{4}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The number of random variables for our interaction model\n",
    "([1](#eq:non_add_model)) has increased to $N_{rv} = 2\\, r$ and we arrange them in a vector in the following manner:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto4\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\mathbf{X}= \\left ( Z_1, Z_2, \\ldots, Z_r, \\Omega_1, \\Omega_2, \\ldots, \\Omega_r \\right)\n",
    "\\label{_auto4} \\tag{5}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below we will demonstrate how the Sobol sensitivity indices may be\n",
    "computed with both the Monte Carlo Method and the Polynomial Chaos\n",
    "method. The python code for the model is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "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": [
    "### Definition of the random input variables\n",
    "\n",
    "In this example of interaction models we set $r=4$ in\n",
    "([1](#eq:non_add_model)) and in the code snippet below we let the\n",
    "numpy-arrays `zm` and `wm` hold the expectations in the first column\n",
    "and the standard deviation in the second column of $Z_i$ and\n",
    "$\\Omega_i$, repectively. You may uncomment the lines in the code\n",
    "snippet to set values of e.g. the standard deviations yourself."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "    # definition of mu and sig for z and w\n",
    "    N_terms = 4\n",
    "    c = 0.5\n",
    "    zm = np.array([[0., i] for i in range(1, N_terms+1)])\n",
    "    wm = np.array([[i * c, i] for i in range(1, N_terms+1)])\n",
    "\n",
    "    # to see the effect of changing the values in zm uncomment and change one of these lines and re-run\n",
    "\n",
    "    # zm[0, 1] = 1\n",
    "    # zm[1, 1] = 20\n",
    "    # zm[2, 1] = 3\n",
    "    # zm[3, 1] = 10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example we choose to set $r=4$ and use the chaospy package to\n",
    "define our random input variables for our model\n",
    "([1](#eq:non_add_model)). The function `generate_distributions` takes\n",
    "two numpy-arrays `zm` and `wm` as arguments, where `wm` is\n",
    "optional. Whenever `wm` is passed to the function (not `None`) it is\n",
    "treated as an array holding the expectation and standard deviations as\n",
    "`zm` and appended to the `xm` array to hold all the expectations and\n",
    "standard deviations for both $Z_i$ and $\\Omega_i$. The marginal\n",
    "distributions (or probability density functions or pdf for short) for each element in `xm` are generated with repeated calls to chaospy for each elmement in `xm`: `[cp.Normal(*mu_sig) for mu_sig in xm]`. Note that the `*-operator` is used to unpack the arguments out of a list or tuple.   \n",
    "\n",
    "Finally, we generate the the joint pdf by unpacking all the elements\n",
    "in the list `marginal distributions` and passing them to `cp.J`, the\n",
    "joint probability density function of Chaospy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# definition of the distributions\n",
    "def generate_distributions(zm, wm=None):\n",
    "    # define marginal distributions\n",
    "    if wm is not None:\n",
    "        zm = np.append(zm, wm, axis=0)\n",
    "    marginal_distributions = [cp.Normal(*mu_sig) for mu_sig in zm]\n",
    "    # define joint distributions\n",
    "    jpdf = cp.J(*marginal_distributions)\n",
    "\n",
    "    return jpdf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Scatter plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "    # generate the joint distribution\n",
    "    jpdf = generate_distributions(zm, wm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "    # Scatter plots of data for visual inspection of sensitivity\n",
    "    N_plot=100\n",
    "    N_prms = len(jpdf)\n",
    "    N_terms = N_prms//2\n",
    "\n",
    "    Xs=jpdf.sample(N_plot,sample_method='R').transpose()  \n",
    "    Zs = Xs[:, :N_terms]  # Split X in two vectors for X and W\n",
    "    Ws = Xs[:, N_terms:]\n",
    "    Ys = linear_model(Ws, Zs)\n",
    "\n",
    "    scatter = plt.figure('Scatter plots')\n",
    "    for k in range(N_terms):\n",
    "        plt.subplot(2, N_terms, k + 1)\n",
    "        _=plt.plot(Zs[:, k], Ys[:], '.')\n",
    "        if k==0 : plt.label('Y')\n",
    "        plt.xlabel('Z{}'.format(k+1))\n",
    "        plt.ylim([-150, 150])\n",
    "\n",
    "        plt.subplot(2, N_terms, k + 1 + N_terms)\n",
    "        _=plt.plot(Ws[:, k], Ys[:], '.')\n",
    "        if k==0 : plt.label('Y')\n",
    "        plt.xlabel('W{}'.format(k+1))\n",
    "        plt.ylim([-150, 150])\n",
    "    scatter.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# First order sensitivity coefficients\n",
    "\n",
    "### Analytical computation of the sensitivity indices\n",
    "\n",
    "The calculation of the sensitivity indices for this model can be done analytically.\n",
    "We will use this as reference for the comparison with the results of the Monte Carlo and polynomial chaos method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "    # calculate the analytic sensitivity coefficients\n",
    "\n",
    "    VarY = np.sum(zm[:, 1]**2 * (wm[:, 0]**2 + wm[:, 1]**2), axis=0)\n",
    "    Sz = wm[:, 0]**2 * zm[:, 1]**2/VarY # first order indices\n",
    "    Sw = np.zeros_like(Sz)\n",
    "    Szw= wm[:, 1]**2 * zm[:, 1]**2/VarY  # second order indices\n",
    "    StZ = (wm[:, 0]**2 * zm[:, 1]**2 + wm[:, 1]**2 * zm[:, 1]**2)/VarY # total indices\n",
    "    Stw = (wm[:, 1]**2 * zm[:, 1]**2)/VarY\n",
    "\n",
    "    # join sensitivity arrays\n",
    "    Sa = np.append(Sz, Sw)\n",
    "    Sta = np.append(StZ, Stw)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we will calculate of the sensitivity indices with both the Monte Carlo and polynomial chaos method:\n",
    "\n",
    "### Define functions for computing the sensitivity indices with the MCM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# calculate sens indices of non additive model\n",
    "def monte_carlo_sens_nonlin(Ns, jpdf, sample_method='R'):\n",
    "\n",
    "    N_prms = len(jpdf)\n",
    "\n",
    "    # 1. Generate sample matrices\n",
    "    XA, XB, XC = generate_sample_matrices_mc(Ns, N_prms, jpdf, sample_method)\n",
    "\n",
    "    # 2. Evaluate the model\n",
    "    Y_A, Y_B, Y_C = evaluate_non_additive_linear_model(XA, XB, XC)\n",
    "\n",
    "    # 3. Approximate the sensitivity indices\n",
    "    S, ST = calculate_sensitivity_indices_mc(Y_A, Y_B, Y_C)\n",
    "\n",
    "    return XA, XB, XC, 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_non_additive_linear_model(X_A, X_B, X_C):\n",
    "\n",
    "    N_prms = X_A.shape[1]\n",
    "    Ns = X_A.shape[0]\n",
    "    N_terms = int(N_prms / 2)\n",
    "    # 1. evaluate sample matrices X_A\n",
    "    Z_A = X_A[:, :N_terms]  # Split X in two vectors for X and W\n",
    "    W_A = X_A[:, N_terms:]\n",
    "    Y_A = linear_model(W_A, Z_A)\n",
    "\n",
    "    # 2. evaluate sample matrices X_B\n",
    "    Z_B = X_B[:, :N_terms]\n",
    "    W_B = X_B[:, N_terms:]\n",
    "    Y_B = linear_model(W_B, Z_B)\n",
    "\n",
    "    # 3. evaluate sample matrices X_C\n",
    "    Y_C = np.empty((Ns, N_prms))\n",
    "    for i in range(N_prms):\n",
    "        x = X_C[i, :, :]\n",
    "        z = x[:, :N_terms]\n",
    "        w = x[:, N_terms:]\n",
    "        Y_C[:, i] = linear_model(w, z)\n",
    "\n",
    "    return Y_A, Y_B, Y_C"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define functions for computing the sensitivity indices with polynomial chaos expansions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# polynomial chaos sensitivity analysis\n",
    "def polynomial_chaos_sens(Ns_pc, jpdf, polynomial_order, poly=None, return_reg=False):\n",
    "    N_terms = int(len(jpdf) / 2)\n",
    "    # 1. generate orthogonal polynomials\n",
    "    poly = poly or cp.orth_ttr(polynomial_order, jpdf)\n",
    "    # 2. generate samples with random sampling\n",
    "    samples_pc = jpdf.sample(size=Ns_pc, rule='R')\n",
    "    # 3. evaluate the model, to do so transpose samples and hash input data\n",
    "    transposed_samples = samples_pc.transpose()\n",
    "    samples_z = transposed_samples[:, :N_terms]\n",
    "    samples_w = transposed_samples[:, N_terms:]\n",
    "    model_evaluations = linear_model(samples_w, samples_z)\n",
    "    # 4. calculate generalized polynomial chaos expression\n",
    "    gpce_regression = cp.fit_regression(poly, samples_pc, model_evaluations)\n",
    "    # 5. get sensitivity indices\n",
    "    Spc = cp.Sens_m(gpce_regression, jpdf)\n",
    "    Stpc = cp.Sens_t(gpce_regression, jpdf)\n",
    "    \n",
    "    if return_reg:\n",
    "        return Spc,Stpc,gpce_regression    \n",
    "    else:\n",
    "        return Spc, Stpc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Perform the computations\n",
    "\n",
    "**MCM**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "    # Monte Carlo\n",
    "    #Ns_mc = 1000000 # Number of samples mc\n",
    "    Ns_mc = 10000 # Number of samples mc\n",
    "    # calculate sensitivity indices with mc\n",
    "    A_s, B_s, C_s, f_A, f_B, f_C, Smc, Stmc = monte_carlo_sens_nonlin(Ns_mc, jpdf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Polynomial chaos**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "    # compute with Polynomial Chaos\n",
    "    Ns_pc = 200\n",
    "    polynomial_order = 3\n",
    "    \n",
    "    # calculate sensitivity indices with gpc\n",
    "    Spc, Stpc, gpce_reg = polynomial_chaos_sens(Ns_pc, jpdf, polynomial_order,return_reg=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Compare the computations**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "    # compare the computations\n",
    "    import pandas as pd\n",
    "    row_labels  = ['X_'+str(x) for x in range(1,N_terms*2+1)]\n",
    "    S=np.column_stack((Sa,Spc,Smc,Sta,Stpc,Stmc))\n",
    "    S_table = pd.DataFrame(S, columns=['Sa','Spc','Smc','Sta','Stpc','Stmc'], index=row_labels)  \n",
    "    print(S_table.round(3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Second and higher order sensitivity coefficients\n",
    "\n",
    "In the [sensitivity introduction notebook](sensitivity_introduction.ipynb) we introduced the first order\n",
    "sensitivity indices as:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- 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}} (Y\\;|\\;Z_i))}{\\text{V}(Y)}\n",
    "\\label{_auto5} \\tag{6}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Based on the seminal works of the Russian mathematician Sobol, we may\n",
    "introduce decompositions of the form:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"eq:Vi\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "V_i = V \\left ( E(Y\\,|\\, X_i) \\right ) \n",
    "\\label{eq:Vi} \\tag{7} \n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"eq:Vij\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation} \n",
    "V_{ij} = V \\left ( E(Y\\,|\\, X_i, X_j) \\right )  - V_i - V_j \n",
    "\\label{eq:Vij} \\tag{8} \n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where we for convenience have dropped the subscripts for the\n",
    "expectation oprerator and the variance operator. Keep in mind\n",
    "that the outer operator always operates on the conditioned factors, in\n",
    "this case $X_i$ and $X_j$.\n",
    "\n",
    "   * The second order conditional variance $V \\left ( E(Y\\,|\\, X_i, X_j) \\right )$ is a measure of the combined effect of the pair $(X_i,X_j)$.\n",
    "\n",
    "   * The $V_{ij}$ is interaction term between the factors $X_i$ and $X_j$ and the combined effect minus the first order effetcs $V_i$ and $V_j$.\n",
    "\n",
    "   * The $V_{ij}$ captures the part of $V(Y)$ which cannot be expressed as a sum of the separate effects of $X_i$ and $X_j$\n",
    "\n",
    "By repeated introduction of progressively higher order conditional\n",
    "variences, the total variance of the output $V(Y)$ may be presented as\n",
    "a ANOVA-HDMR decomposition:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto6\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "V(Y) = \\sum_{i} V_i +\\sum_{i}\\sum_{j>i} V_{ij}+ \\sum_{i}\\sum_{j>i}\\sum_{l>j} V_{ijl}+\\ldots + V_{12...k}\n",
    "\\label{_auto6} \\tag{9}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "which by division of $V(Y)$ yields the very useful relation:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"eq:S_decomp\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\sum_{i} S_i +\\sum_{i}\\sum_{j>i} S_{ij}+ \\sum_{i}\\sum_{j>i}\\sum_{l>j} S_{ijl}+\\ldots + S_{12...k} =1 \n",
    "\\label{eq:S_decomp} \\tag{10}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Higher order sensitivity indices**\n",
    "\n",
    "The variance-based sensitivity indices in ([10](#eq:S_decomp)) represent\n",
    "a way to express 100% of the variance of $Y$ regardless of whether $Y$\n",
    "is additive or non-additive. In other words, ([10](#eq:S_decomp))\n",
    "provides a complete quantification of the sensitivity of the model\n",
    "with respect to its predictive factors. The only problem is that\n",
    "([10](#eq:S_decomp)) may have as many as $2^k-1$ terms for $k$ input factors.\n",
    "\n",
    "  * For $k=3$ this gives $2^3-1=7$ indices: $S_1$, $S_2$, $S_3$,$S_{12}$, $S_{13}$, $S_{23}$, $S_{123}$,\n",
    "\n",
    "  * For $k=10$ we get $2^{10}-1=1023$ indices, which is neither tractable and informative in practice. \n",
    "\n",
    "Of the $2^k-1$ terms we have:\n",
    "  * $k$ first order indices $S_i$\n",
    "\n",
    "  * $ {k\\choose 2} =\\frac{k!}{2!\\,(k-2)!}$ second order indices $S_{ij}$\n",
    "\n",
    "  * $ {k\\choose 3}$ third order indices $S_{ijk}$\n",
    "\n",
    "  * ..and so on\n",
    "\n",
    "The computation of the second order indices for our model example is\n",
    "illustrated in the code snippet below. For our simple model,\n",
    "interaction terms higher then second order will be zero. Consequently,\n",
    "by ([10](#eq:S_decomp)) the sum of the first and second order indicecs\n",
    "should sum up to unity. We test that this is indeed the case by\n",
    "printing the sum. Further, we compute the error for the indices\n",
    "computed with `chaospy` for each individual index."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "    # Second order indices with gpc\n",
    "    \n",
    "    S2 = cp.Sens_m2(gpce_reg, jpdf) # second order indices with gpc\n",
    "    \n",
    "    # print all second order indices\n",
    "    print(pd.DataFrame(S2,columns=row_labels,index=row_labels).round(3))\n",
    "    \n",
    "    # sum all second order indices \n",
    "    SumS2=np.sum(np.triu(S2))\n",
    "    print('\\nSum Sij = {:2.2f}'.format(SumS2))\n",
    "    \n",
    "    # sum all first and second order indices\n",
    "    print('Sum Si + Sij = {:2.2f}\\n'.format(np.sum(Spc)+SumS2))\n",
    "    \n",
    "    # compare nonzero second order indices with analytical indices \n",
    "    Szw_pc=[S2[i,i+N_terms] for i in range(N_terms) ]\n",
    "    Szw_table=np.column_stack((Szw_pc,Szw,(Szw_pc-Szw)/Szw))\n",
    "    print(pd.DataFrame(Szw_table,columns=['Szw','Szw pc','Error%']).round(3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The important thing to realize is that the number of interaction terms\n",
    "increase exponentialy as the number of input factors $k$ increase.  In\n",
    "this case, the remedy offered by variance based analysis is to compute\n",
    "the 'total effect term', which we will outline in the following.\n",
    "\n",
    "\n",
    "\n",
    "# Convergence Analysis\n",
    "## Monte Carlo Method\n",
    "This may take a while. Be patient!\n",
    "(You can also reduce iteration number for averaging.)\n",
    "\n",
    "Convergence analysis for the monte carlo method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "        # Convergence Monte Carlo with random sampling\n",
    "        list_of_samples = np.array([10000, 50000, 100000, 500000, 1000000])\n",
    "        s_mc_err = np.zeros((len(list_of_samples), N_prms))\n",
    "        st_mc_err = np.zeros((len(list_of_samples), N_prms))\n",
    "        # average over\n",
    "        N_iter = 5\n",
    "        print('MC convergence analysis:')\n",
    "        for i, N_smpl in enumerate(list_of_samples):\n",
    "            print('    N_smpl {}'.format(N_smpl))\n",
    "            for j in range(N_iter):\n",
    "                A_s, XB, XC, Y_A, Y_B, Y_C, S, ST = monte_carlo_sens_nonlin(N_smpl,\n",
    "                                                                                jpdf,\n",
    "                                                                                sample_method='R')\n",
    "                s_mc_err[i] += np.abs(S - Sa)\n",
    "                st_mc_err[i] += np.abs(ST - Sta)\n",
    "                print('         finished with iteration {} of {}'.format(1 + j, N_iter))\n",
    "            s_mc_err[i] /= float(N_iter)\n",
    "            st_mc_err[i] /= float(N_iter)\n",
    "        # Plot results for monte carlo\n",
    "        fig_random = plt.figure('Random sampling - average of indices')\n",
    "        fig_random.suptitle('Random sampling - average of indices')\n",
    "\n",
    "        ax = plt.subplot(1, 2, 1)\n",
    "        plt.title('First order sensitivity indices')\n",
    "        _=plt.plot(list_of_samples / 1000, np.sum(s_mc_err, axis=1), '-')\n",
    "        ax.set_yscale('log')\n",
    "        _=plt.ylabel('abs error')\n",
    "        _=plt.xlabel('number of samples [1e3]')\n",
    "\n",
    "        ax1 = plt.subplot(1, 2, 2)\n",
    "        plt.title('Total sensitivity indices')\n",
    "        _=plt.plot(list_of_samples / 1000, np.sum(st_mc_err, axis=1), '-')\n",
    "        ax1.set_yscale('log')\n",
    "        _=plt.ylabel('abs error')\n",
    "        _=plt.xlabel('number of samples [1e3]')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plotting the individual sensitivity indices:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "        # Plot results for monte carlo figure individual\n",
    "        fig_random = plt.figure('Random sampling')\n",
    "        fig_random.suptitle('Random sampling')\n",
    "        for l, (s_e, st_e) in enumerate(zip(s_mc_err.T, st_mc_err.T)):\n",
    "            ax = plt.subplot(1, 2, 1)\n",
    "            plt.title('First order sensitivity indices')\n",
    "            plt.plot(list_of_samples / 1000, s_e, '-', label='S_{}'.format(l))\n",
    "            ax.set_yscale('log')\n",
    "            _=plt.ylabel('abs error')\n",
    "            _=plt.xlabel('number of samples [1e3]')\n",
    "            _=plt.legend()\n",
    "\n",
    "            ax1 = plt.subplot(1, 2, 2)\n",
    "            plt.title('Total sensitivity indices')\n",
    "            _=plt.plot(list_of_samples / 1000, st_e, '-', label='ST_{}'.format(l))\n",
    "            ax1.set_yscale('log')\n",
    "            _=plt.ylabel('abs error')\n",
    "            _=plt.xlabel('number of samples [1e3]')\n",
    "            plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Polynomial Chaos\n",
    "(Note how much faster this is.)\n",
    "Convergence analysis for the polynomial chaos method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "        # Convergence Polynomial Chaos\n",
    "        list_of_samples = np.array([140, 160, 200, 220])\n",
    "        s_pc_err = np.zeros((len(list_of_samples), N_prms))\n",
    "        st_pc_err = np.zeros((len(list_of_samples), N_prms))\n",
    "        polynomial_order = 3\n",
    "        # average over\n",
    "        N_iter = 4\n",
    "        print('PC convergence analysis:')\n",
    "        poly = cp.orth_ttr(polynomial_order, jpdf)\n",
    "        for i, N_smpl in enumerate(list_of_samples):\n",
    "            print('    N_smpl {}'.format(N_smpl))\n",
    "            for j in range(N_iter):\n",
    "                # calculate sensitivity indices\n",
    "                Spc, Stpc = polynomial_chaos_sens(N_smpl, jpdf, polynomial_order, poly)\n",
    "                s_pc_err[i] += np.abs(Spc - Sa)\n",
    "                st_pc_err[i] += np.abs(Stpc - Sta)\n",
    "                print('         finished with iteration {} of {}'.format(1 + j, N_iter))\n",
    "            s_pc_err[i] /= float(N_iter)\n",
    "            st_pc_err[i] /= float(N_iter)\n",
    "\n",
    "        # Plot results for polynomial chaos\n",
    "        fig_random = plt.figure('Polynomial Chaos - average of indices')\n",
    "        fig_random.suptitle('Polynomial Chaos - average of indices')\n",
    "\n",
    "        ax = plt.subplot(1, 2, 1)\n",
    "        plt.title('First order sensitivity indices')\n",
    "        _=plt.plot(list_of_samples, np.sum(s_pc_err, axis=1), '-')\n",
    "        ax.set_yscale('log')\n",
    "        _=plt.ylabel('abs error')\n",
    "        _=plt.xlabel('number of samples [1e3]')\n",
    "\n",
    "        ax1 = plt.subplot(1, 2, 2)\n",
    "        plt.title('Total sensitivity indices')\n",
    "        _=plt.plot(list_of_samples, np.sum(st_pc_err, axis=1), '-')\n",
    "        ax1.set_yscale('log')\n",
    "        _=plt.ylabel('abs error')\n",
    "        _=plt.xlabel('number of samples [1e3]')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plotting the individual sensitivity indices:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "        # Plot results for polynomial chaos individual\n",
    "        fig_random = plt.figure('Polynomial Chaos')\n",
    "        fig_random.suptitle('Polynomial Chaos')\n",
    "        for l, (s_e, st_e) in enumerate(zip(s_pc_err.T, st_pc_err.T)):\n",
    "            ax = plt.subplot(1, 2, 1)\n",
    "            plt.title('First order sensitivity indices')\n",
    "            _=plt.plot(list_of_samples, s_e, '-', label='S_{}'.format(l))\n",
    "            ax.set_yscale('log')\n",
    "            plt.ylabel('abs error')\n",
    "            plt.xlabel('number of samples [1e3]')\n",
    "            plt.legend()\n",
    "\n",
    "            ax1 = plt.subplot(1, 2, 2)\n",
    "            plt.title('Total sensitivity indices')\n",
    "            _=plt.plot(list_of_samples, st_e, '-', label='ST_{}'.format(l))\n",
    "            ax1.set_yscale('log')\n",
    "            plt.ylabel('abs error')\n",
    "            plt.xlabel('number of samples [1e3]')\n",
    "            plt.legend()"
   ]
  }
 ],
 "metadata": {},
 "nbformat": 4,
 "nbformat_minor": 2
}