{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "# An introduction to the Sensitivity Analysis with the Monte Carlo method\n", "\n", " **Vinzenz Gregor Eck**, Expert Analytics\n", "\n", "Date: **Mar 9, 2017**\n", "hello\n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline\n", "%load_ext autoreload\n", "%autoreload 2\n", "import numpy as np\n", "import monte_carlo\n", "from sensitivity_examples_nonlinear import generate_distributions\n", "from sensitivity_examples_nonlinear import calculate_sensitivity_indices_non_additive_model" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Introduction\n", "Previously we conducted sensitivity analysis for the following model:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# start the linear model\n", "def linear_model(w, z):\n", " assert w.shape == z.shape\n", " return np.sum(w*z, axis=1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The analysis was done with the following lines of code:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": true }, "outputs": [], "source": [ " # sensitivity reference values\n", " S_book = np.array([0.0006, 0.009, 0.046, 0.145, 0, 0, 0, 0])\n", " St_book = np.array([0.003, 0.045, 0.229, 0.723, 0.002, 0.036, 0.183, 0.578])\n", "\n", " # get joint distributions\n", " joint_distribution = generate_distributions()\n", "\n", " number_of_samples = 1000000\n", " A, B, C, Y_A, Y_B, Y_C, S, ST = calculate_sensitivity_indices_non_additive_model(number_of_samples,\n", " joint_distribution)\n", "\n", " print(\"First Order Indices Total Sensitivity Indices\\n\")\n", " print(' S_est | S_ref ST_est | ST_ref\\n')\n", " for k, (s, sb, s_t, sb_t) in enumerate(zip(S, S_book, ST, St_book)):\n", " print('S_{} : {:>6.3f} | {:2.3f} ST_{} : {:>6.3f} | {:2.3f}'.format(k + 1, s, sb, k+1, s_t, sb_t))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The actual algorithm calculating the sensitivity analysis was hidden in this function call which did the magic for us:\n", "\n", "\"\"A, B, C, Y_A, Y_B, Y_C, S, ST = calculate_sensitivity_indices_non_additive_model(number_of_samples,\n", " joint_distribution)\"\"\n", "\n", "Now we want to look inside these function and look at the same time at the algorithm used.\n", "\n", "# Saltelli's Algorithm\n", "\n", "Brief declaration Introduction..\n", "\n", "The Algorithm can be split up in three major parts:\n", "\n", "1. Generate sample matrices\n", "2. Evaluate the model\n", "3. Approximate the sensitivity indices\n", "\n", "We implemented these steps as python functions to keep it close to the algorithm steps." ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# calculate sens indices of non additive model\n", "def calculate_sensitivity_indices_non_additive_model(number_of_samples, joint_distribution, sample_method='R'):\n", "\n", " number_of_parameters = len(joint_distribution)\n", "\n", " # 1. Generate sample matrices\n", " A, B, C = generate_sample_matrices_mc(number_of_samples, number_of_parameters, joint_distribution, sample_method)\n", "\n", " # 2. Evaluate the model\n", " Y_A, Y_B, Y_C = evaluate_non_additive_linear_model(A, B, C)\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" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Generate sample matrices" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# sample matrices\n", "def generate_sample_matrices_mc(number_of_samples, number_of_parameters, joint_distribution, sample_method='R'):\n", "\n", " Xtot = joint_distribution.sample(2 * number_of_samples, sample_method).transpose()\n", " A = Xtot[0:number_of_samples, :]\n", " B = Xtot[number_of_samples:, :]\n", "\n", " C = np.empty((number_of_parameters, number_of_samples, 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": [ "## Evaluate the model" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# model evaluation\n", "def evaluate_non_additive_linear_model(A, B, C):\n", "\n", " number_of_parameters = A.shape[1]\n", " number_of_samples = A.shape[0]\n", " number_of_factors = int(number_of_parameters / 2)\n", " # 1. evaluate sample matrices A\n", " Z_A = A[:, :number_of_factors]\n", " W_A = A[:, number_of_factors:]\n", " Y_A = linear_model(W_A, Z_A)\n", "\n", " # 2. evaluate sample matrices B\n", " Z_B = B[:, :number_of_factors]\n", " W_B = B[:, number_of_factors:]\n", " Y_B = linear_model(W_B, Z_B)\n", "\n", " # 3. evaluate sample matrices C\n", " Y_C = np.empty((number_of_samples, number_of_parameters))\n", " for i in range(number_of_parameters):\n", " x = C[i, :, :]\n", " z = x[:, :number_of_factors]\n", " w = x[:, number_of_factors:]\n", " Y_C[:, i] = linear_model(w, z)\n", "\n", " return Y_A, Y_B, Y_C" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Approximate the sensitivity indices" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# mc algorithm for variance based sensitivity coefficients\n", "def calculate_sensitivity_indices_mc(y_a, y_b, y_c):\n", "\n", " number_of_samples, n_parameters = y_c.shape\n", "\n", " # Shift the data around zero\n", " y_mean = np.mean(np.hstack([y_a, y_b]))\n", " y_a -= y_mean\n", " y_b -= y_mean\n", " y_c -= y_mean\n", "\n", " f0sq = np.mean(y_a * y_b)\n", " # f0sq = (sum(Y_A) / number_of_samples) ** 2\n", "\n", " y_a_var = np.sum(y_a ** 2.) / number_of_samples - f0sq\n", " y_b_var = np.sum(y_b ** 2.) / number_of_samples - f0sq\n", "\n", " s = np.zeros(n_parameters)\n", " st = np.zeros(n_parameters)\n", "\n", " for i in range(n_parameters):\n", " cond_var_X = np.sum(y_a * y_c[:, i]) / number_of_samples - f0sq\n", " s[i] = cond_var_X / y_b_var\n", "\n", " cond_var_not_X = np.sum(y_b * y_c[:, i]) / number_of_samples - f0sq\n", " st[i] = 1 - cond_var_not_X / y_a_var\n", "\n", " return s, st" ] } ], "metadata": { "anaconda-cloud": {}, "kernelspec": { "display_name": "Python [default]", "language": "python", "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.13" } }, "nbformat": 4, "nbformat_minor": 2 }