{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Sensitivity Analysis\n", "\n", "### [Neil D. Lawrence](http://inverseprobability.com), University of\n", "\n", "Cambridge\n", "\n", "### 2023-11-07" ], "id": "701fa31f-d310-49af-80b7-90a45602a445" }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Abstract**: This week we introduce sensitivity analysis through\n", "Emukit, showing how Emukit can deliver Sobol indices for understanding\n", "how the output of the system is affected by different inputs." ], "id": "3377b37a-a3ce-4893-a012-2cff14ba6e1e" }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "$$" ], "id": "e26a2108-f8f3-454d-b038-609911ab391c" }, { "cell_type": "markdown", "metadata": {}, "source": [ "::: {.cell .markdown}\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "" ], "id": "e48727e9-25c2-460d-92d6-d2c2b667dd0e" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Setup\n", "\n", "\\[edit\\]" ], "id": "13200db6-bf73-4437-9cd0-56aed8472005" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "plt.rcParams.update({'font.size': 22})" ], "id": "dd0276d0-f23f-46d0-af05-55a6642b72f2" }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ], "id": "42aa0041-1a04-4c67-8ded-ea82b723efd0" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## notutils\n", "\n", "\\[edit\\]\n", "\n", "This small package is a helper package for various notebook utilities\n", "used below.\n", "\n", "The software can be installed using" ], "id": "e41bb598-9b47-4dc6-8bf3-1675a76ba533" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%pip install notutils" ], "id": "86a365af-ca72-4a41-beb3-1ceae2c13641" }, { "cell_type": "markdown", "metadata": {}, "source": [ "from the command prompt where you can access your python installation.\n", "\n", "The code is also available on GitHub:\n", "\n", "\n", "Once `notutils` is installed, it can be imported in the usual manner." ], "id": "19d4a0bd-4956-4a85-afa3-3560faa9b472" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import notutils" ], "id": "f3e8bec1-d938-48c6-94a1-31f2d7c737e0" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## pods\n", "\n", "\\[edit\\]\n", "\n", "In Sheffield we created a suite of software tools for ‘Open Data\n", "Science’. Open data science is an approach to sharing code, models and\n", "data that should make it easier for companies, health professionals and\n", "scientists to gain access to data science techniques.\n", "\n", "You can also check this blog post on [Open Data\n", "Science](http://inverseprobability.com/2014/07/01/open-data-science).\n", "\n", "The software can be installed using" ], "id": "50b79cf0-f8b2-46e3-b734-3db2c9576dfd" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%pip install pods" ], "id": "ab526a7f-47a1-4cd7-9539-9a6c600b54a8" }, { "cell_type": "markdown", "metadata": {}, "source": [ "from the command prompt where you can access your python installation.\n", "\n", "The code is also available on GitHub: \n", "\n", "Once `pods` is installed, it can be imported in the usual manner." ], "id": "a7531e22-dc74-4fb4-8d51-01f2f0c32bda" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import pods" ], "id": "cdb6daef-bd8f-40b0-a0c8-e6cf0ef91767" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## mlai\n", "\n", "\\[edit\\]\n", "\n", "The `mlai` software is a suite of helper functions for teaching and\n", "demonstrating machine learning algorithms. It was first used in the\n", "Machine Learning and Adaptive Intelligence course in Sheffield in 2013.\n", "\n", "The software can be installed using" ], "id": "0d634efa-551b-4e48-9542-12144d84774a" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%pip install mlai" ], "id": "072ac4d2-429d-448f-83b1-ac8911c3af4d" }, { "cell_type": "markdown", "metadata": {}, "source": [ "from the command prompt where you can access your python installation.\n", "\n", "The code is also available on GitHub: \n", "\n", "Once `mlai` is installed, it can be imported in the usual manner." ], "id": "88c59751-9d23-4f4d-9619-30e13ecf1d8c" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import mlai" ], "id": "bcf4dfb1-3d37-4107-ab7d-5be7f87f60d2" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%pip install gpy" ], "id": "daf1dfcc-442e-4d1e-bd53-a4126a121fdc" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%pip install pyDOE" ], "id": "34644e58-64e2-46d7-bb62-b553a786fcb1" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%pip install emukit" ], "id": "be6d7efa-13a6-42e8-8eff-6f175626db8b" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Emukit Sensitivity Analysis\n", "\n", "\\[edit\\]\n", "\n", "This introduction is based on [Introduction to Global Sensitivity\n", "Analysis with\n", "Emukit](https://github.com/EmuKit/emukit/blob/master/notebooks/Emukit-tutorial-sensitivity-montecarlo.ipynb)\n", "written by Mark Pullin, Javier Gonzalez, Juan Emmanuel Johnson and\n", "Andrei Paleyes. Some references include (Kennedy and O’Hagan, 2000;\n", "Saltelli et al., 2010, 2008, 2004; Sobol, 2001, 1990)\n", "\n", "> A possible definition of sensitivity analysis is the following: The\n", "> study of how uncertainty in the output of a model (numerical or\n", "> otherwise) can be apportioned to different sources of uncertainty in\n", "> the model input (Saltelli et al., 2004). A related practice is\n", "> ‘uncertainty analysis’, which focuses rather on quantifying\n", "> uncertainty in model output. Ideally, uncertainty and sensitivity\n", "> analyses should be run in tandem, with uncertainty analysis preceding\n", "> in current practice.\n", ">\n", "> In Chapter 1 of Saltelli et al. (2008)" ], "id": "7b56d93e-4a82-42cb-8237-879e74ee31e1" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "from matplotlib import colors as mcolors\n", "from matplotlib import cm" ], "id": "c3131527-6a93-46b3-9fce-36df14b69ca7" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%pip install pyDOE" ], "id": "0a497bfc-982a-437b-aa8e-e9f2fd072684" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import mlai\n", "import mlai.plot as plot" ], "id": "d659734d-e97a-4f63-8dd2-97cdf6af1b36" }, { "cell_type": "markdown", "metadata": {}, "source": [ "Sensitivity analysis is a statistical technique widely used to test the\n", "reliability of real systems. Imagine a simulator of taxis picking up\n", "customers in a city like the one showed in the [Emukit\n", "playground](https://github.com/amzn/emukit-playground). The profit of\n", "the taxi company depends on factors like the number of taxis on the road\n", "and the price per trip. In this example, a global sensitivity analysis\n", "of the simulator could be useful to decompose the variance of the profit\n", "in a way that can be assigned to the input variables of the simulator.\n", "\n", "There are different ways of doing a sensitivity analysis of the\n", "variables of a simulator. In this notebook we will start with an\n", "approach based on Monte Carlo sampling that is useful when evaluating\n", "the simulator is cheap. If evaluating the simulator is expensive,\n", "emulators can then be used to speed up computations. We will show this\n", "in the last part of the notebook. Next, we start with a few formal\n", "definitions and literature review so we can understand the basics of\n", "Sensitivity Analysis and how it can be performed with Emukit." ], "id": "2e64cc84-3c29-429e-819d-24e777ec24e5" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Local Sensitivity\n", "\n", "Given any function, $g(\\cdot)$, we might be interested in how sensitive\n", "that function is to variations in its input space. One route to\n", "determining this is to compute the partial derivatives of that function\n", "with respect to its inputs, $$\n", "\\frac{\\partial}{\\partial x_i} g(\\mathbf{ x}).\n", "$$ The matrix of all these partial derivatives is known as the Jacobian.\n", "\n", "These types of local sensitivity analysis can be used for determining\n", "the effect of changing an input variable around an operating point. But\n", "they don’t give us an understanding of the response of the target\n", "function to variations in the input across the domain of inputs. For\n", "this, we need to look to *global sensitivity analysis*." ], "id": "135d505a-4e97-47ea-ac08-c77f6ed8f0c2" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Global Sensitivity Analysis\n", "\n", "In global sensitivity analysis, rather than looking around a single\n", "operating point, we’re interested in the overall sensitivity of a\n", "function to its inputs, or combinations of inputs, across its entire\n", "domain. The key tool in determining this sensitivity is known as the\n", "ANOVA decomposition, or the *Hoeffding-Sobol decomposition*.\n", "\n", "For global sensitivity analysis, we need to make an assumption about how\n", "inputs are going to vary to create different values of the function. The\n", "fundamental object we’re interested in is the total variance of the\n", "function, $$\n", "\\text{var}\\left(g(\\mathbf{ x})\\right) = \\left\\langle g(\\mathbf{ x})^2 \\right\\rangle _{p(\\mathbf{ x})} - \\left\\langle g(\\mathbf{ x}) \\right\\rangle _{p(\\mathbf{ x})}^2,\n", "$$ where $$\n", "\\left\\langle h(\\mathbf{ x}) \\right\\rangle _{p(\\mathbf{ x})} = \\int_\\mathbf{ x}h(\\mathbf{ x}) p(\\mathbf{ x}) \\text{d}\\mathbf{ x}\n", "$$ is the expectation of the function $h(\\mathbf{ x})$ under the density\n", "$p(\\mathbf{ x})$, which represents the probability distribution of\n", "inputs we’re interested in.\n", "\n", "The total variance of the function gives us the overall variation of the\n", "function across the domain of inputs, as represented by the probability\n", "density, $p(\\mathbf{ x})$. Normally, we perform analysis by assuming\n", "that, $$\n", "p(\\mathbf{ x}) = \\prod_{i=1}^pp(x_i)\n", "$$ and that each $p(x_i)$ is *uniformly distributed* across its input\n", "domain. Assuming we scale the input domain down to the interval\n", "$[0, 1]$, that gives us $$\n", "x_i \\sim \\mathcal{U}\\left(0,1\\right).\n", "$$" ], "id": "a198b13e-dcf0-443e-8721-4a761885c040" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Hoeffding-Sobol Decomposition\n", "\n", "The Hoeffding-Sobol, or ANOVA, decomposition of a function allows us to\n", "write it as, $$\n", "\\begin{align*}\n", "g(\\mathbf{ x}) = & g_0 + \\sum_{i=1}^pg_i(x_i) + \\sum_{i\\[edit\\]\n", "\n", "The Ishigami function (Ishigami and Homma, 1989) is a well-known test\n", "function for uncertainty and sensitivity analysis methods because of its\n", "strong nonlinearity and peculiar dependence on $x_3$. More details of\n", "this function can be found in (Sobol and Levitan, 1999).\n", "\n", "Mathematically, the form of the Ishigami function is $$\n", "g(\\textbf{x}) = \\sin(x_1) + a \\sin^2(x_2) + b x_3^4 \\sin(x_1). \n", "$$ We will set the parameters to be $a = 5$ and $b=0.1$ . The input\n", "variables are sampled randomly\n", "$x_i \\sim \\mathcal{U}\\left(-\\pi,\\pi\\right)$.\n", "\n", "Next, we create the function object and visualize its shape marginally\n", "for each one of its three inputs.\n", "\n", "Load the Ishigami function" ], "id": "c0f21145-97f8-400b-af00-f2b1ebd8eb94" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from emukit.test_functions.sensitivity import Ishigami" ], "id": "9baa9de0-0652-4a91-ac0d-0416fe954ffd" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ishigami = Ishigami(a=5, b=0.1)\n", "target_function = ishigami.fidelity1" ], "id": "d46b36f0-637c-4dd8-a289-fffc2232a6cb" }, { "cell_type": "markdown", "metadata": {}, "source": [ "That gives us the target function, next we define the input space for\n", "the simulator." ], "id": "2385663d-4f22-4062-960c-95635cabf2a9" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "from emukit.core import ContinuousParameter, ParameterSpace" ], "id": "a6fb119e-e5dc-485a-985a-7d0e20820e2a" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "variable_domain = (-np.pi,np.pi)\n", " \n", "space = ParameterSpace(\n", " [ContinuousParameter('x1', *variable_domain), \n", " ContinuousParameter('x2', *variable_domain),\n", " ContinuousParameter('x3', *variable_domain)])" ], "id": "738b77b6-ee50-49eb-b69c-268e83d34635" }, { "cell_type": "markdown", "metadata": {}, "source": [ "Before moving to any further analysis, we first plot the non-zero\n", "components $g(\\mathbf{ x})$. These components are $$\n", "\\begin{align*}\n", "g_1(x_1) & = \\sin(x_1) \\\\\n", "g_2(x_2) & = a \\sin^2 (x_2) \\\\\n", "g_{13}(x_1,x_3) & = b x_3^4 \\sin(x_1) \n", "\\end{align*}\n", "$$" ], "id": "59b480eb-6592-4e47-a8aa-aebdb330193c" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x_grid = np.linspace(*variable_domain,100)\n", "target_simulator = ishigami.fidelity1\n", "f1 = ishigami.f1(x_grid)\n", "f2 = ishigami.f2(x_grid)\n", "F13 = ishigami.f13(np.array([x_grid,x_grid]).T)[:,np.newaxis]" ], "id": "6e0929a3-9dfe-4b80-9467-d683bfae1ca1" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from mpl_toolkits.mplot3d import Axes3D" ], "id": "217f5f27-a939-47b4-97b4-612c1f612db6" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, axs = plt.subplots(2, 2, figsize=plot.big_wide_figsize)\n", "gs = axs[1, 1].get_gridspec()\n", "for ax in axs[1, 0:]:\n", " ax.remove()\n", "\n", "ax2 = fig.add_subplot(gs[1, 0:], projection='3d')\n", "\n", "axs[0,0].plot(x_grid, f1,'-r')\n", "axs[0,0].set_xlabel('$x_1$')\n", "axs[0,0].set_ylabel('$f_1$')\n", "\n", "axs[0,1].plot(x_grid,f2,'-r')\n", "axs[0,1].set_xlabel('$x_2$')\n", "axs[0,1].set_ylabel('$f_2$')\n", "\n", "X, Y = np.meshgrid(x_grid, x_grid)\n", "surf = ax2.plot_surface(X, Y, F13, cmap=cm.coolwarm, linewidth=0, antialiased=False)\n", "\n", "ax2.set_xlabel('$x_1$')\n", "ax2.set_ylabel('$x_3$')\n", "ax2.set_zlabel('$f_{13}$')\n", "\n", "mlai.write_figure(filename='non-zero-sobol-ishigami.svg', directory='./uq')" ], "id": "82c9066c-c201-4d14-a304-fb53aff3ab7f" }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "Figure: The non-zero components of the Ishigami function." ], "id": "bdb9f2c0-a0ee-429d-9496-10754e8aa999" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Total Variance\n", "\n", "The total variance $\\text{var}(y)$ in this example is" ], "id": "03694d5e-66e9-4b29-be39-8a3f33a652c7" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(ishigami.variance_total)" ], "id": "53e5b6ea-1d27-4749-9cf8-a6f0392303a3" }, { "cell_type": "markdown", "metadata": {}, "source": [ "which is the sum of the variance of $\\text{var}\\left(g_1(x_1)\\right)$,\n", "$\\text{var}\\left(g_2(x_2)\\right)$ and\n", "$\\text{var}\\left(g_{1,3}(x_{1,3})\\right)$" ], "id": "b8ee62c1-6d82-4154-8695-dc521dce188f" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(ishigami.variance_x1, ishigami.variance_x2, ishigami.variance_x13)\n", "print(ishigami.variance_x1 + ishigami.variance_x2 + ishigami.variance_x13)" ], "id": "15df70d6-cca4-446f-a4df-bfedae7e3ac0" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## First Order Sobol Indices using Monte Carlo\n", "\n", "The first order Sobol indices are a measure of “first order sensitivity”\n", "of each input variable. They account for the proportion of variance of\n", "$y$ explained by changing each variable alone while marginalizing over\n", "the rest. Recall that the Sobol index of the $i$th variable is computed\n", "as $$\n", "S_i = \\frac{\\text{var}\\left(g_i(x_i)\\right)}{\\text{var}\\left(g(\\mathbf{ x})\\right)}.\n", "$$ This value is standardized using the total variance, so it is\n", "possible to account for a fractional contribution of each variable to\n", "the total variance of the output.\n", "\n", "The Sobol indices for higher order interactions $S_{i,j}$ are computed\n", "similarly. Due to the normalization by the total variance, the the sum\n", "of all Sobol indices equals to one.\n", "\n", "In most cases we are interested in the first order indices. The Ishigami\n", "function has the benefit that these can be computed analytically. In\n", "`EmuKit` you can extract these values with the code." ], "id": "31672759-bf93-4980-9d8f-ab0effc9f557" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ishigami.main_effects" ], "id": "74c1d95b-2f6d-4219-af40-f2424416059b" }, { "cell_type": "markdown", "metadata": {}, "source": [ "But in general, these indices need to be sampled using Monte Carlo or\n", "one of the quasi-Monte Carlo methods we’ve seen in the model-free\n", "experimental design. Details are given in (Sobol, 2001).\n", "\n", "With Emukit, the first-order Sobol indices can be easily computed. We\n", "first need to define the space where the target simulator is analyzed." ], "id": "16caf190-2c7e-43f4-a333-1b5f1d6733a7" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from emukit.sensitivity.monte_carlo import ModelFreeMonteCarloSensitivity" ], "id": "f6e29d34-bfbb-42bb-862d-738cb28040a5" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "np.random.seed(10) # for reproducibility\n", "\n", "num_monte_carlo_points = 10000 # Number of MC samples\n", "senstivity_ishigami = ModelFreeMonteCarloSensitivity(target_simulator, space)\n", "main_effects, total_effects, _ = senstivity_ishigami.compute_effects(num_monte_carlo_points = num_monte_carlo_points)\n", "print(main_effects)" ], "id": "08d54dbd-aa6d-48df-bbe4-7bda91cd03f3" }, { "cell_type": "markdown", "metadata": {}, "source": [ "We compare the true effects with the Monte Carlo effects in a bar-plot.\n", "The total effects are discussed later." ], "id": "1ed849d6-8286-4f0c-b0f3-9576b3bc2eeb" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import pandas as pd" ], "id": "f81474b4-e466-473c-a1b2-e613a1a6da5f" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=plot.big_wide_figsize)\n", "\n", "d = {'Sobol True': ishigami.main_effects,\n", " 'Monte Carlo': main_effects}\n", "\n", "pd.DataFrame(d).plot(kind='bar', ax=ax)\n", "ax.set_title('First-order Sobol indices - Ishigami')\n", "ax.set_ylabel('% of explained output variance')\n", "\n", "mlai.write_figure(filename='first-order-sobol-indices-ishigami.svg', directory='./uq')" ], "id": "fe17692a-e326-492f-922b-829e4e66cfd3" }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "Figure: The non-zero components of the Ishigami function." ], "id": "d35b5bfc-0e14-485c-a6e9-7eadf6aea9f6" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Total Effects Using Monte Carlo\n", "\n", "Computing high order sensitivity indices can be computationally very\n", "demanding in high dimensional scenarios and measuring the total\n", "influence of each variable on the variance of the output is infeasible.\n", "To solve this issue the *total* indices are used which account for the\n", "contribution to the output variance of $x_i$ including all variance\n", "caused by the variable alone and all its interactions of any order.\n", "\n", "The total effect for $x_i$ is given by: $$ \n", "S_{Ti} = \\frac{\\left\\langle \\text{var}_{x_i} (y\\mid \\mathbf{ x}_{\\sim i}) \\right\\rangle _{p(\\mathbf{ x}_{\\sim i})}}{\\text{var}\\left(g(\\mathbf{ x})\\right)} = 1 - \\frac{\\text{var}_{\\mathbf{ x}_{\\sim i}} \\left\\langle y\\mid \\mathbf{ x}_{\\sim i} \\right\\rangle _{p(\\mathbf{ x}_{\\sim i})}}{\\text{var}\\left(g(\\mathbf{ x})\\right)}\n", "$$\n", "\n", "Note that the sum of $S_{Ti}$ is not necessarily one in this case unless\n", "the model is additive. In the Ishigami example the value of the total\n", "effects is" ], "id": "6b646ff8-43e5-46d5-916e-5a68232b8417" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ishigami.total_effects" ], "id": "41ac647c-cac9-4e27-81ea-79f5c522c331" }, { "cell_type": "markdown", "metadata": {}, "source": [ "As in the previous example, the total effects can be computed with Monte\n", "Carlo. In the next plot we show the comparison with the true total\n", "effects." ], "id": "e8b64438-da09-4011-822e-0d656a3e9e6c" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=plot.big_wide_figsize)\n", "\n", "d = {'Sobol True': ishigami.total_effects,\n", " 'Monte Carlo': total_effects}\n", "\n", "pd.DataFrame(d).plot(kind='bar', ax=ax)\n", "ax.set_title('Total effects - Ishigami')\n", "ax.set_ylabel('Effects value')\n", "\n", "mlai.write_figure(filename='total-effects-ishigami.svg', directory='./uq')" ], "id": "c9bacd63-5075-41cd-9a67-a14317c1a5d4" }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "Figure: The total effects from the Ishigami function as computed via\n", "Monte Carlo estimate alongside the true total effects for the Ishigami\n", "function." ], "id": "cc48fc62-c3d6-4ed7-ba5c-c835fdfb2a82" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Computing the Sensitivity Indices Using the Output of a Model\n", "\n", "In the example used above the Ishigami function is very cheap to\n", "evaluate. However, in most real scenarios the functions of interest are\n", "expensive, and we need to limit ourselves to a few number of\n", "evaluations. Using Monte Carlo methods is infeasible in these scenarios\n", "as a large number of samples are typically required to provide good\n", "estimates of the Sobol indices.\n", "\n", "An alternative in these cases is to use Gaussaian process emulator of\n", "the function of interest trained on a few inputs and outputs (Marrel et\n", "al., 2009). If the model is properly trained, its mean prediction which\n", "is cheap to evaluate, can be used to compute the Monte Carlo estimates\n", "of the Sobol indices, the variance from the GP emulator can also be used\n", "to assess our uncertainty about the Sobol indices. Let’s see how we can\n", "do this in Emukit.\n", "\n", "We start by generating 100 samples in the input domain. Note that this a\n", "just 1% of the number of samples that we used to compute the Sobol\n", "coefficients using Monte Carlo." ], "id": "1a0b4bf1-fc9b-4e8c-af52-4c6c7e47cc57" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from emukit.core.initial_designs import RandomDesign" ], "id": "1f650170-42fa-4055-bfea-6a21208d677e" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "design = RandomDesign(space)\n", "x = design.get_samples(500)\n", "y = ishigami.fidelity1(x)[:,np.newaxis]" ], "id": "c2a2b4d0-6f72-4514-a455-bfa9402a0466" }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, we fit a standard Gaussian process to the samples, and we wrap it\n", "as an Emukit model." ], "id": "57edc552-1675-4fd2-bb9b-dfd58a574a39" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from GPy.models import GPRegression\n", "from emukit.model_wrappers import GPyModelWrapper\n", "from emukit.sensitivity.monte_carlo import MonteCarloSensitivity" ], "id": "00f494a9-7fe7-48ec-8ea8-4ac6efa7d3f2" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "model_gpy = GPRegression(x,y)\n", "model_emukit = GPyModelWrapper(model_gpy)\n", "model_emukit.optimize()" ], "id": "bda6d5b8-d306-48fc-8870-e339d892f931" }, { "cell_type": "markdown", "metadata": {}, "source": [ "The final step is to compute the coefficients using the class\n", "`ModelBasedMonteCarloSensitivity` which directly calls the model and\n", "uses its predictive mean to compute the Monte Carlo estimates of the\n", "Sobol indices. We plot the true estimates, those computed using 10000\n", "direct evaluations of the object using Monte Carlo and those computed\n", "using a Gaussian process model trained on 100 evaluations." ], "id": "572e763b-3d5c-460e-af82-eb8b0ab855b8" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "num_mc = 10000\n", "senstivity_ishigami_gpbased = MonteCarloSensitivity(model = model_emukit, input_domain = space)\n", "main_effects_gp, total_effects_gp, _ = senstivity_ishigami_gpbased.compute_effects(num_monte_carlo_points = num_mc)" ], "id": "7a3ab658-f42f-4744-8b0d-33dcbfe3bec7" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=plot.big_wide_figsize)\n", "\n", "main_effects_gp = {ivar: main_effects_gp[ivar][0] for ivar in main_effects_gp}\n", "\n", "d = {'Sobol True': ishigami.main_effects,\n", " 'Monte Carlo': main_effects,\n", " 'GP Monte Carlo':main_effects_gp}\n", "\n", "pd.DataFrame(d).plot(kind='bar', ax=ax)\n", "plt.title('First-order Sobol indices - Ishigami')\n", "plt.ylabel('% of explained output variance')\n", "\n", "mlai.write_figure(filename='first-order-sobol-indices-gp-ishigami.svg', directory='./uq')" ], "id": "592a33ae-a2da-4b05-a6f1-effd3dfea244" }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "Figure: First order Sobol indices as estimated by Monte Carlo and\n", "GP-emulator based Monte Carlo." ], "id": "5b3d284a-5084-442b-b2c4-45f35253cf0d" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=plot.big_wide_figsize)\n", "\n", "total_effects_gp = {ivar: total_effects_gp[ivar][0] for ivar in total_effects_gp}\n", "\n", "d = {'Sobol True': ishigami.total_effects,\n", " 'Monte Carlo': total_effects,\n", " 'GP Monte Carlo':total_effects_gp}\n", "\n", "pd.DataFrame(d).plot(kind='bar', ax=ax)\n", "ax.set_title('Total effects - Ishigami')\n", "ax.set_ylabel('% of explained output variance')\n", "\n", "mlai.write_figure(filename='total-effects-sobol-indices-gp-ishigami.svg', directory='./uq')" ], "id": "c99a0f5a-4231-464d-ad87-4c2f60c1e1b3" }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "Figure: Total effects as estimated by Monte Carlo and GP based Monte\n", "Carlo.\n", "\n", "We observe some discrepancies with respect to the real value of the\n", "Sobol index when using the Gaussian process, but we get a fairly good\n", "approximation with a very reduced number of evaluations of the original\n", "target function." ], "id": "cccac528-bf33-432b-9071-7496de3e6cab" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Conclusions\n", "\n", "The Sobol indices are a tool for explaining the variance of the output\n", "of a function as components of the input variables. Monte Carlo is an\n", "approach for computing these indices if the function is cheap to\n", "evaluate. Other approaches are needed when $g(\\cdot)$ is expensive to\n", "compute." ], "id": "b2cdb6bf-ab5b-4b1e-9188-af25b425e92d" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Catapult Simulation\n", "\n", "\\[edit\\]\n", "\n", "\n", "\n", " \n", "\n", "\n", "\n", " \n", "\n", "\n", "\n", "Nicolas Durrande\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "As a worked example we’re going to introduce a catapult simulation\n", "written by Nicolas Durrande, .\n", "\n", "\n", "\n", "Figure: A catapult simulation for experimenting with surrogate\n", "models, kindly provided by Nicolas Durrande\n", "\n", "The simulator allows you to set various parameters of the catapult\n", "including the axis of rotation, `roation_axis`, the position of the arm\n", "stop, `arm_stop`, and the location of the two bindings of the catapult’s\n", "spring, `spring_binding_1` and `spring_binding_2`.\n", "\n", "These parameters are then collated in a vector, $$\n", "\\mathbf{ x}_i = \\begin{bmatrix}\n", "\\texttt{rotation_axis} \\\\\n", "\\texttt{arm_stop} \\\\\n", "\\texttt{spring_binding_1} \\\\\n", "\\texttt{spring_binding_2}\n", "\\end{bmatrix}\n", "$$\n", "\n", "Having set those parameters, you can run an experiment, by firing the\n", "catapult. This will show you how far it goes.\n", "\n", "Because you will need to operate the catapult yourself, we’ll create a\n", "function to query you about the result of an individual firing." ], "id": "e5448d31-4336-4ea0-9668-1cb58d62e36d" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np" ], "id": "fa429977-d4eb-40fc-b7fa-bdadcf37c1f4" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def catapult_distance(x):\n", " \"\"\"Request user input for the catapult.\"\"\"\n", " y = np.zeros((x.shape[0], 1))\n", " for i in range(x.shape[0]):\n", " rotation_axis=x[i, 0]\n", " arm_stop=x[i, 1]\n", " spring_binding_1=x[i, 2]\n", " spring_binding_2=x[i, 3]\n", " \n", " print('Please set the following values:')\n", " print('x_1 = {rotation_axis:.2f} (rotation axis)'.format(rotation_axis=rotation_axis))\n", " print('x_2 = {arm_stop:.2f} (arm stop)'.format(arm_stop=arm_stop))\n", " print('x_3 = {spring_binding_1:.2f} (spring binding 1)'.format(spring_binding_1=spring_binding_1))\n", " print('x_4 = {spring_binding_2:.2f} (spring binding 2)'.format(spring_binding_2=spring_binding_2))\n", " y[i, 0] = float(input('What is the distance? '))\n", " return y" ], "id": "f5c9331d-109c-43f7-888d-c69e57a22f60" }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can also set the parameter space for the model. Each of these\n", "variables is scaled to operate $\\in [0, 1]$." ], "id": "fc32f90e-a55d-4dc0-8136-fc240cc73dc4" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from emukit.core import ContinuousParameter, ParameterSpace" ], "id": "875e736e-7a82-4d88-941a-ca77b607036a" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "variable_domain = [0,1]\n", " \n", "space = ParameterSpace(\n", " [ContinuousParameter('rotation_axis', *variable_domain), \n", " ContinuousParameter('arm_stop', *variable_domain),\n", " ContinuousParameter('spring_binding_1', *variable_domain),\n", " ContinuousParameter('spring_binding_2', *variable_domain)])" ], "id": "b486ebfe-313c-4b88-93ac-e9f3e5655263" }, { "cell_type": "markdown", "metadata": {}, "source": [ "Before we perform sensitivity analysis, we need to build an emulator of\n", "the catapulter, which we do using our experimental design process." ], "id": "561646b6-f91a-4d6b-a19b-61500194ee26" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Experimental Design for the Catapult\n", "\n", "\\[edit\\]\n", "\n", "Now we will build an emulator for the catapult using the experimental\n", "design loop.\n", "\n", "We’ll start with a small model-free design, we’ll use a random design\n", "for initializing our model." ], "id": "d8a8fe61-052e-4435-9ad3-8458bc84b2c7" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from emukit.core.initial_designs import RandomDesign" ], "id": "75a3e8a4-b700-4711-b768-38db16aba4f1" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "design = RandomDesign(space)\n", "x = design.get_samples(5)\n", "y = catapult_distance(x)" ], "id": "f9390779-4c58-42e3-a2c8-08245ab655f0" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from GPy.models import GPRegression\n", "from emukit.model_wrappers import GPyModelWrapper\n", "from emukit.sensitivity.monte_carlo import MonteCarloSensitivity" ], "id": "62b93669-6604-4ef2-838b-cc02aaae03f9" }, { "cell_type": "markdown", "metadata": {}, "source": [ "Set up the GPy model. The variance of the RBF kernel is set to $150^2$\n", "because that’s roughly the square of the range of the catapult. We set\n", "the noise variance to a small value." ], "id": "832e5928-7ffc-4e33-9303-77dfa100dc6b" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "model_gpy = GPRegression(x,y)\n", "model_gpy.kern.variance = 150**2\n", "model_gpy.likelihood.variance.fix(1e-5)\n", "display(model_gpy)" ], "id": "fa3d6df6-d2cb-4723-91d0-1cff243600f9" }, { "cell_type": "markdown", "metadata": {}, "source": [ "Wrap the model for EmuKit." ], "id": "f401b436-0885-41e5-867a-c0aae2a59283" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "model_emukit = GPyModelWrapper(model_gpy)\n", "model_emukit.optimize()" ], "id": "1cc97693-e450-45c4-b1eb-e268e8ab0ba2" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "display(model_gpy)" ], "id": "96a671cb-cca8-49fa-b0ef-1fc8e610ba5d" }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we set up the model loop. We’ll use integrated variance reduction as\n", "the acquisition function for our model-based design loop.\n", "\n", "*Warning*: This loop runs much slower on Google `colab` than on a local\n", "machine." ], "id": "09ce8df0-1f60-40a7-ba6d-fd1a3056a003" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from emukit.experimental_design.experimental_design_loop import ExperimentalDesignLoop" ], "id": "4102f11e-0dda-465d-b279-46a5c7a6b7a3" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from emukit.experimental_design.acquisitions import IntegratedVarianceReduction, ModelVariance" ], "id": "f1ab785a-f8bc-434f-abf5-ee61743208f2" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "integrated_variance = IntegratedVarianceReduction(space=space,\n", " model=model_emukit)\n", "ed = ExperimentalDesignLoop(space=space, \n", " model=model_emukit, \n", " acquisition = integrated_variance)\n", "ed.run_loop(catapult_distance, 10)" ], "id": "7c4638e1-1a2a-49ce-9e8b-c46f2376e4d7" }, { "cell_type": "markdown", "metadata": {}, "source": [], "id": "4dda42bf-e1af-4db2-9665-f5375e9da511" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Sensitivity Analysis of a Catapult Simulation\n", "\n", "The final step is to compute the coefficients using the class\n", "`ModelBasedMonteCarloSensitivity` which directly calls the model and\n", "uses its predictive mean to compute the Monte Carlo estimates of the\n", "Sobol indices. We plot the estimates of the Sobol indices computed using\n", "a Gaussian process model trained on the observations we’ve acquired." ], "id": "cb308fe2-8ebe-4af7-a6e2-24ac4c414713" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "num_mc = 10000\n", "senstivity = MonteCarloSensitivity(model = model_emukit, input_domain = space)\n", "main_effects_gp, total_effects_gp, _ = senstivity.compute_effects(num_monte_carlo_points = num_mc)" ], "id": "b7cb0f2b-448c-42d3-a43c-992a0673af9a" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "import mlai.plot as plot\n", "import mlai" ], "id": "a3ab15b3-add7-4127-bd91-40f5a24459f4" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import pandas as pd" ], "id": "6324735b-481b-49fe-a451-177032f6b6e9" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=plot.big_wide_figsize)\n", "\n", "main_effects_gp_plot = {ivar: main_effects_gp[ivar][0] for ivar in main_effects_gp}\n", "\n", "d = {'GP Monte Carlo':main_effects_gp_plot}\n", "\n", "pd.DataFrame(d).plot(kind='bar', ax=ax)\n", "plt.ylabel('% of explained output variance')\n", "\n", "mlai.write_figure(filename='first-order-sobol-indices-gp-catapult.svg', directory='./uq')" ], "id": "781af768-a9d1-4f92-bd38-e1ad495a6f18" }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "Figure: First Order sobol indices as estimated by GP-emulator based\n", "Monte Carlo on the catapult." ], "id": "35eb2695-6cf1-44c4-898f-cc2ac8826576" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=plot.big_wide_figsize)\n", "\n", "total_effects_gp_plot = {ivar: total_effects_gp[ivar][0] for ivar in total_effects_gp}\n", "\n", "d = {'GP Monte Carlo':total_effects_gp_plot}\n", "\n", "pd.DataFrame(d).plot(kind='bar', ax=ax)\n", "ax.set_ylabel('% of explained output variance')\n", "\n", "mlai.write_figure(filename='total-effects-sobol-indices-gp-catapult.svg', directory='./uq')" ], "id": "fa2c57d5-ee49-4f80-aa83-825ff17d7615" }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "Figure: Total effects as estimated by GP based Monte Carlo on the\n", "catapult." ], "id": "ff6f9f67-f2a9-4b9e-8f58-7a56c8624767" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Thanks!\n", "\n", "For more information on these subjects and more you might want to check\n", "the following resources.\n", "\n", "- twitter: [@lawrennd](https://twitter.com/lawrennd)\n", "- podcast: [The Talking Machines](http://thetalkingmachines.com)\n", "- newspaper: [Guardian Profile\n", " Page](http://www.theguardian.com/profile/neil-lawrence)\n", "- blog:\n", " [http://inverseprobability.com](http://inverseprobability.com/blog.html)" ], "id": "ad607633-7760-4cd2-9f49-ef8b8a031b2e" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## References" ], "id": "5018c88d-1aef-4497-982c-8948a20cd092" }, { "cell_type": "markdown", "metadata": {}, "source": [ "Durrande, N., Ginsbourger, D., Olivier, Carraro, L., 2013. ANOVA kernels\n", "and RKHS of zero mean functions for model-based sensitivity analysis.\n", "Journal of Multivariate Analysis 115, 57–67.\n", "https://doi.org/\n", "\n", "Ishigami, T., Homma, T., 1989. An importance quantification technique in\n", "uncertainty analysis for computer models. \\[1990\\] Proceedings. First\n", "International Symposium on Uncertainty Modeling and Analysis 398–403.\n", "\n", "Kennedy, M.C., O’Hagan, A., 2000. [Predicting the output from a complex\n", "computer code when fast approximations are\n", "available](http://www.jstor.org/stable/2673557). Biometrika 87, 1–13.\n", "\n", "Marrel, A., Iooss, B., Laurent, B., Roustant, O., 2009. Calculations of\n", "Sobol indices for the Gaussian process metamodel. Reliability\n", "Engineering & System Safety 94, 742–751.\n", "https://doi.org/\n", "\n", "Saltelli, A., Annoni, P., Azzini, I., Campolongo, F., Ratto, M.,\n", "Tarantola, S., 2010. Variance based sensitivity analysis of model\n", "output. Design and estimator for the total sensitivity index. Computer\n", "Physics Communications 181, 259–270.\n", "\n", "\n", "Saltelli, A., Ratto, M., Andres, T., Campolongo, F., Cariboni, J.,\n", "Gatelli, D., Saisana, M., Tarantola, S., 2008. Global sensitivity\n", "analysis: The primer. wiley.\n", "\n", "Saltelli, A., Tarantola, S., Campolongo, F., Ratto, M., 2004.\n", "Sensitivity analysis in practice: A guide to assessing scientific\n", "methods. wiley.\n", "\n", "Sobol, I.M., 2001. Global sensitivity indices for nonlinear mathematical\n", "models and their Monte Carlo estimates. Mathematics and Computers in\n", "Simulation 55, 271–280. \n", "\n", "Sobol, I.M., 1990. On sensitivity estimation for nonlinear mathematical\n", "models. Matematicheskoe Modelirovanie 2, 112–118.\n", "\n", "Sobol, I.M., Levitan, Y.L., 1999. On the use of variance reducing\n", "multipliers in Monte Carlo computations of a global sensitivity index.\n", "Computer Physics Communications 117, 52–61.\n", "" ], "id": "04bf02e4-bc28-491b-8ba4-aa83d0e298ee" } ], "nbformat": 4, "nbformat_minor": 5, "metadata": {} }