{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "# A practical introduction to polynomial chaos with the chaospy package\n", "\n", " \n", "**Vinzenz Gregor Eck**, Expert Analytics, Oslo \n", "\n", " **Jacob Sturdy**, Department of Structural Engineering, NTNU\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 chaospy as cp\n", "import numpy as np\n", "from numpy import linalg as LA" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Introduction\n", "
\n", "\n", "To conduct UQSA with the polynomial chaos (PC) method we use the python package `chaospy`.\n", "This package includes many features and methods to apply non-intrusive polynomial chaos to any model with\n", "few lines of code.\n", "Since, `chaospy` allows for convenient definition of random variables, calculation of joint distributions and generation of samples,\n", "we apply `chaospy` as well for the Monte Carlo method.\n", "\n", "In the following we will briefly describe the theory of polynomial chaos and give a practical step by step guide\n", "for the application of `chaospy`.\n", "\n", "# Polynomial chaos\n", "
\n", "\n", "For a more in depth introduction to the polynomial chaos method see ([[xiu_numerical_2010;@ smith_uncertainty_2013]](#xiu_numerical_2010;@ smith_uncertainty_2013)).\n", "\n", "## Concept\n", "\n", "The main concept of polynomial chaos is to approximate the stochastic model response of a model with stochastic input as polynomial expansion relating the response to the inputs.\n", "For simplicity, we consider a model function $f$ which takes random variables $\\mathbf{Z}$ and non-random variables as\n", "spatial coordinates $x$ or time $t$, as input:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "Y = y{}(x, t, \\mathbf{Z})\n", "\\label{_auto1} \\tag{1}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then we seek to generate a polynomial expansion such that we can approximate the model response." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "Y = y{}(x, t, \\mathbf{Z}) \\approx y_{N} = \\sum_{i=1}^{P} v_i(x, t) \\Phi_i(\\mathbf{Z}),\n", "\\label{_auto2} \\tag{2}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "where $v_i$ are the expansion coefficients, which are only dependent on the non-stochastic parameters, and $\\Phi_i$ are orthogonal polynomials, which depend\n", "solely on the random input parameter $\\mathbf{Z}$.\n", "Once the polynomial expansion is calculated, we can analytically calculate the statistics of approximated model output.\n", "\n", "Since polynomial chaos generally requires fewer model evaluations to estimates statistics than Monte Carlo methods, this method is preferable applied for models which need long computational time.\n", "\n", "### Applicability\n", "\n", "Polynomial chaos may be applied in a wide variety of situations. Requirements for convergence are simply that the function be square integrable with respect to the inner product associated with the orthogonal polynomials. In otherwords if the function has a finite variance, polynomial chaos may be applied without worrying too much. The caveat to this fact is that convergence may not necessarily be fast! The smoother the function the faster the convergence. Polynomial choas can handle discontinuities, however, it may be advisable to reformulate the problem in these situations (see Feinberg, Eck and Langtangen 2016...).\n", "\n", "## Orthogonal polynomials\n", "As state above the polynomial chaos expansion consists of a sum of basis polynomials in the input parameters $\\mathbf{Z}$. \n", "By using orthogonal polynomials as the basis polynomials the efficiency of the convergence may be improved, and the use of the expansion in uncertainty quantification and sensitivity analysis is simplified.\n", "Orthogonality of functions is a general concept developed for functional analysis within an inner product space. Typically the inner product of two functions is defined as a weighted integral of the product of the functions over a domain of interest.\n", "\n", "\n", "\n", "
\n", "In particular for the purposes of polynomial chaos the inner product of two polynomials is defined as the expected value of their product, i.e. the integral of their product weighted with respect to the distribution of random input parameters.\n", "Two polynomials are orthogonal if their inner product is 0, and a set of orthogonal polynomials is orthogonal if every pair of distinct polynomials are orthogonal.\n", "\n", "The following equalities hold for the orthogonal basis polynomials used in polynomial chaos" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{align*}\n", " \\mathbb{E}(\\Phi_i(\\mathbf{Z})\\Phi_j(\\mathbf{Z})) &= \\langle \\Phi_i(\\mathbf{Z}),\\Phi_j(\\mathbf{Z}) \\rangle\\\\\n", " &= \\int_{\\Omega} \\Phi_i(z)\\Phi_j(z)w(z) dz \\\\\n", " &= \\int_{\\Omega} \\Phi_i(z)\\Phi_j(z)f(z) dz \\\\\n", " &= {\\int_{\\Omega} \\Phi_i(z)\\Phi_j(z) dF_Z(\\mathbf{z})} = h_j \\delta_{ij},\n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "where $h_j$ is equal to the normalisation factor of the used polynomials, and $\\Phi_i(\\mathbf{Z})$ indicates the substitution of the random variable $\\mathbf{Z}$ as the polynomial's variable.\n", "Note that $\\Phi_0$ is a polynomial of degree zero, i.e. a constant thus $\\mathbb{E}(\\Phi_0\\Phi_j) = \\Phi_0 \\mathbb{E}(\\Phi_j) = 0$ which implies that $\\mathbb{E}(\\Phi_j) = 0$ for $j \\neq 0$.\n", "
\n", "\n", "\n", "\n", "Once the random inputs $\\mathbf{Z}$ are properly defined with marginal distributions, orthogonal polynomials can be constructed.\n", "For most univariate distributions, polynomial bases functions are defined, and listed in the Asker Wilkeys scheme.\n", "A set of orthogonal polynomials can be created from those basis polynomials with orthogonalization methods as Cholesky decomposition, three terms recursion, and modified Gram-Schmidt.\n", "\n", "\n", "## Expansion coefficients\n", "\n", "There are two non-intrusive ways of approximating a polynomial chaos expansion coefficients:\n", "\n", "### Regression\n", "\n", "Supposing a polynomial expansion approximation $y_{N} = \\sum_i v_i \\Phi_i$,\n", "then *stochastic collocation* specifies a set of nodes, $\\Theta_{M} = \\left\\{\\mathbf{z}^{(s)}\\right\\}_{s=1}^{P}$, where the deterministic values of $y_{N}=y$.\n", "The task is thus to find suitable coefficients $v_i$ such that this condition is satisfied.\n", "Considering the existence of a set of collocation nodes $\\left\\{\\mathbf{z}^{(s)}\\right\\}_{s=1}^{P}$, then $y_{N} = \\sum_i v_i \\Phi_i$ can be\n", " formed into linear system of equations for the coefficients $v_i$ at these nodes:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "\\begin{bmatrix}\n", "\\Phi_0(\\mathbf{z}^{(1)}) & \\cdots & \\Phi_P(\\mathbf{z}^{(1)}) \\\\\n", "\\vdots & & \\vdots \\\\\n", "\\Phi_0(\\mathbf{z}^{(N)}) & \\cdot & \\Phi_P(\\mathbf{z}^{(n)})\n", "\\end{bmatrix}\n", "\\begin{bmatrix}\n", "v_{0}\\\\\n", "\\vdots \\\\\n", "v_{P}\n", "\\end{bmatrix}\n", "= \\begin{bmatrix}\n", "y(\\mathbf{z}^{(1)}) \\\\\n", "\\vdots \\\\\n", "y (\\mathbf{z}^{(N)})\n", "\\end{bmatrix}\n", "\\end{equation}\n", "\\label{eq:stochColl} \\tag{3}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we can use regression to achieve the relaxed condition that $y_{N}$ is \"sufficiently close\" to $y$ at $\\left\\{\\mathbf{z}^{(s)}\\right\\}_{s=1}^{P}$.\n", "This is done by choosing a larger number of samples so that ([3](#eq:stochColl)) is over determined and\n", "then minimizing the appropriate error ${\\lvert\\lvert {y_{N} \\rvert\\rvert}-y}_{R}$ over $\\left\\{\\mathbf{z}^{(s)}\\right\\}_{s=1}^{P}$.\n", "Ordinary least squares, Ridge-Regression, and Tikhonov-Regularization are all regression methods that may be applied to this problem.\n", "\n", "### Pseudo-spectral projection\n", "\n", "Discrete projection refers to the approximation of coefficients for $y_{N} = \\sum_{i=1}^{P} v_i \\Phi_i$, by directly approximating the orthogonal projection coefficients" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "v_i = \\frac{1}{h_i} \\mathbb{E}(y \\Phi_i) = \\frac{1}{h_i} \\int_{\\Omega} y \\Phi_i dF_z,\n", "\\label{_auto3} \\tag{4}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "using a quadrature scheme to calculate the integral $\\int_{\\Omega} y \\Phi_i dF_z$ as a sum $\\sum_{i=1}^{P} w_i y(\\mathbf{z}^{(i)})$, where $w_i$ are the quadrature weights.\n", "\n", "This results in an approximation $\\tilde{v}_{i}$ of $v_{i}$ the error between the final approximation may be split" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "\\lVert{y_{N} - y}\\rVert \\leq \\lVert \\sum_{i=1}^{P} v_i \\Phi_i - y{} \\rVert + \\lVert {\\sum_{i=1}^{P} \\left( v_i - \\tilde{v}_i \\right)\\Phi_i} \\rVert\n", "\\label{eq:ps_error} \\tag{5}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "where the first term is called the truncation error and the second term the quadrature error.\n", "Thus one may consider the maximal accuracy for a given polynomial order $P$ and see this should be achieved as the quadrature error is reduced to almost 0 by increased number of collocation nodes.\n", "\n", "## Statistics\n", "\n", "Once the expansion was generated, it can be used directly to calculated statistics for the uncertainty and sensitivity analysis.\n", "The two most common measures for uncertainty quantification, **expected value** and **variance**, can be calculated by\n", "inserting the expansion into the definition of the measures.\n", "\n", "The expected value is equal to the first expansion coefficient:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", " \\begin{aligned}\n", " {\\mathbb{E}}(Y) \\approx \\int_{\\Omega} \\sum_{i=1}^{P} v_i \\Phi_i(\\mathbf{z}{}) dF_Z(\\mathbf{z}) = v_1.\n", " \\end{aligned}\n", " \\label{eq:expectedValue_gPCE} \\tag{6}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The variance is the sum of squared expansion coefficients multiplied by normalisation constants of the polynomials:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", " \\begin{aligned}\n", " \\operatorname{Var}(Y) &\\approx {\\mathbb{E}}{({v(x,t,\\mathbf{Z})}(\\mathbf{Z})-{\\mathbb{E}}(Y))^2} = \\int_{\\Omega}({v(x,t,\\mathbf{Z})}(\\mathbf{z}) - v_1)^2 dF_Z(\\mathbf{z}) \\\\\n", " &= \\int_{\\Omega}\\left(\\sum_{i=1}^{P} v_i \\Phi_i(\\mathbf{z}) \\right)^2 dF_Z(\\mathbf{z}) - v_1^2 = \\sum_{i=1}^{P} v_i^2 \\int_{\\Omega}{\\Phi^2_i(\\mathbf{z})}dF_Z(\\mathbf{z}) - v_1^2 \\\\\n", " &= \\sum_{i=1}^{P} v_i^2 h_i - v_1^2 = \\sum_{i=2}^{P} v_i^2 h_i\n", " \\end{aligned}\n", " \\label{eq:variance_gPCE} \\tag{7}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "(Note the orthogonality of individual terms implies their covariance is zero, thus the variance is simply the sum of the variances of the terms.)\n", "\n", "The Sobol indices may be calculated quite simply from the expansion terms due to the fact that the ANOVA decomposition is unique. Thus the main effect of a parameter $z_i$ is simply the variance of all terms only in $z_i$. Explicitly let $\\mathcal{A}_{i} = \\{k | \\Phi_{k}(\\mathbf{z}) = \\Phi_{k}(z_i)\\} $ ,i.e. $\\mathcal{A}_{i}$ is the set of all indices of basis functions depending only on $z_i$ then" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", " \\begin{aligned}\n", " f_i &= \\sum_{k\\in \\mathcal{A}_{i}} v_{k} \\Phi_{k} \\implies \\\\\n", " S_i &= \\frac{1}{\\operatorname{Var}(Y)} \\sum_{k\\in \\mathcal{A}_{i}} \\operatorname{Var}(v_{k} \\Phi_{k}) \\\\\n", " S_i &= \\frac{1}{\\operatorname{Var}(Y)} \\sum_{k\\in \\mathcal{A}_{i}} v_{k}^2 h_{k} \n", " \\end{aligned}\n", " \\label{eq:sensitivity_gPCE} \\tag{8}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "and similarly one may define $\\mathcal{A}_{ij}$ for pairwise combinations of inputs and further to calculate all orders of sensitivities.\n", "\n", "\n", "\n", "\n", "# Chaospy\n", "
\n", "\n", "The python package `chaospy` an introductory paper to the package including a comparison to other software packages is\n", "presented here ([[feinberg_2015]](#feinberg_2015)).\n", "You can find an introduction, tutorials and the source code at the projects homepage: \n", "\n", "The installation of the package can be done via `pip`:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " pip install chaospy\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the following we will use the import naming convention of the package creator to import the package in python:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " import chaospy as cp\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Therefore it will be convenient to see whenever a method of the package is applied.\n", "\n", "The package `chaospy` is doc-string annotated which means that every method provides a short help text with small examples.\n", "To show the method documentation simply type a `?` after the method name in a ipython console or notebook.\n", "As shown in the following two examples:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# show help for uniform distributions\n", "cp.Uniform?" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# show help for sample generation\n", "cp.samplegen?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Steps for polynomial chaos analysis with chaospy\n", "
\n", "\n", "To conduct UQSA analysis with polynomial chaos we need to follow the following steps:\n", "\n", "* Definition of the marginal and joint distributions\n", "\n", "* Generation of the orthogonal polynomials\n", "\n", "* Linear regression\n", "\n", " * Generation of samples\n", "\n", " * `Evaluation of the model for all samples`\n", "\n", " * Generation of the polynomial chaos expansion\n", "\n", "\n", "* Pseudo-spectral projection\n", "\n", " * Generation of integration nodes and weights\n", "\n", " * `Evaluation of the model for all nodes`\n", "\n", " * Generation of the polynomial chaos expansion\n", "\n", "\n", "* Calculations of all statistics\n", "\n", "Note, that steps **3 Linear regression** and **4 Pseudo-spectral projection** are interchangeable. They are simply different methods of cacluating the expansion coefficients. In both cases\n", "generate a set of points in the parameter space where the model must be evaluated (steps 3.b and 4.b, respectively).\n", "\n", "\n", "## Step 1: Definition of marginal and joint distributions\n", "
\n", "\n", "The analysis of a each model starts with the definition of the marginal distributions for each random model input, i.e. describing it as\n", "random variable.\n", "Univariate random variables can be defined with `chaospy` by calling the class-constructor of a distribution type, e.g `cp.Normal()`, with arguments to describe the particular distribution, e.g. mean value and standard deviation for `cp.Normal`.\n", "The help function can be used to find out more about the required arguments, e.g. `help(cp.Normal)`.\n", "\n", "In the following an example for 3 random variables with uniform, normal and log-normal distribution:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# simple distributions\n", "rv1 = cp.Uniform(0, 1)\n", "rv2 = cp.Normal(0, 1)\n", "rv3 = cp.Lognormal(0, 1, 0.2, 0.8)\n", "print(rv1, rv2, rv3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "After all random input variables are defined with univariate random variables a multi-variate random variable and its joint distribution\n", "can be constructed with the following command:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# joint distributions\n", "joint_distribution = cp.J(rv1, rv2, rv3)\n", "print(joint_distribution)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It is also possible to construct independent identical distributed random variables from any univariate variable:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# creating iid variables\n", "X = cp.Normal()\n", "Y = cp.Iid(X, 4)\n", "print(Y)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "All 64 distributions available in the chaospy package can be found in the following table:\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
Distributions implemented in chaospy
Alpha Birnbaum-Sanders Laplace Power log-normal
Anglit Fisher-Snedecor Levy Power normal
Arcsinus Fisk/log-logistic Log-gamma Raised cosine
Beta Folded Cauchy Log-laplace Rayleigh
Brandford Folded normal Log-normal Reciprocal
Burr Frechet Log-uniform Right-skewed Gumbel
Cauchy Gamma Logistic Student-T
Chi Gen. exponential Lomax Triangle
Chi-square Gen. extreme value Maxwell Truncated exponential
Double Gamma Gen. gamma Mielke's beta-kappa Truncated normal
Double Weibull Gen. half-logistic Nakagami Tukey-Lambda
Epanechnikov Gilbrat Non-central chi-squared Uniform
Erlang Truncated Gumbel Non-central Student-T Wald
Exponential Gumbel Non-central F Weibull
Exponential power Hypergeometric secant Normal Wigner
Exponential Weibull Kumaraswamy Pareto (first kind) Wrapped Cauchy
\n", "\n", "## Step 2: Orthogonal Polynomials\n", "
\n", "\n", "The orthogonal polynomials can be generated with different methods, in `chaospy` there are 4 methods implemented. The\n", "most stable method, and therefore most advised is the *three terms recursion* method.\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
Orthogonalization Method
Cholesky decomposition cp.orth\\_chol
Three terms recursion cp.orth\\_ttr
Modified Gram-Schmidt cp.orth\\_gs
\n", "\n", "Regarding the *three terms recursion* method:\n", "For the distributions Normal, Uniform, Gamma, Log-normal, Triangle, Beta and stochastic independent variable combinations of those,\n", "the three terms recursion coefficients are known.\n", "For all other distributions the coefficients are estimated numerically.\n", "The *three terms recursion* method is then also called **discretized stieltjes method**.\n", "\n", "The most stable method and therefore most applied method is the **three terms recursion** (**discretized stieltjes method**) method.\n", "\n", "We will look at all in a small example, try to increase the polynomial order and the instabilities of the methods become visible." ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# example orthogonalization schemes\n", "# a normal random variable\n", "n = cp.Normal(0, 1)\n", "\n", "x = np.linspace(0,1, 50)\n", "# the polynomial order of the orthogonal polynomials\n", "polynomial_order = 3\n", "\n", "poly = cp.orth_chol(polynomial_order, n, normed=True)\n", "print('Cholesky decomposition {}'.format(poly))\n", "ax = plt.subplot(131)\n", "ax.set_title('Cholesky decomposition')\n", "_=plt.plot(x, poly(x).T)\n", "_=plt.xticks([])\n", "\n", "poly = cp.orth_ttr(polynomial_order, n, normed=True)\n", "print('Discretized Stieltjes / Three terms reccursion {}'.format(poly))\n", "ax = plt.subplot(132)\n", "ax.set_title('Discretized Stieltjes ')\n", "_=plt.plot(x, poly(x).T)\n", "\n", "poly = cp.orth_gs(polynomial_order, n, normed=True)\n", "print('Modified Gram-Schmidt {}'.format(poly))\n", "ax = plt.subplot(133)\n", "ax.set_title('Modified Gram-Schmidt')\n", "_=plt.plot(x, poly(x).T)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 3.: Linear regression\n", "\n", "The linear regression method requires to conduct the three following steps:\n", "\n", "1. Generation of samples\n", "\n", "2. `Evaluation of the model for all samples`\n", "\n", "3. Generation of the polynomial chaos expansion\n", "\n", "In the following we will not consider the model evaluation.\n", "\n", "### Step 3.a: Sampling\n", "\n", "
\n", "\n", "Once a random variable is defined or a joint random variable, also referred as distribution here, the following\n", "method can be used to generate as set of samples:" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# sampling in chaospy\n", "u = cp.Uniform(0,1)\n", "u.sample?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The method takes the arguments **size** which is the number of samples and **rule** which is the applied sampling scheme.\n", "The following example shows the creation of 2 set of samples for the sampling schemes *(Pseudo-)Random* and *Hammersley*." ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# example sampling\n", "u1 = cp.Uniform(0,1)\n", "u2 = cp.Uniform(0,1)\n", "joint_distribution = cp.J(u1, u2)\n", "number_of_samples = 350\n", "samples_random = joint_distribution.sample(size=number_of_samples, rule='R')\n", "samples_hammersley = joint_distribution.sample(size=number_of_samples, rule='M')\n", "\n", "fig1, ax1 = plt.subplots()\n", "ax1.set_title('Random')\n", "ax1.scatter(*samples_random)\n", "ax1.set_xlabel(\"Uniform 1\")\n", "ax1.set_ylabel(\"Uniform 2\")\n", "ax1.axis('equal')\n", "\n", "fig2, ax2 = plt.subplots()\n", "ax2.set_title('Hammersley sampling')\n", "ax2.scatter(*samples_hammersley)\n", "ax2.set_xlabel(\"Uniform 1\")\n", "ax2.set_ylabel(\"Uniform 2\")\n", "ax2.axis('equal')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "All sampling schemes implemented in chaospy are listed in the following table:\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
Key Name Nested
C Chebyshev nodes no
NC Nested Chebyshev yes
K Korobov no
R (Pseudo-)Random no
RG Regular grid no
NG Nested grid yes
L Latin hypercube no
S Sobol yes
H Halton yes
M Hammersley yes
\n", "### Step 3.a.i: Importing and exporting samples\n", "\n", "It may be useful to export the samples from `chaospy` for use in another program. The most useful format for exporting the samples likely depends on the external program, but it is quite simple to save the samples as a `CSV` file with a delimeter of your choice:" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# example save samples to file\n", "# Creates a csv file where each row corresponds to the sample number and each column with teh variables in the joint distribution\n", "csv_file = \"csv_samples.csv\"\n", "sep = '\\t'\n", "header = [\"u1\", \"u2\"]\n", "header = sep.join(header)\n", "np.savetxt(csv_file, samples_random, delimiter=sep, header=header)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Each row of the csv file now contains a single sample from the joint distribution with the columns corresponding to each component.\n", "\n", "Now you may evaluate these samples with an external program and save the resulting data into a similarly formatted `CSV` file. Again each row should correspond to a single sample value and each column to different components of the model output." ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# example load samples from file\n", "# loads a csv file where the samples/or model evaluations for each sample are saved\n", "# with one sample per row. Multiple components ofoutput can be stored as separate columns \n", "filepath = \"external_evaluations.csv\"\n", "data = np.loadtxt(filepath)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Step 3.c: Polynomial Chaos Expansion\n", "\n", "After the model is evaluated for all samples, the polynomial chaos expansion can be generated with the following method:" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# linear regression in chaospy\n", "cp.fit_regression?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the following we show a complete example for polynomial chaos expansion using the linear regression.\n", "The model applied the very simple mathematical expression:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", " y(z_1, z_2) = z_1 + z_1 z_2\n", "\\label{eq:dummy_model} \\tag{9}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The random variables for $Z_1, Z_2$ are defined as simple uniform random variables:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "Z_1 = \\mbox{U}(0,1), \\quad Z_2 = \\mbox{U}(0,1)\n", "\\label{eq:dummy_rv} \\tag{10}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The mean of this should be $\\frac{3}{4}$, the variance should be $\\frac{31}{144}$ and the sensitivites to $Z_1$ and $Z_2$ are respectively $\\frac{3}{31}$ and $\\frac{27}{31}$.\n", "\n", "Here is the annotated example code with all steps required to generate a polynomial chaos expansion with\n", "linear regression:" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# example linear regression\n", "# 1. define marginal and joint distributions\n", "u1 = cp.Uniform(0,1)\n", "u2 = cp.Uniform(0,1)\n", "joint_distribution = cp.J(u1, u2)\n", "\n", "# 2. generate orthogonal polynomials\n", "polynomial_order = 3\n", "poly = cp.orth_ttr(polynomial_order, joint_distribution)\n", "\n", "# 3.1 generate samples\n", "number_of_samples = cp.bertran.terms(polynomial_order, len(joint_distribution))\n", "samples = joint_distribution.sample(size=number_of_samples, rule='R')\n", "\n", "# 3.2 evaluate the simple model for all samples\n", "model_evaluations = samples[0]+samples[1]*samples[0]\n", "\n", "# 3.3 use regression to generate the polynomial chaos expansion\n", "gpce_regression = cp.fit_regression(poly, samples, model_evaluations)\n", "print(\"Success\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 4: Pseudo-spectral projection\n", "\n", "### Step 4.a: Quadrature nodes and weights\n", "\n", "
\n", "\n", "Once a random variable is defined or joint random variables, also referred as distribution here, the following\n", "method can be used to generate nodes and weights for different quadrature methods:" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# quadrature in polychaos\n", "cp.generate_quadrature?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will look at the following arguments of the method: **order** is the order of the quadrature, **domain** is\n", "the , **rule** is the *name* or *key* of the quadrature rule to apply.\n", "\n", "In the following example we look at some quadrature nodes for the same uniform variables as for the sampling,\n", "for Optimal Gaussian quadrature and Clenshaw-Curtis quadrature." ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# example quadrature\n", "u1 = cp.Uniform(0,1)\n", "u2 = cp.Uniform(0,1)\n", "joint_distribution = cp.J(u1, u2)\n", "\n", "order = 5\n", "\n", "nodes_gaussian, weights_gaussian = cp.generate_quadrature(order=order, domain=joint_distribution, rule='G')\n", "nodes_clenshaw, weights_clenshaw = cp.generate_quadrature(order=order, domain=joint_distribution, rule='C')\n", "\n", "print('Number of nodes gaussian quadrature: {}'.format(len(nodes_gaussian[0])))\n", "print('Number of nodes clenshaw-curtis quadrature: {}'.format(len(nodes_clenshaw[1])))\n", "\n", "\n", "fig1, ax1 = plt.subplots()\n", "ax1.scatter(*nodes_gaussian, marker='o', color='b')\n", "ax1.scatter(*nodes_clenshaw, marker= 'x', color='r')\n", "ax1.set_xlabel(\"Uniform 1\")\n", "ax1.set_ylabel(\"Uniform 2\")\n", "ax1.axis('equal')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the following all quadrature rules implemented in chaospy are highlighted:\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
Collection of quadrature rules Name Key
Optimal Gaussian quadrature Gaussian G
Gauss-Legendre quadrature Legendre E
Clenshaw-Curtis quadrature Clenshaw C
Leja quadrature Leja J
Hermite Genz-Keizter 16 rule Genz Z
Gauss-Patterson quadrature rule Patterson P
\n", "It is also possible to use sparse grid quadrature. For this purpose Clenshaw-Curtis method is advised since it is nested.\n", "\n", "In the following example we show sparse vs. normal quadrature nodes:" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# example sparse grid quadrature\n", "u1 = cp.Uniform(0,1)\n", "u2 = cp.Uniform(0,1)\n", "joint_distribution = cp.J(u1, u2)\n", "\n", "order = 2\n", "# sparse grid has exponential growth, thus a smaller order results in more points\n", "nodes_clenshaw, weights_clenshaw = cp.generate_quadrature(order=order, domain=joint_distribution, rule='C')\n", "nodes_clenshaw_sparse, weights_clenshaw_sparse = cp.generate_quadrature(order=order, domain=joint_distribution, rule='C', sparse=True)\n", "\n", "print('Number of nodes normal clenshaw-curtis quadrature: {}'.format(len(nodes_clenshaw[0])))\n", "print('Number of nodes clenshaw-curtis quadrature with sparse grid : {}'.format(len(nodes_clenshaw_sparse[0])))\n", "\n", "fig1, ax1 = plt.subplots()\n", "ax1.scatter(*nodes_clenshaw, marker= 'x', color='r')\n", "ax1.scatter(*nodes_clenshaw_sparse, marker= 'o', color='b')\n", "ax1.set_xlabel(\"Uniform 1\")\n", "ax1.set_ylabel(\"Uniform 2\")\n", "ax1.axis('equal')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Step 4.c: Polynomial Chaos Expansion\n", "\n", "After the model is evaluated for all integration nodes, the polynomial chaos expansion can be generated with the following method:" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# spectral projection in chaospy\n", "cp.fit_quadrature?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the following we show again a complete example for polynomial chaos expansion using the pseudo spectral approach to calculate the\n", "expansion coefficients.\n", "The model applied the same simple mathematical expression as before:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", " y(z_1, z_2) = z_1 + z_1 z_2\n", "\\label{eq:dummy_model_repeat} \\tag{11}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The random variables for $Z_1, Z_2$ are defined as simple uniform random variables:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "Z_1 = \\mbox{U}(0,1), \\quad Z_2 = \\mbox{U}(0,1)\n", "\\label{eq:dummy_rv_repeat} \\tag{12}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# example spectral projection\n", "# 1. define marginal and joint distributions\n", "u1 = cp.Uniform(0,1)\n", "u2 = cp.Uniform(0,1)\n", "joint_distribution = cp.J(u1, u2)\n", "\n", "# 2. generate orthogonal polynomials\n", "polynomial_order = 3\n", "poly = cp.orth_ttr(polynomial_order, joint_distribution)\n", "\n", "# 4.1 generate quadrature nodes and weights\n", "order = 5\n", "nodes, weights = cp.generate_quadrature(order=order, domain=joint_distribution, rule='G')\n", "\n", "# 4.2 evaluate the simple model for all nodes\n", "model_evaluations = nodes[0]+nodes[1]*nodes[0]\n", "\n", "# 4.3 use quadrature to generate the polynomial chaos expansion\n", "gpce_quadrature = cp.fit_quadrature(poly, nodes, weights, model_evaluations)\n", "print(\"Success\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 5: Statistical Analysis\n", "
\n", "\n", "Once the polynomial chaos expansion is created either with **pseudo-spectral projection** or with **regression** method\n", "The calculation of statistics is straight forward.\n", "The following listing gives an overview of all available methods take all the same input parameter the **polynomial-expansion** and\n", "the **joint-distribution** (see also example below).\n", "\n", "Note, that one can also calculate uncertainty statistics on distributions only as well.\n", "\n", "### Uncertainty quantification\n", "\n", "* Expected value: `cp.E`\n", "\n", "* Variance: `cp.Var`\n", "\n", "* Standard deviation: `cp.Std`\n", "\n", "* Curtosis: `cp.Kurt`\n", "\n", "* Skewness: `cp.Skew`\n", "\n", "* Distribution of Y: `cp.QoI_Dist`\n", "\n", "* Prediction intervals: `cp.Perc`, which is a method to calculate percentiles: an additional argument defining the percentiles needs to be passed.\n", "\n", "If multiple quantities of interest are available:\n", "\n", "* Covariance matrix: `cp.Cov`\n", "\n", "* Correlation matrix: `cp.Corr`\n", "\n", "* Spearman correlation: `cp.Spearman`\n", "\n", "* Auto-correlation function: `cp.Acf`" ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# example uq\n", "exp_reg = cp.E(gpce_regression, joint_distribution)\n", "exp_ps = cp.E(gpce_quadrature, joint_distribution)\n", "\n", "std_reg = cp.Std(gpce_regression, joint_distribution)\n", "str_ps = cp.Std(gpce_quadrature, joint_distribution)\n", "\n", "prediction_interval_reg = cp.Perc(gpce_regression, [5, 95], joint_distribution)\n", "prediction_interval_ps = cp.Perc(gpce_quadrature, [5, 95], joint_distribution)\n", "\n", "print(\"Expected values Standard deviation 90 % Prediction intervals\\n\")\n", "print(' E_reg | E_ps std_reg | std_ps pred_reg | pred_ps')\n", "print(' {} | {} {:>6.3f} | {:>6.3f} {} | {}'.format(exp_reg,\n", " exp_ps,\n", " std_reg,\n", " str_ps,\n", " [\"{:.3f}\".format(p) for p in prediction_interval_reg],\n", " [\"{:.3f}\".format(p) for p in prediction_interval_ps]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Sensitivity analysis\n", "\n", "The variance bases sensitivity indices can be calculated directly from the expansion.\n", "The `chaospy` package provides the following methods:\n", "\n", "* first order indices: `cp.Sens_m`\n", "\n", "* second order indices: `cp.Sens_m2`\n", "\n", "* total indices: `cp.Sens_t`\n", "\n", "Here is an example for the first and total indices for both expansions:" ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# example sens\n", "sensFirst_reg = cp.Sens_m(gpce_regression, joint_distribution)\n", "sensFirst_ps = cp.Sens_m(gpce_quadrature, joint_distribution)\n", "\n", "sensT_reg = cp.Sens_t(gpce_regression, joint_distribution)\n", "sensT_ps = cp.Sens_t(gpce_quadrature, joint_distribution)\n", "\n", "print(\"First Order Indices Total Sensitivity Indices\\n\")\n", "print(' S_reg | S_ps ST_reg | ST_ps \\n')\n", "for k, (s_reg, s_ps, st_reg, st_ps) in enumerate(zip(sensFirst_reg, sensFirst_ps, sensT_reg, sensT_ps)):\n", " print('S_{} : {:>6.3f} | {:>6.3f} ST_{} : {:>6.3f} | {:>6.3f}'.format(k, s_reg, s_ps, k, st_reg, st_ps))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# References\n", "\n", " 1.
**D. Xiu**. \n", " *Numerical Methods for Stochastic Computations: a Spectral Method Approach*,\n", " Princeton University Press,\n", " 2010.\n", "\n", " 2.
**R. C. Smith**. \n", " *Uncertainty Quantification: Theory, Implementation, and Applications*,\n", " *Computational science and engineering series*,\n", " Society for Industrial and Applied Mathematics,\n", " 2013.\n", "\n", " 3.
**J. Feinberg and H. P. Langtangen**. \n", " Chaospy: An open source tool for designing methods of uncertainty quantification,\n", " *Journal of Computational Science*,\n", " 11,\n", " pp. 46-57,\n", " 2015." ] } ], "metadata": {}, "nbformat": 4, "nbformat_minor": 2 }