Emulation

p

We may regard the present state of the universe as the effect of its past and the cause of its future. An intellect which at a certain moment would know all forces that set nature in motion, and all positions of all items of which nature is composed, …

… if this intellect were also vast enough to submit these data to analysis, it would embrace in a single formula the movements of the greatest bodies of the universe and those of the tiniest atom; for such an intellect nothing would be uncertain and the future just like the past would be present before its eyes.

— Pierre Simon Laplace (Laplace, 1814)

Game of Life

Loneliness

loneliness

Crowding

overcrowding

Birth

birth

Glider

Glider (1969)

Loafer

Loafer (2013)

data + model \overset{compute}{\to} prediction

$\text{data} + \text{model} \stackrel{\text{compute}}{\rightarrow} \text{prediction}$

Laplace’s Gremlin

Philosophical Essay on Probabilities Laplace (1814) pg 5

Simulation System

Monolithic System

Service Oriented Architecture

Experiment, Analyze, Design

A Vision

We don’t know what science we’ll want to do in five years’ time, but we won’t want slower experiments, we won’t want more expensive experiments and we won’t want a narrower selection of experiments.

What do we want?

Faster, cheaper and more diverse experiments.

Objective

A two order of magnitude increase in number of experiments run on supply chain three years’ time.

Statistical Emulation

Emulation

GPy and Emulation

f (x) = a (x_{2} - b x_{1}^{2} + c x_{1} - r)^{2} + s (1 - t \cos (x_{1})) + s

$f(\mathbf{ x}) = a(x_2 - bx_1^2 + cx_1 - r)^2 + s(1-t \cos(x_1)) + s$

Computation of $\left\langle f(U)\right\rangle$

Compute the expectation of the Branin function.
Sample at 25 points on a grid.
Compute the mean (Reiemann sum approximation).

Emulate and do MC Samples

Alternatively, we can fit a GP model and compute the integral of the best predictor by Monte Carlo sampling.

Branin Function Fit

Uncertainty Quantification

History of interest, see e.g. McKay et al. (1979)
The review:
- Random Sampling
- Stratified Sampling
- Latin Hypercube Sampling
As approaches for Monte Carlo estimates

Random Sampling

Let the input values $\mathbf{ x}_1, \dots, \mathbf{ x}_n$ be a random sample from $f(\mathbf{ x})$ . This method of sampling is perhaps the most obvious, and an entire body of statistical literature may be used in making inferences regarding the distribution of $Y(t)$ .

Stratified Sampling

Using stratified sampling, all areas of the sample space of $\mathbf{ x}$ are represented by input values. Let the sample space $S$ of $\mathbf{ x}$ be partitioned into $I$ disjoint strata $S_t$ . Let $\pi = P(\mathbf{ x}\in S_i)$ represent the size of $S_i$ . Obtain a random sample $\mathbf{ x}_{ij}$ , $j = 1, \dots, n$ from $S_i$ . Then of course the $n_i$ sum to $n$ . If $I = 1$ , we have random sampling over the entire sample space.

Latin Hypercube Sampling

The same reasoning that led to stratified sampling, ensuring that all portions of $S$ were sampled, could lead further. If we wish to ensure also that each of the input variables $\mathbf{ x}_k$ has all portions of its distribution represented by input values, we can divide the range of each $\mathbf{ x}_k$ into $n$ strata of equal marginal probability $1/n$ , and sample once from each stratum. Let this sample be $\mathbf{ x}_{kj}$ , $j = 1, \dots, n$ . These form the $\mathbf{ x}_k$ component, $k = 1, \dots , K$ , in $\mathbf{ x}_i$ , $i = 1, \dots, n$ . The components of the various $\mathbf{ x}_k$ ’s are matched at random. This method of selecting input values is an extension of quota sampling (Steinberg 1963), and can be viewed as a $K$ -dimensional extension of Latin square sampling (Raj 1968).

Emukit

https://emukit.github.io/emukit/

Emukit

Led by: Javier Gonzalez and Andrei Paleyes
- Team: Mark Pullin, Maren Mahsereci, Alex Gessner, Aaron Klein, Henry Moss and David-Elias Künstle.
- Management: Cliff McCollum & Neil
Available on Github
- Example sensitivity notebook, documentation https://emukit.readthedocs.io/en/latest/

Modular Design

Introduce your own surrogate models.

from emukit.model_wrappers import GPyModelWrapper

To building your own model see this notebook.

from emukit.model_wrappers import YourModelWrapperHere

Emukit Playground

Work Adam Hirst, Software Engineering Intern and Cliff McCollum.
Tutorial on emulation.

Emukit Playground

Emukit Sensitivity Analysis

A possible definition of sensitivity analysis is the following: The study of how uncertainty in the output of a model (numerical or otherwise) can be apportioned to different sources of uncertainty in the model input (Saltelli et al., 2004). A related practice is ‘uncertainty analysis,’ which focuses rather on quantifying uncertainty in model output. Ideally, uncertainty and sensitivity analyses should be run in tandem, with uncertainty analysis preceding in current practice.

In Chapter 1 of Saltelli et al. (2008)

Local Sensitivity

\frac{\partial}{\partial x_{i}} g (x) .

$\frac{\partial}{\partial x_i} g(\mathbf{ x}).$

Global Sensitivity Analysis

var (g (x)) = {⟨ g (x)^{2} ⟩}_{p (x)} - {⟨ g (x) ⟩}_{p (x)}^{2},

$\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,$

{⟨ h (x) ⟩}_{p (x)} = \int_{x} h (x) p (x) d x

$\left\langle h(\mathbf{ x}) \right\rangle _{p(\mathbf{ x})} = \int_\mathbf{ x}h(\mathbf{ x}) p(\mathbf{ x}) \text{d}\mathbf{ x}$

Input Density

p (x) = \prod_{i = 1}^{p} p (x_{i})

$p(\mathbf{ x}) = \prod_{i=1}^pp(x_i)$

x_{i} \sim U (0, 1) .

$x_i \sim \mathcal{U}\left(0,1\right).$

Hoeffding-Sobol Decomposition

\begin{aligned} g (x) = & g_{0} + \sum_{i = 1}^{p} g_{i} (x_{i}) + \sum_{i < j}^{p} g_{i j} (x_{i}, x_{j}) + \dots \\ + g_{1, 2, \dots, p} (x_{1}, x_{2}, \dots, x_{p}), \end{aligned}

$\begin{align*} g(\mathbf{ x}) = & g_0 + \sum_{i=1}^pg_i(x_i) + \sum_{i<j}^{p} g_{ij}(x_i,x_j) + \cdots \\ & + g_{1,2,\dots,p}(x_1,x_2,\dots,x_p), \end{align*}$

g_{0} = {⟨ g (x) ⟩}_{p (x)}

$g_0 = \left\langle g(\mathbf{ x}) \right\rangle _{p(\mathbf{ x})}$

g_{i} (x_{i}) = {⟨ g (x) ⟩}_{p (x_{\sim i})} - g_{0},

$g_i(x_i) = \left\langle g(\mathbf{ x}) \right\rangle _{p(\mathbf{ x}_{\sim i})} - g_0,$

p (x_{\sim i}) = \int p (x) d x_{i}

$p(\mathbf{ x}_{\sim i}) = \int p(\mathbf{ x}) \text{d}x_i$

g_{i, j} (x_{i}, x_{j}) = {⟨ g (x) ⟩}_{p (x_{\sim i, j})} - g_{i} (x_{i}) - g_{j} (x_{j}) - g_{0}

$g_{i,j}(x_i, x_j) = \left\langle g(\mathbf{ x}) \right\rangle _{p(\mathbf{ x}_{\sim i,j})} - g_i(x_i) - g_j(x_j) - g_0$

\begin{aligned} var (g) = & {⟨ g (x)^{2} ⟩}_{p (x)} - {⟨ g (x) ⟩}_{p (x)}^{2} \\ = & {⟨ g (x)^{2} ⟩}_{p (x)} - g_{0}^{2} \\ = & \sum_{i = 1}^{p} var (g_{i} (x_{i})) + \sum_{i < j}^{p} var (g_{i j} (x_{i}, x_{j})) + \dots \\ + var (g_{1, 2, \dots, p} (x_{1}, x_{2}, \dots, x_{p})) . \end{aligned}

$\begin{align*} \text{var}(g) = & \left\langle g(\mathbf{ x})^2 \right\rangle _{p(\mathbf{ x})} - \left\langle g(\mathbf{ x}) \right\rangle _{p(\mathbf{ x})}^2 \\ = & \left\langle g(\mathbf{ x})^2 \right\rangle _{p(\mathbf{ x})} - g_0^2\\ = & \sum_{i=1}^p\text{var}\left(g_i(x_i)\right) + \sum_{i<j}^{p} \text{var}\left(g_{ij}(x_i,x_j)\right) + \cdots \\ & + \text{var}\left(g_{1,2,\dots,p}(x_1,x_2,\dots,x_p)\right). \end{align*}$

Sobol Indices

S_{ℓ} = \frac{var (g (x_{ℓ}))}{var (g (x))},

$S_\ell = \frac{\text{var}\left(g(\mathbf{ x}_\ell)\right)}{\text{var}\left(g(\mathbf{ x})\right)},$

Example: the Ishigami function

Ishigami Function

g (x) = \sin (x_{1}) + a \sin^{2} (x_{2}) + b x_{3}^{4} \sin (x_{1}) .

$g(\textbf{x}) = \sin(x_1) + a \sin^2(x_2) + b x_3^4 \sin(x_1).$

Total Variance

First Order Sobol Indices using Monte Carlo

S_{i} = \frac{var (g_{i} (x_{i}))}{var (g (x))} .

$S_i = \frac{\text{var}\left(g_i(x_i)\right)}{\text{var}\left(g(\mathbf{ x})\right)}.$

Total Effects Using Monte Carlo

Computing the Sensitivity Indices Using the Output of a Model

Conclusions

Sobol indices tool for explaining variance of output as coponents of input variables.

Catapult Simulation

x_{i} = [\begin{matrix} rotation_axis \\ arm_stop \\ spring_binding_1 \\ spring_binding_2 \end{matrix}]

$\mathbf{ x}_i = \begin{bmatrix} \texttt{rotation_axis} \\ \texttt{arm_stop} \\ \texttt{spring_binding_1} \\ \texttt{spring_binding_2} \end{bmatrix}$

Experimental Design for the Catapult

First build an experimental design loop.
Start with model free design.
- Random design
- Latin hypercube sambling
- Sobol sequences
- Orthogonal design

Model Based Design

Use integrated variance reduction.
Try afterwards with uncertainty sampling.

Sensitivity Analysis of a Catapult Simulation

Thanks!

References

Laplace, P.S., 1814. Essai philosophique sur les probabilités, 2nd ed. Courcier, Paris.

McKay, M.D., Beckman, R.J., Conover, W.J., 1979. A comparison of three methods for selecting values of input variables in the analysis of output from a computer code. Technometrics 21, 239–245.

Saltelli, A., Ratto, M., Andres, T., Campolongo, F., Cariboni, J., Gatelli, D., Saisana, M., Tarantola, S., 2008. Global sensitivity analysis: The primer. wiley.

Saltelli, A., Tarantola, S., Campolongo, F., Ratto, M., 2004. Sensitivity analysis in practice: A guide to assessing scientific methods. wiley.