Vinzenz Gregor Eck, Expert Analytics, Oslo
**Jacob Sturdy**, Department of Structural Engineering, NTNUDate: Jul 13, 2018
# ipython magic
%matplotlib notebook
%load_ext autoreload
%autoreload 2
# plot configuration
import matplotlib
import matplotlib.pyplot as plt
plt.style.use("ggplot")
# import seaborn as sns # sets another style
matplotlib.rcParams['lines.linewidth'] = 3
fig_width, fig_height = (7.0,5.0)
matplotlib.rcParams['figure.figsize'] = (fig_width, fig_height)
# font = {'family' : 'sans-serif',
# 'weight' : 'normal',
# 'size' : 18.0}
# matplotlib.rc('font', **font) # pass in the font dict as kwar
import chaospy as cp
import numpy as np
from numpy import linalg as LA
To conduct UQSA with the polynomial chaos (PC) method we use the python package chaospy
.
This package includes many features and methods to apply non-intrusive polynomial chaos to any model with
few lines of code.
Since, chaospy
allows for convenient definition of random variables, calculation of joint distributions and generation of samples,
we apply chaospy
as well for the Monte Carlo method.
In the following we will briefly describe the theory of polynomial chaos and give a practical step by step guide
for the application of chaospy
.
For a more in depth introduction to the polynomial chaos method see ([[xiu_numerical_2010;@ smith_uncertainty_2013]](#xiu_numerical_2010;@ smith_uncertainty_2013)).
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. For simplicity, we consider a model function $f$ which takes random variables $\mathbf{Z}$ and non-random variables as spatial coordinates $x$ or time $t$, as input:
Then we seek to generate a polynomial expansion such that we can approximate the model response.
where $v_i$ are the expansion coefficients, which are only dependent on the non-stochastic parameters, and $\Phi_i$ are orthogonal polynomials, which depend solely on the random input parameter $\mathbf{Z}$. Once the polynomial expansion is calculated, we can analytically calculate the statistics of approximated model output.
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.
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...).
As state above the polynomial chaos expansion consists of a sum of basis polynomials in the input parameters $\mathbf{Z}$. 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. 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.
The following equalities hold for the orthogonal basis polynomials used in polynomial chaos
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. 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$.
Once the random inputs $\mathbf{Z}$ are properly defined with marginal distributions, orthogonal polynomials can be constructed. For most univariate distributions, polynomial bases functions are defined, and listed in the Asker Wilkeys scheme. 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.
There are two non-intrusive ways of approximating a polynomial chaos expansion coefficients:
Supposing a polynomial expansion approximation $y_{N} = \sum_i v_i \Phi_i$, 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$. The task is thus to find suitable coefficients $v_i$ such that this condition is satisfied. 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 formed into linear system of equations for the coefficients $v_i$ at these nodes:
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}$. This is done by choosing a larger number of samples so that (3) is over determined and then minimizing the appropriate error ${\lvert\lvert {y_{N} \rvert\rvert}-y}_{R}$ over $\left\{\mathbf{z}^{(s)}\right\}_{s=1}^{P}$. Ordinary least squares, Ridge-Regression, and Tikhonov-Regularization are all regression methods that may be applied to this problem.
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
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.
This results in an approximation $\tilde{v}_{i}$ of $v_{i}$ the error between the final approximation may be split
where the first term is called the truncation error and the second term the quadrature error. 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.
Once the expansion was generated, it can be used directly to calculated statistics for the uncertainty and sensitivity analysis. The two most common measures for uncertainty quantification, expected value and variance, can be calculated by inserting the expansion into the definition of the measures.
The expected value is equal to the first expansion coefficient:
The variance is the sum of squared expansion coefficients multiplied by normalisation constants of the polynomials:
(Note the orthogonality of individual terms implies their covariance is zero, thus the variance is simply the sum of the variances of the terms.)
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
and similarly one may define $\mathcal{A}_{ij}$ for pairwise combinations of inputs and further to calculate all orders of sensitivities.
The python package chaospy
an introductory paper to the package including a comparison to other software packages is
presented here ([feinberg_2015]).
You can find an introduction, tutorials and the source code at the projects homepage: https://github.com/jonathf/chaospy.
The installation of the package can be done via pip
:
pip install chaospy
In the following we will use the import naming convention of the package creator to import the package in python:
import chaospy as cp
Therefore it will be convenient to see whenever a method of the package is applied.
The package chaospy
is doc-string annotated which means that every method provides a short help text with small examples.
To show the method documentation simply type a ?
after the method name in a ipython console or notebook.
As shown in the following two examples:
# show help for uniform distributions
cp.Uniform?
# show help for sample generation
cp.samplegen?
To conduct UQSA analysis with polynomial chaos we need to follow the following steps:
Definition of the marginal and joint distributions
Generation of the orthogonal polynomials
Linear regression
Generation of samples
Evaluation of the model for all samples
Generation of the polynomial chaos expansion
Pseudo-spectral projection
Generation of integration nodes and weights
Evaluation of the model for all nodes
Generation of the polynomial chaos expansion
Calculations of all statistics
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 generate a set of points in the parameter space where the model must be evaluated (steps 3.b and 4.b, respectively).
The analysis of a each model starts with the definition of the marginal distributions for each random model input, i.e. describing it as
random variable.
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
.
The help function can be used to find out more about the required arguments, e.g. help(cp.Normal)
.
In the following an example for 3 random variables with uniform, normal and log-normal distribution:
# simple distributions
rv1 = cp.Uniform(0, 1)
rv2 = cp.Normal(0, 1)
rv3 = cp.Lognormal(0, 1, 0.2, 0.8)
print(rv1, rv2, rv3)
After all random input variables are defined with univariate random variables a multi-variate random variable and its joint distribution can be constructed with the following command:
# joint distributions
joint_distribution = cp.J(rv1, rv2, rv3)
print(joint_distribution)
It is also possible to construct independent identical distributed random variables from any univariate variable:
# creating iid variables
X = cp.Normal()
Y = cp.Iid(X, 4)
print(Y)
All 64 distributions available in the chaospy package can be found in the following table:
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 |
The orthogonal polynomials can be generated with different methods, in chaospy
there are 4 methods implemented. The
most stable method, and therefore most advised is the three terms recursion method.
Orthogonalization Method | |
---|---|
Cholesky decomposition | cp.orth\_chol |
Three terms recursion | cp.orth\_ttr |
Modified Gram-Schmidt | cp.orth\_gs |
Regarding the three terms recursion method: For the distributions Normal, Uniform, Gamma, Log-normal, Triangle, Beta and stochastic independent variable combinations of those, the three terms recursion coefficients are known. For all other distributions the coefficients are estimated numerically. The three terms recursion method is then also called discretized stieltjes method.
The most stable method and therefore most applied method is the three terms recursion (discretized stieltjes method) method.
We will look at all in a small example, try to increase the polynomial order and the instabilities of the methods become visible.
# example orthogonalization schemes
# a normal random variable
n = cp.Normal(0, 1)
x = np.linspace(0,1, 50)
# the polynomial order of the orthogonal polynomials
polynomial_order = 3
poly = cp.orth_chol(polynomial_order, n, normed=True)
print('Cholesky decomposition {}'.format(poly))
ax = plt.subplot(131)
ax.set_title('Cholesky decomposition')
_=plt.plot(x, poly(x).T)
_=plt.xticks([])
poly = cp.orth_ttr(polynomial_order, n, normed=True)
print('Discretized Stieltjes / Three terms reccursion {}'.format(poly))
ax = plt.subplot(132)
ax.set_title('Discretized Stieltjes ')
_=plt.plot(x, poly(x).T)
poly = cp.orth_gs(polynomial_order, n, normed=True)
print('Modified Gram-Schmidt {}'.format(poly))
ax = plt.subplot(133)
ax.set_title('Modified Gram-Schmidt')
_=plt.plot(x, poly(x).T)
The linear regression method requires to conduct the three following steps:
Generation of samples
Evaluation of the model for all samples
Generation of the polynomial chaos expansion
In the following we will not consider the model evaluation.
Once a random variable is defined or a joint random variable, also referred as distribution here, the following method can be used to generate as set of samples:
# sampling in chaospy
u = cp.Uniform(0,1)
u.sample?
The method takes the arguments size which is the number of samples and rule which is the applied sampling scheme. The following example shows the creation of 2 set of samples for the sampling schemes (Pseudo-)Random and Hammersley.
# example sampling
u1 = cp.Uniform(0,1)
u2 = cp.Uniform(0,1)
joint_distribution = cp.J(u1, u2)
number_of_samples = 350
samples_random = joint_distribution.sample(size=number_of_samples, rule='R')
samples_hammersley = joint_distribution.sample(size=number_of_samples, rule='M')
fig1, ax1 = plt.subplots()
ax1.set_title('Random')
ax1.scatter(*samples_random)
ax1.set_xlabel("Uniform 1")
ax1.set_ylabel("Uniform 2")
ax1.axis('equal')
fig2, ax2 = plt.subplots()
ax2.set_title('Hammersley sampling')
ax2.scatter(*samples_hammersley)
ax2.set_xlabel("Uniform 1")
ax2.set_ylabel("Uniform 2")
ax2.axis('equal')
All sampling schemes implemented in chaospy are listed in the following table:
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 |
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:
# example save samples to file
# Creates a csv file where each row corresponds to the sample number and each column with teh variables in the joint distribution
csv_file = "csv_samples.csv"
sep = '\t'
header = ["u1", "u2"]
header = sep.join(header)
np.savetxt(csv_file, samples_random, delimiter=sep, header=header)
Each row of the csv file now contains a single sample from the joint distribution with the columns corresponding to each component.
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.
# example load samples from file
# loads a csv file where the samples/or model evaluations for each sample are saved
# with one sample per row. Multiple components ofoutput can be stored as separate columns
filepath = "external_evaluations.csv"
data = np.loadtxt(filepath)
After the model is evaluated for all samples, the polynomial chaos expansion can be generated with the following method:
# linear regression in chaospy
cp.fit_regression?
In the following we show a complete example for polynomial chaos expansion using the linear regression. The model applied the very simple mathematical expression:
The random variables for $Z_1, Z_2$ are defined as simple uniform random variables:
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}$.
Here is the annotated example code with all steps required to generate a polynomial chaos expansion with linear regression:
# example linear regression
# 1. define marginal and joint distributions
u1 = cp.Uniform(0,1)
u2 = cp.Uniform(0,1)
joint_distribution = cp.J(u1, u2)
# 2. generate orthogonal polynomials
polynomial_order = 3
poly = cp.orth_ttr(polynomial_order, joint_distribution)
# 3.1 generate samples
number_of_samples = cp.bertran.terms(polynomial_order, len(joint_distribution))
samples = joint_distribution.sample(size=number_of_samples, rule='R')
# 3.2 evaluate the simple model for all samples
model_evaluations = samples[0]+samples[1]*samples[0]
# 3.3 use regression to generate the polynomial chaos expansion
gpce_regression = cp.fit_regression(poly, samples, model_evaluations)
print("Success")
# quadrature in polychaos
cp.generate_quadrature?
We will look at the following arguments of the method: order is the order of the quadrature, domain is the , rule is the name or key of the quadrature rule to apply.
In the following example we look at some quadrature nodes for the same uniform variables as for the sampling, for Optimal Gaussian quadrature and Clenshaw-Curtis quadrature.
# example quadrature
u1 = cp.Uniform(0,1)
u2 = cp.Uniform(0,1)
joint_distribution = cp.J(u1, u2)
order = 5
nodes_gaussian, weights_gaussian = cp.generate_quadrature(order=order, domain=joint_distribution, rule='G')
nodes_clenshaw, weights_clenshaw = cp.generate_quadrature(order=order, domain=joint_distribution, rule='C')
print('Number of nodes gaussian quadrature: {}'.format(len(nodes_gaussian[0])))
print('Number of nodes clenshaw-curtis quadrature: {}'.format(len(nodes_clenshaw[1])))
fig1, ax1 = plt.subplots()
ax1.scatter(*nodes_gaussian, marker='o', color='b')
ax1.scatter(*nodes_clenshaw, marker= 'x', color='r')
ax1.set_xlabel("Uniform 1")
ax1.set_ylabel("Uniform 2")
ax1.axis('equal')
In the following all quadrature rules implemented in chaospy are highlighted:
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 |
In the following example we show sparse vs. normal quadrature nodes:
# example sparse grid quadrature
u1 = cp.Uniform(0,1)
u2 = cp.Uniform(0,1)
joint_distribution = cp.J(u1, u2)
order = 2
# sparse grid has exponential growth, thus a smaller order results in more points
nodes_clenshaw, weights_clenshaw = cp.generate_quadrature(order=order, domain=joint_distribution, rule='C')
nodes_clenshaw_sparse, weights_clenshaw_sparse = cp.generate_quadrature(order=order, domain=joint_distribution, rule='C', sparse=True)
print('Number of nodes normal clenshaw-curtis quadrature: {}'.format(len(nodes_clenshaw[0])))
print('Number of nodes clenshaw-curtis quadrature with sparse grid : {}'.format(len(nodes_clenshaw_sparse[0])))
fig1, ax1 = plt.subplots()
ax1.scatter(*nodes_clenshaw, marker= 'x', color='r')
ax1.scatter(*nodes_clenshaw_sparse, marker= 'o', color='b')
ax1.set_xlabel("Uniform 1")
ax1.set_ylabel("Uniform 2")
ax1.axis('equal')
After the model is evaluated for all integration nodes, the polynomial chaos expansion can be generated with the following method:
# spectral projection in chaospy
cp.fit_quadrature?
In the following we show again a complete example for polynomial chaos expansion using the pseudo spectral approach to calculate the expansion coefficients. The model applied the same simple mathematical expression as before:
The random variables for $Z_1, Z_2$ are defined as simple uniform random variables:
# example spectral projection
# 1. define marginal and joint distributions
u1 = cp.Uniform(0,1)
u2 = cp.Uniform(0,1)
joint_distribution = cp.J(u1, u2)
# 2. generate orthogonal polynomials
polynomial_order = 3
poly = cp.orth_ttr(polynomial_order, joint_distribution)
# 4.1 generate quadrature nodes and weights
order = 5
nodes, weights = cp.generate_quadrature(order=order, domain=joint_distribution, rule='G')
# 4.2 evaluate the simple model for all nodes
model_evaluations = nodes[0]+nodes[1]*nodes[0]
# 4.3 use quadrature to generate the polynomial chaos expansion
gpce_quadrature = cp.fit_quadrature(poly, nodes, weights, model_evaluations)
print("Success")
Once the polynomial chaos expansion is created either with pseudo-spectral projection or with regression method The calculation of statistics is straight forward. The following listing gives an overview of all available methods take all the same input parameter the polynomial-expansion and the joint-distribution (see also example below).
Note, that one can also calculate uncertainty statistics on distributions only as well.
Expected value: cp.E
Variance: cp.Var
Standard deviation: cp.Std
Curtosis: cp.Kurt
Skewness: cp.Skew
Distribution of Y: cp.QoI_Dist
Prediction intervals: cp.Perc
, which is a method to calculate percentiles: an additional argument defining the percentiles needs to be passed.
If multiple quantities of interest are available:
Covariance matrix: cp.Cov
Correlation matrix: cp.Corr
Spearman correlation: cp.Spearman
Auto-correlation function: cp.Acf
# example uq
exp_reg = cp.E(gpce_regression, joint_distribution)
exp_ps = cp.E(gpce_quadrature, joint_distribution)
std_reg = cp.Std(gpce_regression, joint_distribution)
str_ps = cp.Std(gpce_quadrature, joint_distribution)
prediction_interval_reg = cp.Perc(gpce_regression, [5, 95], joint_distribution)
prediction_interval_ps = cp.Perc(gpce_quadrature, [5, 95], joint_distribution)
print("Expected values Standard deviation 90 % Prediction intervals\n")
print(' E_reg | E_ps std_reg | std_ps pred_reg | pred_ps')
print(' {} | {} {:>6.3f} | {:>6.3f} {} | {}'.format(exp_reg,
exp_ps,
std_reg,
str_ps,
["{:.3f}".format(p) for p in prediction_interval_reg],
["{:.3f}".format(p) for p in prediction_interval_ps]))
The variance bases sensitivity indices can be calculated directly from the expansion.
The chaospy
package provides the following methods:
first order indices: cp.Sens_m
second order indices: cp.Sens_m2
total indices: cp.Sens_t
Here is an example for the first and total indices for both expansions:
# example sens
sensFirst_reg = cp.Sens_m(gpce_regression, joint_distribution)
sensFirst_ps = cp.Sens_m(gpce_quadrature, joint_distribution)
sensT_reg = cp.Sens_t(gpce_regression, joint_distribution)
sensT_ps = cp.Sens_t(gpce_quadrature, joint_distribution)
print("First Order Indices Total Sensitivity Indices\n")
print(' S_reg | S_ps ST_reg | ST_ps \n')
for k, (s_reg, s_ps, st_reg, st_ps) in enumerate(zip(sensFirst_reg, sensFirst_ps, sensT_reg, sensT_ps)):
print('S_{} : {:>6.3f} | {:>6.3f} ST_{} : {:>6.3f} | {:>6.3f}'.format(k, s_reg, s_ps, k, st_reg, st_ps))