{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "# A brief introduction to UQ and SA with the Monte Carlo method\n", "\n", " \n", "**Vinzenz Gregor Eck**, Expert Analytics \n", "\n", " **Leif Rune Hellevik**, 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": true }, "outputs": [], "source": [ "import numpy as np\n", "import chaospy as cp\n", "import monte_carlo\n", "from sensitivity_examples_nonlinear import generate_distributions\n", "from sensitivity_examples_nonlinear import monte_carlo_sens_nonlin\n", "from sensitivity_examples_nonlinear import analytic_sensitivity_coefficients\n", "from sensitivity_examples_nonlinear import polynomial_chaos_sens" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Monte Carlo\n", "\n", "The Monte Carlo method (MCM) is probably the most widely applied method for\n", "variance based uncertainty quantification and sensitivity\n", "analysis. Monte carlo methods are generally straight forward to use\n", "and may be applied to a wide variety of problems as they require few\n", "assumptions about the model or quantity of interest and require no\n", "modifications of the model itself, i.e. the model may be used as a\n", "black box. The basic idea is to calculate statistics (mean, standard\n", "deviation, variance, sobol indices) of $Y$ directly from large amount\n", "of sample evaluations from the black box model $y$.\n", "\n", "\n", "\n", "
\n", "**Monte Carlo approach.**\n", "\n", "1. Sample a set of input samples $\\mathbf{z}^{(s)}$ from the input space $\\Omega_\\mathbf{Z}$ that is defined by the joint probability density function ${F_Z}$.\n", "\n", "2. Evaluate the deterministic model $y(\\mathbf{z})$ for each sample in $\\mathbf{z}^{(s)}$ to produce a set of model outputs $y^{(s)}$. \n", "\n", "3. Estimate all uncertainty measures and sensitivity indices from $y^{(s)}$.\n", "
\n", "\n", "\n", "\n", "For demonstration purposes we will use the same model as before:" ] }, { "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": [ "### Expectation and variance\n", "\n", "Once the model outputs have been computed the expectation and variance\n", "of the output are computed with the normal estimators." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", " {\\mathbb{E}}(Y) \\approx \\frac{1}{N} \\sum_{s=1}^{N} y^{(s)} \\qquad \\text{and} \\qquad \\operatorname{Var}(Y) \\approx \\frac{1}{N\\!-\\!1} \\sum_{s=1}^{N} \\left( y^{(s)} - {\\mathbb{E}}(Y)\\right)^2.\n", " \\label{eq:expected_value_MonteCarlo} \\tag{1}\n", " \\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Below we demonstrate how `chaospy` may be used for sampling and `numpy` for the statistics." ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [], "source": [ " # start uq\n", " # generate the distributions for the problem\n", " Nrv = 4\n", " c = 0.5\n", " zm = np.array([[0., i] for i in range(1, Nrv + 1)])\n", " wm = np.array([[i * c, i] for i in range(1, Nrv + 1)])\n", " jpdf = generate_distributions(zm, wm)\n", "\n", " # 1. Generate a set of Xs\n", " Ns = 20000\n", " Xs = jpdf.sample(Ns, rule='R').T # <- transform the sample matrix\n", "\n", " # 2. Evaluate the model\n", " Zs = Xs[:, :Nrv]\n", " Ws = Xs[:, Nrv:]\n", " Ys = linear_model(Ws, Zs)\n", "\n", " # 3. Calculate expectation and variance\n", " EY = np.mean(Ys)\n", " VY = np.var(Ys, ddof=1) # NB: use ddof=1 for unbiased variance estimator, i.e /(Ns - 1)\n", "\n", " print('E(Y): {:2.5f} and Var(Y): {:2.5f}'.format(EY, VY))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Variance based sensitivity measures\n", "\n", "In our [sensitivity_introduction notebook](sensitivity_introduction.ipynb) model we calculated the sensitivity\n", "coefficients with the MCM in the following manner:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [], "source": [ " # sensitivity analytical values\n", " Sa, Szw, Sta = analytic_sensitivity_coefficients(zm, wm)\n", " \n", "\n", " # 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)\n", "\n", " # 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)\n", "\n", " # 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))\n", "\n", " # 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))\n", " \n", " # end second order\n", " convergence_analysis = False\n", " if convergence_analysis:\n", " # Convergence analysis\n", " # 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]')\n", "\n", " # 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()\n", "\n", " # 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]')\n", "\n", " # 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()\n", "\n", " # # Convergence Monte Carlo with sobol 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 = 10\n", " # for i, N_smpl in enumerate(list_of_samples):\n", " # for j in range(N_iter):\n", " # A_s, XB, XC, Y_A, Y_B, Y_C, S, ST = monte_carlo_sens(N_smpl,\n", " # jpdf,\n", " # sample_method='S')\n", " # s_mc_err[i] += np.abs(S - Sa)\n", " # st_mc_err[i] += np.abs(ST - Sta)\n", " #\n", " # print('MC convergence analysis: N_smpl {} - finished with iteration {} of {}'.format(N_smpl, 1 + j, N_iter))\n", " # s_mc_err[i] /= float(N_iter)\n", " # st_mc_err[i] /= float(N_iter)\n", " #\n", " # fig_sobol = plt.figure('Sobol sampling')\n", " # fig_sobol.suptitle('Sobol 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()\n", " #\n", " # fig_random = plt.figure('Sobol sampling - average of indices')\n", " # fig_random.suptitle('Sobol 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]')\n", "\n", " plt.show()\n", " plt.close()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The actual algorithm calculating the sensitivity analysis was hidden in this function call which did the magic for us: `A_s, B_s, C_s, f_A, f_B, f_C, Smc, Stmc = monte_carlo_sens_nonlin(Ns_mc, jpdf)` \n", "\n", "Below we explain in greater detail Saltelli's algorithm which is used to compute the Sobol indices. \n", "\n", "### Saltelli's algorithm for Sobol indices estimation\n", "\n", "Calculating the sensitivity coefficients with MCM directly is\n", "computationally very expensive. To see this, consider how on would\n", "estimate $\\operatorname{Var}\\mathbb{E}(Y|Z_i))$ which is the numerator in the Sobol\n", "indices, in a direct brute force, manner. Let $M$ be the evaluations\n", "needed to estimate the inner, conditional expected value $\\mathbb{E}(Y|Z_i)$\n", "for a fixed $Z_i$. To get an approxiamation of the outer variance, one\n", "would have to repeat this process for the whole range of $Z_i$, which\n", "could also amount to $\\propto M$. Finally, this would have to be done\n", "for all $r$ input random variables of $Y$. Consecquently, the number\n", "of evalutations amounts to $\\mathcal{O}(M^2 \\;r)$. To get a impression\n", "of what this could to, note that in many cases a reasonable $M$ could\n", "be $5000$ which would results in $M^2 =25 000 000$ necessary\n", "evaluations!\n", "\n", "Luckily Saltelli came up with an algorithm to approximate of the sensitivity first order coefficients using $M(p+2)$ evaluations in total\n", "There are many adaptations and improvements of the algorithm available, here we will present the basic idea of the algorithm.\n", "\n", "\n", "\n", "
\n", "**Saltelli's algorithm.**\n", "\n", "1. Use a sampling method to draw a set of input samples $\\mathbf{z}^{(s)}$\n", "\n", "2. Evaluate the deterministic model $y(\\mathbf{z})$ for each sample\n", "\n", "3. Estimate all sensitivity indices from $y^{(s)}$.\n", "
\n", "\n", "\n", "\n", "Thus, the blackbox function mentioned above, follows these steps:" ] }, { "cell_type": "code", "execution_count": 7, "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" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Saltelli's algorithm step by step\n", "### Step 1: sample matrix creation\n", "\n", "For Saltellis Algorithm we need to create two different sample matrices $A,B$ each of the size $M\\times P$:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{align*}\n", "\\mathbf{A} =\n", "\\begin{bmatrix}\n", "z_1^{(A,1)} & \\cdots z_i^{(A,1)} \\cdots & z_P^{(A,1)} \\\\\n", "\\vdots &\t\t & \\vdots \\\\\n", "z_i^{(A,M)} & \\cdots z_i^{(A,M)} \\cdots & z_P^{(A,M)}\n", "\\end{bmatrix}\n", ", \\quad\n", "\\mathbf{B} =\n", "\\begin{bmatrix}\n", "z_1^{(B,1)} & \\cdots z_i^{(B,1)} \\cdots & z_P^{(B,1)} \\\\\n", "\\vdots &\t\t & \\vdots \\\\\n", "z_i^{(B,M)} & \\cdots z_i^{(B,M)} \\cdots & z_P^{(B,M)}\n", "\\end{bmatrix}\n", ".\n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In addition we create $P$ additional matrices $C_i$ of the size $M\\times P$ compound of matrix $A$ and matrix $B$. In a matrix $C_i$ all colums will be have the same values as the $B$ matrix, except the $i$-th column, which will have the values of $A$:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{align*}\n", "\\mathbf{C}_i =\n", "\\begin{bmatrix}\n", "z_1^{(B,1)} & \\cdots z_i^{(A,1)} \\cdots & z_P^{(B,1)} \\\\\n", "\\vdots &\t\t & \\vdots \\\\\n", "z_i^{(B,M)} & \\cdots z_i^{(A,M)} \\cdots & z_P^{(B,M)}\n", "\\end{bmatrix}\n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This was implemented in the method:\n", "`A, B, C = generate_sample_matrices_mc(number_of_samples, number_of_parameters, joint_distribution, sample_method)`." ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# sample matrices\n", "def generate_sample_matrices_mc(Ns, number_of_parameters, jpdf, sample_method='R'):\n", "\n", " Xtot = jpdf.sample(2*Ns, sample_method).transpose()\n", " A = Xtot[0:Ns, :]\n", " B = Xtot[Ns:, :]\n", "\n", " C = np.empty((number_of_parameters, Ns, number_of_parameters))\n", " # create C sample matrices\n", " for i in range(number_of_parameters):\n", " C[i, :, :] = B.copy()\n", " C[i, :, i] = A[:, i].copy()\n", "\n", " return A, B, C" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Step 2: evaluate the model for samples\n", "\n", "In the second step we evaluate the model for samples in the matrices\n", "and save the results in vectors $Y_{\\mathbf{A}}$, $Y_{\\mathbf{B}}$ and\n", "$Y_{\\mathbf{C_i}}$:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{align*}\n", "Y_{\\mathbf{A}} = y(\\mathbf{A}), \\qquad Y_{\\mathbf{B}} = y(\\mathbf{B}), \\qquad Y_{\\mathbf{C_i}} = y(\\mathbf{C_i}),\n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The corresponding python code for our example:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# 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": [ "### Step 3: approximate the sensitivity indices\n", "\n", "In the final step the first order and total Sobol indices are estimated.\n", "Since the numerical approximation of all indices are quite demanding, approximations are used to speed up the process.\n", "For both, the first and total sensitivity index, exist several approximations, which the most common can be found in ([[saltelli2010]](#saltelli2010)).\n", "\n", "### The first order sensitivity indices\n", "\n", "The first order indices are defined as:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{align*}\n", "S_i = \\frac{\\operatorname{Var}\\left(\\mathbb{E}(Y)| Z_i \\right)}{\\operatorname{Var}(Y)}\n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Both, the nominator and denominator are now approximated numerically, whereas the variance (nominator) is defined with:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{align*}\n", "\\operatorname{Var}(Y) = \\left(\\frac{1}{M-1} \\sum_{j=1}^M \\left(y_{\\mathbf{B}}^j\\right)^2\\right) - f_0^2\n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "with $f_0^2$ which is $\\left(\\mathbb{E}(Y)\\right)^2$.\n", "For $f_0^2$ exist several approximations, two common are:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{align*}\n", "f_0^2 = \\frac{1}{M^2} \\left(\\sum_{j=1}^M y_{\\mathbf{A}}^j \\right) \\left( \\sum_{j=1}^M y_{\\mathbf{B}}^j \\right)\n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The conditional variance is approximated as:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{align*}\n", "\\operatorname{Var}\\left(\\mathbb{E}(Y)| Z_i \\right) = \\frac{1}{M-1} \\sum_{j=1}^M y_{\\mathbf{A}}^j y_{\\mathbf{C_i}}^j - f_0^2\n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### The total indices" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{align*}\n", "S_{Ti} = \\frac{\\mathbb{E}\\left(\\operatorname{Var}(Y)| \\mathbf{Z}_{-i} \\right)}{\\operatorname{Var}(Y)} = 1 - \\frac{\\operatorname{Var}\\left(\\mathbb{E}(Y)| \\mathbf{Z}_{-i} \\right)}{\\operatorname{Var}(Y)}\n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here the variance is estimated accordingly, but taking the matrix A:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{align*}\n", "\\operatorname{Var}(Y) = \\left(\\frac{1}{M-1} \\sum_{j=1}^M \\left(y_{\\mathbf{A}}^j\\right)^2\\right) - f_0^2\n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "here $f_0^2$ is approximated with:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{align*}\n", "f_0^2 = \\frac{1}{M^2} \\left(\\sum_{j=1}^M y_{\\mathbf{A}}^j \\right)\n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And the conditional variance of not given $Z_i$ is approximated with:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{align*}\n", "\\operatorname{Var}\\left(\\mathbb{E}(Y)| \\mathbf{Z}_{-i} \\right) = \\left(\\frac{1}{M-1} \\sum_{j=1}^M y_{\\mathbf{B}}^j y_{\\mathbf{C_i}}^j\\right) - f_0^2\n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Those equations are implemented in the following way:" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# mc algorithm for variance based sensitivity coefficients\n", "def calculate_sensitivity_indices_mc(y_a, y_b, y_c):\n", "\n", " # single output value y_a for one set of samples\n", " if len(y_c.shape) == 2:\n", " Ns, n_parameters = y_c.shape\n", "\n", " # for the first order index\n", " f0sq_first = np.sum(y_a*y_b)/ Ns \n", " y_var_first = np.sum(y_b**2.)/(Ns-1) - f0sq_first\n", "\n", " # for the total index\n", " f0sq_total = (sum(y_a)/Ns)**2\n", " y_var_total = np.sum(y_a**2.)/(Ns-1) - f0sq_total\n", "\n", " s = np.zeros(n_parameters)\n", " st = np.zeros(n_parameters)\n", "\n", " for i in range(n_parameters):\n", " # first order index\n", " cond_var_X = np.sum(y_a*y_c[:, i])/(Ns - 1) - f0sq_first\n", " s[i] = cond_var_X/y_var_first\n", "\n", " # total index\n", " cond_exp_not_X = np.sum(y_b*y_c[:, i])/(Ns - 1) - f0sq_total\n", " st[i] = 1 - cond_exp_not_X/y_var_total\n", "\n", " # vector output value y_a,.. for one set of samples\n", " elif len(y_c.shape) == 3:\n", " n_y, Ns, n_parameters = y_c.shape\n", " # for the first order index\n", " f0sq_first = np.sum(y_a*y_b, axis=1) / Ns\n", " y_var_first = np.sum(y_b ** 2., axis=1) / (Ns - 1) - f0sq_first\n", "\n", " # for the total index\n", " f0sq_total = (np.sum(y_a, axis=1) / Ns) ** 2\n", " y_var_total = np.sum(y_a ** 2., axis=1) / (Ns - 1) - f0sq_total\n", "\n", " s = np.zeros((n_parameters, n_y))\n", " st = np.zeros((n_parameters, n_y))\n", "\n", " for i in range(n_parameters):\n", " # first order index\n", " cond_var_X = np.sum(y_a * y_c[:, :, i], axis=1) / (Ns - 1) - f0sq_first\n", "\n", " s[i, :] = cond_var_X / y_var_first\n", "\n", " # total index\n", " cond_exp_not_X = np.sum(y_b * y_c[:, :, i], axis=1) / (Ns - 1) - f0sq_total\n", " st[i, :] = 1 - cond_exp_not_X / y_var_total\n", "\n", " return s, st" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# References\n", "\n", " 1.
**A. Saltelli, P. Annoni, I. Azzini, F. Campolongo, M. Ratto and S. Tarantola**. \n", " Variance based sensitivity analysis of model output. Design and estimator for the total sensitivity index,\n", " *Computer Physics Communications*,\n", " 181(2),\n", " pp. 259-270,\n", " 2010." ] } ], "metadata": {}, "nbformat": 4, "nbformat_minor": 2 }