Sequential Computer Experimental Design for Estimating an Extreme Probability or Quantile

TL;DR

This paper develops a sequential experimental design method to efficiently estimate extreme probabilities or quantiles of a system's output using a limited number of computer simulations, improving accuracy and practicality.

Contribution

It introduces a sequential design approach for estimating tail probabilities with fewer simulations, comparing criteria and addressing initial design and input modeling.

Findings

Effective sequential design criteria identified

Method reduces number of simulations needed

Guidelines for initial design and input modeling provided

Abstract

A computer code can simulate a system's propagation of variation from random inputs to output measures of quality. Our aim here is to estimate a critical output tail probability or quantile without a large Monte Carlo experiment. Instead, we build a statistical surrogate for the input-output relationship with a modest number of evaluations and then sequentially add further runs, guided by a criterion to improve the estimate. We compare two criteria in the literature. Moreover, we investigate two practical questions: how to design the initial code runs and how to model the input distribution. Hence, we close the gap between the theory of sequential design and its application.

Tables

Table 1. Table 1 : RMSE of the final ( n = 40 𝑛 40 n=40 ) probability estimate for the short column function.

Initial design	EI	Discrepancy
Random	$0.00439$	$0.00014$
Uniform	$0.00094$	$0.00011$

Table 2. Table 2 : RMSE of the final ( n = 40 𝑛 40 n=40 ) quantile estimate for the short column function.

Initial design	EI	Discrepancy
Random	$0.07592$	$0.01287$
Uniform	$0.04999$	$0.00991$

Table 3. Table 3 : Functional ANOVA for the floor-system application showing the percentage contributions from the eight main effects to the total variation of the GP prediction over the 8 8 8 -dimensional input space. The main effects explain 99.2 % percent 99.2 99.2\% of the total variance in the GP prediction. Interactions between two input factors are not shown as none of them contributes more than 0.5%.

Input (MOE)	% variation
$x_1$	0.0
$x_2$	5.3
$x_3$	12.8
$x_4$	26.2
$x_5$	37.5
$x_6$	16.1
$x_7$	1.0
$x_8$	0.2

Equations

$$\hat{F}(w)=\frac{1}{N}\sum _{i=1}^N\text{mathds}1(y({\mathbf{x}}_{\text{MC}}^{(i)})<w),$$

$$\mu +Z(\mathbf{x}).$$

$$R(\mathbf{x},{\mathbf{x}}^{\prime })=\exp \left(-\sum _{j=1}^d{\theta }_j|x_j-x_j^{\prime }{|}^{p_j}\right),$$

$$N\left(m({\mathbf{x}}^{\ast }),v({\mathbf{x}}^{\ast })\right),$$

$$m({\mathbf{x}}^{\ast })=\mu +{\mathbf{r}}^T({\mathbf{x}}^{\ast }){\mathbf{R}}^{-1}(\mathbf{y}-1\mu ),$$

$$v({\mathbf{x}}^{\ast })={\sigma }^2\left(1-{\mathbf{r}}^T({\mathbf{x}}^{\ast }){\mathbf{R}}^{-1}\mathbf{r}({\mathbf{x}}^{\ast })\right).$$

$$t_{n-1}\left({\hat{m}}_{\mathbf{\psi }}({\mathbf{x}}^{\ast }),{\hat{v}}_{\mathbf{\psi }}({\mathbf{x}}^{\ast })\right),$$

$${\hat{m}}_{\mathbf{\psi }}({\mathbf{x}}^{\ast })={\hat{\mu }}_{\mathbf{\psi }}+{\mathbf{r}}^T({\mathbf{x}}^{\ast }){\mathbf{R}}^{-1}(\mathbf{y}-1{\hat{\mu }}_{\mathbf{\psi }}),$$

$${\hat{\mu }}_{\mathbf{\psi }}=\frac{1^T{\mathbf{R}}^{-1}\mathbf{r}({\mathbf{x}}^{\ast })}{{\mathbf{r}}^T({\mathbf{x}}^{\ast }){\mathbf{R}}^{-1}\mathbf{r}({\mathbf{x}}^{\ast })},$$

$${\hat{v}}_{\mathbf{\psi }}({\mathbf{x}}^{\ast })={\widehat{{\sigma }^2}}_{\mathbf{\psi }}\left(1-{\mathbf{r}}^T({\mathbf{x}}^{\ast }){\mathbf{R}}^{-1}\mathbf{r}({\mathbf{x}}^{\ast })+\frac{{\left(1-1^T{\mathbf{R}}^{-1}\mathbf{r}({\mathbf{x}}^{\ast })\right)}^2}{1^T{\mathbf{R}}^{-1}1}\right),$$

$${\widehat{{\sigma }^2}}_{\mathbf{\psi }}=\frac{(\mathbf{y}-1{\hat{\mu }}_{\mathbf{\psi }}{)}^T{\mathbf{R}}^{-1}(\mathbf{y}-1{\hat{\mu }}_{\mathbf{\psi }})}{n-1}.$$

$$\hat{m}({\mathbf{x}}^{\ast })=\frac{1}{M}\sum _{i=1}^M{\hat{m}}_{{\mathbf{\psi }}^{(i)}}({\mathbf{x}}^{\ast }),$$

$$\hat{v}({\mathbf{x}}^{\ast })=\frac{1}{M}\sum _{i=1}^M{\hat{v}}_{{\mathbf{\psi }}^{(i)}}({\mathbf{x}}^{\ast })+\frac{1}{M-1}\sum _{i=1}^M{\left({\hat{m}}_{{\mathbf{\psi }}^{(i)}}({\mathbf{x}}^{\ast })-\hat{m}({\mathbf{x}}^{\ast })\right)}^2.$$

$$I({\bf x})=\begin{cases}\alpha^{2}v_{{\boldsymbol{\psi}},\sigma^{2}}({\bf x})-(y({\bf x})-y_{f})^{2}&\mbox{if $|y({\bf x})-y_{f}|<\alpha\sqrt{v_{{\boldsymbol{\psi}},\sigma^{2}}({\bf x})}$}\\ 0&\mbox{otherwise},\end{cases}$$

$$v_{\mathbf{\psi },{\sigma }^2}({\mathbf{x}}^{\ast })={\sigma }^2\left(1-{\mathbf{r}}^T({\mathbf{x}}^{\ast }){\mathbf{R}}^{-1}\mathbf{r}({\mathbf{x}}^{\ast })+\frac{{\left(1-1^T{\mathbf{R}}^{-1}\mathbf{r}({\mathbf{x}}^{\ast })\right)}^2}{1^T{\mathbf{R}}^{-1}1}\right).$$

$$\begin{array}[]{lcl}E(I({\bf x}))&=&\left(\alpha^{2}v_{{\boldsymbol{\psi}},\sigma^{2}}({\bf x})-(\hat{m}_{{\boldsymbol{\psi}}}({\bf x})-y_{f})^{2}\right)(\Phi(u_{2})-\Phi(u_{1}))\\ &+&v_{{\boldsymbol{\psi}},\sigma^{2}}({\bf x})\left((u_{2}\phi(u_{2})-u_{1}\phi(u_{1}))-(\Phi(u_{2})-\Phi(u_{1}))\right)\\ &+&2(\hat{m}_{{\boldsymbol{\psi}}}({\bf x})-y_{f})\sqrt{v_{{\boldsymbol{\psi}},\sigma^{2}}({\bf x})}(\phi(u_{2})-\phi(u_{1})),\end{array}$$

$${\mathbf{x}}^{\ast }=\underset{\mathbf{x}\in {\mathcal{X}}_{\text{cand}}}{\mathrm{argmax}}E(I(\mathbf{x})),$$

$$D({\bf x})=\begin{cases}\frac{(y({\bf x})-y_{f})^{2}+\epsilon}{v_{{\boldsymbol{\psi}},\sigma^{2}}}&\mbox{if $v_{{\boldsymbol{\psi}},\sigma^{2}}>0$}\\ \infty&\mbox{otherwise},\end{cases}$$

$${\mathbf{x}}^{\ast }=\underset{\mathbf{x}\in {\mathcal{X}}_{\text{cand}}}{\mathrm{argmin}}E(D(\mathbf{x})),$$

$$E(D({\bf x}))=\begin{cases}\left(\frac{\hat{m}_{{\boldsymbol{\psi}}}({\bf x})-y_{f}}{\sqrt{v_{{\boldsymbol{\psi}},\sigma^{2}}}}\right)^{2}+\frac{\epsilon}{v_{{\boldsymbol{\psi}},\sigma^{2}}}+1,&\mbox{if $v_{{\boldsymbol{\psi}},\sigma^{2}}\neq 0$}\\ \infty&\mbox{otherwise},\end{cases}$$

$${\mathbf{x}}^{\ast }=\underset{\mathbf{x}\in {\mathcal{X}}_{\text{cand}}}{\mathrm{argmin}}\left(\frac{|{\hat{m}}_{\mathbf{\psi }}(\mathbf{x})-y_f|}{\sqrt{v_{\mathbf{\psi },{\sigma }^2}(\mathbf{x})}}\right).$$

$${\mathbf{x}}^{\ast }=\underset{\mathbf{x}\in {\mathcal{X}}_{\text{cand}}}{\mathrm{argmax}}\text{Pr}\left(Z<-\frac{|{\hat{m}}_{\mathbf{\psi }}(\mathbf{x})-y_f|}{\sqrt{v_{\mathbf{\psi },{\sigma }^2}(\mathbf{x})}}\right),$$

$$y(\mathbf{x})=1-\frac{4x_m}{b{h}^2{x}_z}-\frac{x_p^2}{b^2{h}^2{x}_z^2},$$

$$y(\mathbf{x})=1-\frac{4}{b{h}^2}{t}_1-\frac{1}{b^2{h}^2}{t}_2.$$

$$H(x)=p_1\times G_1(x)+p_2\times G_2(x),$$

$${\hat{p}}_f=\sum _{h=1}^{256}{w}_h\frac{1}{50}\sum _{i=1}^{50}\text{mathds}1(\hat{y}({\mathbf{x}}_{\text{MC}}^{(hi)})>y_f),$$

Full text

\StopCensoring

Sequential Computer Experimental Design for Estimating an Extreme Probability or Quantile

\censor

HAO CHEN

\censorDepartment of Statistics, University of British Columbia, Vancouver, BC, Canada V6T 1Z4

\censorhao.chen@stat.ubc.ca*

\censor

WILLIAM J. WELCH111Corresponding author

\censorDepartment of Statistics, University of British Columbia, Vancouver, BC, Canada V6T 1Z4

\censorwill@stat.ubc.ca*

A computer code can simulate a system’s propagation of variation from random inputs to output measures of quality. Our aim here is to estimate a critical output tail probability or quantile without a large Monte Carlo experiment. Instead, we build a statistical surrogate for the input-output relationship with a modest number of evaluations and then sequentially add further runs, guided by a criterion to improve the estimate. We compare two criteria in the literature. Moreover, we investigate two practical questions: how to design the initial code runs and how to model the input distribution. Hence, we close the gap between the theory of sequential design and its application.

Key Words: Computer Experiments; Expected Improvement; Gaussian Process; Sequential Method.

1 Introduction

Consider a deterministic computer model of an engineering system with $d$ inputs ${\bf x}=(x_{1},\ldots,\allowbreak x_{d})$ and output $y=y({\bf x})$ measuring performance. Although any one run of the computer code is deterministic, to model manufacturing or end-use variation, the values of the inputs may be drawn from random variables ${\bf X}=(X_{1},\ldots,\allowbreak X_{d})$. Hence, we can define $Y=y({\bf X})$ as a random output. This propagation of variation from ${\bf X}$ to $Y$ is our focus. Specifically, without loss of generality, suppose $y$ is a smaller-the better output, and $y_{f}$ is a given critical value, where failure occurs say, and we want to estimate the (small) failure probability $p_{f}=\textrm{Pr}(y({\bf X})>y_{f})$. Quantile estimation is a similar inverse problem: for given small probability $p_{f}$, find $y_{f}$ such that $\textrm{Pr}(y({\bf X})>y_{f})=p_{f}$. The focus of the paper, then, is to provide a good estimate of an extreme tail probability $p_{f}$ or an extreme quantile $y_{f}$. We have in mind a computer model that is too expensive to run a simple Monte Carlo experiment.

The motivating computer model in the paper is a model of a floor system. The inputs are $d$ values of the modulus of elasticity (MOE) of the $d$ floor joists and the output is the maximum deflection under a static load. If the maximum deflection exceeds a critical cut-off, the system will fail. The application will be described further in Section 2.

All methods considered are based on the Monte Carlo (MC) method but with implicit values of the output. First, consider the straightforward MC problem, where it is feasible to evaluate the computer model many times and directly estimate $p_{f}$ or $y_{f}$. There is a large MC set ${\bf\mathcal{X}_{\text{MC}}}=\{{\bf x}_{\text{MC}}^{(1)},\ldots,{\bf x}_{\text{MC}}^{(N)}\}$ sampled from the distribution of ${\bf X}$. Plugging the ${\bf\mathcal{X}_{\text{MC}}}$ points into the computer model, one can obtain the $N$ outputs ${\bf y}_{\text{MC}}=\{y({\bf x}_{\text{MC}}^{(1)}),\ldots,y({\bf x}_{\text{MC}}^{(N)})\}$. Then the empirical distribution of the output $Y$ is given by

[TABLE]

where $\mathds{1}(E)=1$ if event $E$ is true and 0 otherwise. The failure probability or quantile are then the (approximate) solutions of $p_{f}=1-\hat{F}(y_{f})$. According to the Glivenko-Cantelli theorem (Cantelli, 1933; Glivenko, 1933), the supremum of the difference between the empirical distribution and the true CDF converges to 0 almost surely.

However, the computer code may be too computationally expensive to obtain a large MC sample. Consider a limited experimental budget, for example $40$ runs of the code. The empirical distribution of $Y$ based directly on $40$ runs would be inaccurate, especially in the tails. The solution is to use a Gaussian process (GP) model (Sacks et al., 1989; Santner et al., 2003), trained with a modest number of runs, as a fast statistical proxy for the expensive computer model. The key ideas are as follows. First, generate a discrete training set ${\bf\mathcal{X}}=\{{\bf x}^{(1)},\ldots,{\bf x}^{(n)}\}$, where $n$ is relatively small, and obtain ${\bf y}=(y({\bf x}^{(1)}),\ldots,y({\bf x}^{(n)}))^{T}$ by running the computer model. The available data enable one to build a GP model. Second, obtain predicted values in a much larger MC set from the trained GP model: $\hat{{\bf y}}_{\text{MC}}=\{\hat{y}({\bf x}_{\text{MC}}^{(1)}),\ldots,\hat{y}({\bf x}_{\text{MC}}^{(N)})\}$. The failure probability $p_{f}$ is estimated by $1-\hat{F}(y_{f})$, where the true outputs ${\bf y}_{\text{MC}}$ in (1) are replaced by $\hat{{\bf y}}_{\text{MC}}$. The estimate of the quantile can be obtained in an analogous way.

There are several ways to choose the training set for fitting the GP. The first, which we call a fixed design strategy, is to use up all of the design budget $n$, to train one statistical surrogate, which is employed to predict the outputs of the points in the MC set. The fixed strategy is simple, but the estimation accuracy may not be good enough, especially for estimating an extreme tail probability or quantile. The second way, a sequential design strategy, is to use part of the design budget at the beginning to train an initial GP and sequentially add points into the design space (one at a time), guided by a search criterion until the budget has been exhausted. Each time a new point is added, both the surrogate model and the probability/quantile estimate are updated. Compared with the fixed design strategy, the second strategy is a dynamic process that allows new information to be added “on the fly”, and thus should intuitively provide more accurate estimation. A better estimation outcome using sequential methodology is well recognized; see Jones et al. (1998) and Ranjan et al. (2008) for examples.

The core of a sequential method is the search criterion. Jones et al. (1998) proposed the expected improvement (EI) criterion based on an improvement function for global optimization, and EI has been popular during the past 20 years. In general, there are different improvement functions for different statistical objectives. Ranjan et al. (2008) introduced an improvement function for contour estimation, i.e., search for ${\bf x}=(x_{1},\ldots,x_{d})$ such that $y({\bf x})=a$, where $a$ is the target level. Roy and Notz (2014) used two different criteria for estimating a quantile where $p_{f}=0.2$: the first one is also based on EI and used nearly the same improvement function as Ranjan et al. (2008) for estimating the region defining a quantile. The second criterion of Roy and Notz (2014), the so-called hypothesis testing-based criterion, optimizes a “discrepancy” between the current prediction and the target quantile at any untried input set. The details will be reviewed in Section 4.

Taking such a sequential design approach, the focus of this article is to provide recommendations for some practical questions arising when implementing the sequential strategy. The questions, along with summaries of some findings, are as follows.

• •

There are two search criteria for quantile and probability estimation: EI and the hypothesis testing-based criterion. Which yields a more accurate result, especially for an extreme probability or quantile? Roy and Notz (2014) compared the performance of sequential designs and fixed designs and showed sequential designs produced more accurate quantile estimates. However, they did not consider probability estimation. None of their examples covered extreme quantile estimation. In this paper, we explicitly contrast the performances of EI and the hypothesis testing-based criterion for quantile and probability estimation. We conclude that the hypothesis testing-based criterion has a faster convergence to its target and hence is preferred.

• •

How to specify the input distribution for ${\bf X}$? The distribution of $Y=y({\bf X})$ is determined by the distribution of ${\bf X}$ by propagating variation through the deterministic computer model. We explore in Section 6.2 several different ways of modelling the input distribution.

• •

How to generate the MC set? A simple random sample from the input distribution is straightforward, but it might be computationally inefficient even with GP prediction. In Section 6.3, we adopt a stratified sampling scheme to generate a smaller MC set such that an extreme probability or quantile is well estimated.

• •

How to select the initial design for training the starting GP? In Section 5.1 we show that over-sampling the tails of the ${\bf X}$ distribution may be preferable.

• •

What is a suitable stopping criterion for a sequential algorithm? It is critically important to know when to stop in practice when the “true” probability or quantile are unknown. Based on the hypothesis testing-based criterion, we generate diagnostic plots in Section 7 to check if the sequence has converged or not.

• •

How to estimate the unknown parameters of a GP for these purposes? We use the Bayesian method of Chen et al. (2017) instead of maximum likelihood estimation (MLE). EI and the hypothesis testing-based criterion both require the predictive standard deviation of the response at an untried point. Bayesian methods are able to fully quantify parameter estimation uncertainty and are therefore preferred for training a GP model and predicting in sequential design.

Thus, the major contributions are to close the gap between the theory of sequential design and its practical use by answering the above questions. The rest of the paper is organized as follows. In Section 2, we introduce the computer model used to illustrate the methodology. Gaussian processes are reviewed in Section 3. Section 4 defines the sequential algorithm and the aforementioned search criteria in more detail. Section 5 presents an empirical study based on failure of a column, and Section 6 revisits the motivating floor system. Section 7 proposes a diagnostic for stopping. Finally, Section 8 makes some concluding remarks.

2 Computer Model of a Floor System

The application that motivates this research is a numerical model of a floor system (McCutcheon, 1984). The performance characteristic of interest is the floor’s maximum deflection under a static load. The $d$ inputs $x_{1},\ldots,\allowbreak x_{d}$ to the model are the modulus of elasticity (MOE) values of $d$ supporting joists in units of psi. For instance, Figure 1 shows a system with $d=8$ joists. The load could be a further input, but it will be kept constant in the analysis of Section 6 because we are more interested in the relationship between joist MOEs and the floor’s maximum deflection. Given the inputs, the computer model outputs the $d$ deflections of the $d$ joists, and the maximum deflection is taken as the performance measure of interest, $y$ (units in); see Figure 2.

The computer code will be treated as a black-box.

The code is deterministic, but if we think of the inputs $x_{1},\ldots,\allowbreak x_{d}$ as realizations of random variables $X_{1},\ldots,\allowbreak X_{d}$, then the output deflection $y$ is a realization of a random variable $Y$. Lumber is heterogenous relative to many other engineering materials, even within the same nominal grade, and characterizing the induced uncertainty is important. The reliability problem is to estimate either:

• •

The failure probability $p_{f}=\textrm{Pr}(Y>y_{f})$, for a given maximum allowable deflection $y_{f}$; or

• •

The quantile $y_{f}$ such that $\textrm{Pr}(Y>y_{f})=p_{f}$, where $p_{f}$ is a given small probability.

3 Gaussian Processes and Parameter Estimation

3.1 Gaussian Process Model

We now briefly review the GP model. For the analysis of computer experiments Sacks et al. (1989) treated a deterministic output function $y({\bf x})$ as a realization of a GP,

[TABLE]

Here, $\mu$ is the mean and $Z({\bf x})$ is a stationary GP with zero mean and unknown variance $\sigma^{2}$. The mean could also be a regression function in the inputs ${\bf x}$, but Chen et al. (2016) showed that a constant mean model is often sufficient. Given two input vectors ${\bf x}$ and ${\bf x}^{\prime}$, a power-exponential structure is assumed for the correlation between $Z({\bf x})$ and $Z({\bf x}^{\prime})$:

[TABLE]

where $\theta_{j}>0$ and $1\leq p_{j}\leq 2$ for all $j$.

Conditional on the training data ${\bf y}$, the correlation parameters ${\boldsymbol{\psi}}=(\theta_{1},\ldots,\theta_{d},p_{1},\ldots,p_{d})$, $\mu$ and $\sigma^{2}$, the predictive distribution for $y({\bf x}^{*})$ at a new point ${\bf x}^{*}$ is Gaussian, i.e.,

[TABLE]

with predictive mean

[TABLE]

and predictive variance

[TABLE]

Here, ${\bf r}({\bf x}^{*})$ is an $n\times 1$ vector obtained from (3) with element $i$ given by $R({\bf x}^{*},{\bf x}^{(i)})$ for $i=1,\ldots,n$, the $n\times n$ matrix ${\bf R}$ has element $i,j$ given by $R({\bf x}^{(i)},{\bf x}^{(j)})$ for $1\leq i,j\leq n$, and ${\bf 1}$ is an $n\times 1$ vector of 1’s.

3.2 Parameter Estimation

The quantification of the total uncertainty of prediction is critical for the search criteria in Section 4. Hence, a Bayesian method accounting for parameter estimation uncertainty (Chen et al., 2017) is used and briefly described here.

At any stage the computer model has been evaluated at $n$ input configurations ${\bf x}^{(i)}$ for $i=1,\ldots,n$. With a vague uniform prior on $\mu$ and an inverse-gamma prior on $\sigma^{2}$ with both the shape and scale parameters tending to zero, the predictive distribution given all correlation parameters ${\boldsymbol{\psi}}$ becomes

[TABLE]

i.e., a shifted and scaled $t$ distribution with $n-1$ degrees of freedom (Handcock and Stein, 1993). The predictive mean in (5) becomes

[TABLE]

where

[TABLE]

and the predictive variance in (6) becomes

[TABLE]

where

[TABLE]

Next, the correlation parameters ${\boldsymbol{\psi}}$ are handled by Markov chain Monte Carlo (MCMC). The MCMC algorithm takes independent priors on the smoothness parameters $p_{j}$ as well as the $\theta_{j}$ as described by Chen et al. (2017) in their fully Bayesian implementation of the power-exponential structure. Let ${\boldsymbol{\psi}}^{(i)}$ for $i=1,\ldots,M$ be the MCMC sample of size $M$ of the correlation parameters. For each ${\boldsymbol{\psi}}^{(i)}$, the conditional predictive mean and variance in (8) and (9), respectively, are computed. The overall predictive mean is the average of these means,

[TABLE]

and the overall predictive variance is obtained through the law of total variance,

[TABLE]

These Bayesian methods are used to compute the predictive mean and variance for all the sequential-search criteria in Section 4.

4 Sequential Experimental Design

Here we describe the sequential design strategies for estimation of a tail probability or quantile. The heart of a sequential algorithm is the criterion for adding new evaluations to the search, and we contrast existing criteria.

4.1 Sequential Algorithms

The basic idea is to apply MC using the GP predictions of $y({\bf x}_{\text{MC}}^{(i)})$ for all points in the MC set ${\bf\mathcal{X}_{\text{MC}}}$ rather than evaluating the computer model. Given initial training data of $n_{0}$ evaluations, ${\bf y}=y({\bf x}^{(1)}),\ldots,y({\bf x}^{(n_{0})})$, and hence a trained GP, Algorithms 1 and 2 compute tail probability and quantile estimates, respectively. They add $n_{+}$ points sequentially according to the search criteria in Sections 4.2 and 4.3.

Note that in Section 6, the probability estimate in step 6 of Algorithm 1 will be replaced by a stratified random sampling estimate. Step 6 in Algorithm 2 for quantile estimation will be replaced in an analogous way. Both algorithms are based on the method of Ranjan et al. (2008) for contour estimation up to the search criterion. Mapping out a given contour where $y({\bf x})=y_{f}$ is equivalent to separating the MC points into those where $y({\bf x}_{\text{MC}}^{(i)})>y_{f}$ versus those where $y({\bf x}_{\text{MC}}^{(i)})\leq y_{f}$. In Algorithm 1 it is straightforward to estimate $p_{f}$ (Bichon et al., 2009) with $\hat{y}({\bf x}_{\text{MC}}^{(i)})$ equal to $\hat{m}({\bf x}_{\text{MC}}^{(i)})$ in (10). For quantile estimation, Algorithm 2 estimates the unknown quantile at each iteration and then chooses the next point as if the estimated quantile is the known contour of interest (Roy and Notz, 2014)

4.2 Expected Improvement Criterion

For mapping out a contour where $y({\bf x})$ equals the constant $y_{f}$, the improvement function proposed by Ranjan et al. (2008) can be written as

[TABLE]

where $v_{{\boldsymbol{\psi}},\sigma^{2}}({\bf x})$ is given by

[TABLE]

For the derivation of EI the predictive distribution $y({\bf x})$ is taken to be normal, with known predictive variance. Hence, the two subscripts ${\boldsymbol{\psi}}$ and $\sigma^{2}$ emphasize the predictive variance is conditional on these two GP parameters. The constant $\alpha$ determines the level of confidence. Thus, a large improvement would result from a new evaluation at ${\bf x}$ where the uncertainty measured by $v_{{\boldsymbol{\psi}},\sigma^{2}}({\bf x})$ is currently large and $y({\bf x})$ turns out to be close to the target $y_{f}$. Taking the expectation of $I({\bf x})$ with respect to the predictive distribution of $y({\bf x})$ gives EI,

[TABLE]

where $\hat{m}_{{\boldsymbol{\psi}}}$ is the predictive mean in (8), $u_{1}=(y_{f}-\hat{m}_{{\boldsymbol{\psi}}}({\bf x}))/\sqrt{v_{{\boldsymbol{\psi}},\sigma^{2}}({\bf x})}-\alpha$, $u_{2}=(y_{f}-\hat{m}_{{\boldsymbol{\psi}}}({\bf x}))/\sqrt{v_{{\boldsymbol{\psi}},\sigma^{2}}({\bf x})}+\alpha$, and $\phi(\cdot)$ and $\Phi(\cdot)$ are the PDF and CDF of the standard normal distribution, respectively. In practice, the conditional mean $\hat{m}_{{\boldsymbol{\psi}}}({\bf x})$ and variance $v_{{\boldsymbol{\psi}},\sigma^{2}}({\bf x})$ have to be replaced by their unconditional estimates in (10) and (11), respectively. The sequential-design criterion for selecting the next point, ${\bf x}^{*}$, based on EI is

[TABLE]

where ${\bf\mathcal{X}_{\text{cand}}}$ can be either a discrete or continuous candidate set.

4.3 Discrepancy Criterion

Roy and Notz (2014) defined the “discrepancy” between $y({\bf x})$ and $y_{f}$ at untried input vector ${\bf x}$ to be

[TABLE]

where $\epsilon>0$. If $\epsilon=0$, the expression looks like an $F$-statistic. Roy and Notz (2014) proposed choosing the next code evaluation at

[TABLE]

where the expectation is respect to the predictive distribution of $y({\bf x})$. The expectation can be written as

[TABLE]

and, as we work with $\epsilon=0$, an equivalent optimization is

[TABLE]

where $v_{{\boldsymbol{\psi}},\sigma^{2}}({\bf x})>0$ if the search omits points already in the training set. Another equivalent formulation is

[TABLE]

where $Z\sim N(0,1)$. We see that the discrepancy criterion is related to EI in (13), but it does not require a confidence-level $\alpha$. Like the maximization of EI in (14), maximizing (17) balances choosing a point with predictive mean close to $y_{f}$ (local search) and choosing a point with large predictive variance (global search).

5 Example: Short Column Function

To illustrate the sequential criteria and to start to address some questions about their practical implementation, we take a model of a short column with uncertain material properties under an uncertain load. It was used by Kuschel and Rackwitz (1997) to study the trade-off between structural reliability and cost; we consider only the reliability component of the function. For a given input ${\bf x}=(x_{m},x_{p},x_{z})$ it is

[TABLE]

where the output $y({\bf x})$ is the limit state, i.e., the difference between resistance and load, and $y({\bf x})<0$ means the system fails. The three inputs ($x_{m},x_{p},x_{z}$) are sampled from their respective independent distributions (Kuschel and Rackwitz, 1997):

• •

$X_{m}\sim N(\text{mean}=2000,\text{sd}=400)$ is the bending moment;

• •

$X_{p}\sim N(\text{mean}=500,\text{sd}=100)$ is the axial force;

• •

$X_{z}\sim\text{Lognormal}(\text{mean}=5,\text{sd}=0.5)$ is the yield stress.

The parameters $b$ and $h$ are the width and depth of the cross-section, respectively, both in mm. Because the function is trivial to compute in this illustration, it is possible to simulate $10$ million points from the input distributions and establish that taking $b=3$ and $h=10$ gives $\textrm{Pr}(y({\bf X})<0)=0.0025$ with negligible binomial standard error. Therefore, the true probability of system failure is $0.0025$, which we aim to estimate using sequential methodology.

Note that the probability of interest in the short column function is actually a lower-tail probability. Hence step 6 of Algorithm 1 will change to $\hat{p}_{f}=(1/N)\sum_{i=1}^{N}\mathds{1}(\hat{y}({\bf x}_{\text{MC}}^{(i)})<0)$. The quantile estimation in Algorithm 2 will change in a similar way.

5.1 Probability Estimation

Suppose the total experimental budget is $n=40$. The size of the initial design is kept as $n_{0}=20$, but two choices are considered for the type of design. The first option simply generates 20 points at random from the independent input distributions. The second choice starts with a random Latin hypercube design (LHD, McKay et al., 1979) in three variables on $[0,1]^{3}$ and then maps them uniformly onto the ranges given by the mean $\pm 3$ standard deviations of the actual distributions. Hence, the 20 values of the first input, $x_{m}$, are uniform on $[2000-3\times 400,2000+3\times 400]$; similarly the second input. The values of the lognormal third input just require exponentiating after analogous operations. The two types of initial design are called “random” and “uniform”.

We can see from (18) that large $x_{m}$, large $x_{p}$, and small $x_{z}$ tend to lead to failure. Hence, we speculate that the uniform initial design will lead to better performance in estimation of the failure probability, as it over-samples the tails of each distribution where failures occur. We suspect that over-sampling of extremes and combinations of extremes in the training set is desirable for many applications, to identify failure modes.

The MC set is $100,000$ points independently generated at random from the input distributions. The optimization in (14) or (17) is done by additionally generating $10,000$ points from the input distributions to form a finite candidate set and adding the point that maximizes EI or the discrepancy criterion in the candidate set. The whole process is repeated 10 times. These sample sizes are meaningfully chosen such that they keep a balance between the estimation accuracy and the computational time.

Figure 3 shows the results starting from random initial designs (the same 10 initial designs are used for the EI and discrepancy criteria). Estimates from the EI search barely show any sign of convergence. In contrast, after adding a further $11$ points, the estimates based on the discrepancy search criterion essentially converge to the true probability. They stay there and do not diverge for all 10 repeats of the experiment. The root mean squared error (RMSE) over the 10 repeats summarized in the first row of Table 1 confirms that the discrepancy criterion performs better than EI for the random initial design. Further calculations, not reported here for brevity, show that the EI method requires about 30 extra points after the 20-point initial design, to converge.

Figure 4 tells a similar story for uniform initial designs: the discrepancy criterion outperforms EI. The second row of Table 1 confirms this. In addition, if we compare Figures 3 and 4, the uniform initial designs perform better than the random initial designs, in agreement with our conjecture. Taking a uniform initial design is especially helpful for the EI criterion.

The superior convergence using the discrepancy criterion relative to EI is next explored by examining the points chosen by the two criteria. Visualization is facilitated by defining $t_{1}=x_{m}/x_{z}$ and $t_{2}=(x_{p}/x_{z})^{2}$ and noting that the function in (18) can be rewritten as

[TABLE]

Thus the failure boundary is a function of just two variables, $t_{1}$ and $t_{2}$, and is linear in them. For one of the repeat experiments adding 20 points, Figure 5 shows the (common) initial points in the $(t_{1},t_{2})$ subspace, the failure boundary of interest, and the $20$ extra points. We observe that the initial points are far from the contour of interest, but the sequential methods explore the space and, in particular, the contour. Among the $20$ points added, however, the discrepancy criterion places $13$ points near the failure boundary, whereas the EI method chooses the majority of the points in less relevant regions.

5.2 Quantile Estimation

Following the same settings as in Section 5.1, we investigate extreme quantile estimation. Results are reported in Figure 6 for random initial designs, where we see the discrepancy criterion again providing much faster convergence, similarly in Figure 7 for initial uniform designs. While comparison of Figures 6 and 7 shows that the initial uniform designs perform poorly at the outset and when just a few points are added sequentially, the RMSE values in Table 2 indicate that the uniform designs ultimately provide a better solution for both search criteria at $n=40$.

Thus, the results for probability estimation and quantile estimation lead to the same conclusions. The discrepancy criterion gives estimates that converge faster than those for the EI method for both initial designs. A uniform initial design tends to outperform a random initial design. Based on these observations, we recommend the discrepancy criterion starting from a uniform initial design.

6 Application to a Computer Model of a Floor System

6.1 Preliminary Analysis

Consider the computer model introduced in Section 2 with $d=8$. Before doing any formal modelling, we carry out a preliminary sensitivity analysis by fitting a constant mean Gaussian process with training data from an initial random LHD with $n=20$. We use functional analysis of variance (ANOVA) which takes each variable’s main effect as well as all of the two level interaction effects into account (Schonlau and Welch, 2006). All higher order interactions are ignored. The ANOVA decomposition results are reported in Table 3.

From Table 3, we observe that the four middle beams have the most important effects on the predicted response, $y$. The relationship between them and the response is illustrated in Figure 8, again using the methods of Schonlau and Welch (2006). The relationship between the four side beams and $y$ are similar to those of the four middle beams but weaker. It is clear from Figure 8 that a small MOE value results in a larger deflection. Therefore, the lower tail of the input distribution is important when estimating the upper probability. This is valuable prior information obtained from the preliminary analysis.

In the short column function of Section 5, the function is explicitly specified in (18), and we can determine immediately the critical part of the input space. In practice, however, when the functional relationship between inputs and output is unknown, we can carry out a preliminary study as in this section to examine the effects of the inputs.

6.2 Modelling the Input Distribution

We need to specify the input distribution. From the analysis in Section 6.1, we know that a small MOE leads to a large deflection. Therefore, the lower tail of each marginal input distribution is critical. A dataset of $580$ MOE values measured from $580$ boards collected from production is available to us. One way to proceed is to fit a parametric input distribution, for example the Weibull. However, training a parametric distribution based on all of the $580$ MOE data values does not emphasize the importance of the lower tail. Therefore, initially we considered the following two semi-parametric input distributions:

• •

a mixture distribution, with probability $p_{1}=0.1$ to sample from a 2-parameter censored Weibull distribution and $p_{2}=0.9$ to sample with replacement from the upper 90% of the empirical distribution of the available data. The 2-parameter censored Weibull distribution is trained using complete data from the lower 10% of the available data with the rest (90% of the data) treated as right censored (Liu et al., 2018).

• •

As above but using a 3-parameter Weibull distribution for the first stratum.

The unknown parameters were estimated by MLE. The resulting estimated density curves for the 2- and 3-parameter Weibull distributions are virtually identical and added to the empirical histogram of the lower 10% data shown in Figure 9.

It is clear that there is a non-negligible bump (3 boards) in the lower tail of the dataset, which neither of the semi-parametric distributions is able to characterize.

Hence, we propose to use instead the following non-parametric input distribution, $H(x)$:

[TABLE]

where $G_{1}(x)$ is the empirical distribution of the lower 10% of the MOE data, and $G_{2}(x)$ is the empirical distribution of the upper 90% of the data. The weights $p_{1}=p_{2}=0.5$ deliberately over-sample the lower tail of the input distribution. They will later be corrected using methods for weighted stratified random sampling when computing the probability/quantile estimates. In a nutshell, we use the empirical distribution of the 580 MOE data values as the input distribution, but deliberately over-sample the lower tail when generating the MC set and the finite candidate set.

6.3 MC Set

The MC set is comprised of $12,800$ points generated as follows. Each dimension has two strata for the mixture distribution: either sampling from $G_{1}(x)$ or sampling from $G_{2}(x)$. Hence, in total there are $2^{8}=256$ strata. For each of the $256$ strata, we generate $50$ different points from the relevant combination of $G_{1}(x)$ and $G_{2}(x)$ distributions and thus we have $256\times 50=12,800$ points in total to form the MC set. The above sample sizes are meaningfully chosen such that they keep a balance between the approximation accuracy and the computational time.

Working with this MC set, step 6 in Algorithm 1 for probability estimation will change to

[TABLE]

where ${\bf x}_{\text{MC}}^{(hi)})$ is MC sample point $i$ in stratum $h$, and the stratum weights $w_{h}$ sum to 1. For instance, the combination with all inputs coming from $G_{1}(x)$ has $w_{h}=(0.1)^{8}$ as the stratum weight. For quantile estimation, step 6 in Algorithm 2 also changes to finding $\hat{y}_{f}$ such that $\hat{p}_{f}$ equals a pre-specified probability, $a$. In practice, any trial value of $y_{f}$ gives a $\hat{p}_{f}$ in (20). We minimize $|\hat{p}_{f}-a|$ numerically with respect to $y_{f}$ using the optimize() function in the MASS package of R to find $\hat{y}_{f}$ with $10^{-6}$ tolerance.

6.4 True Probability and Quantile

With $p_{f}=0.001$, we simulated 10 different MC sets from the mixture distribution in (19) and ran the computer model. With 10 repeats the mean of the estimates of the $0.999$ quantile is $3.88057$ inches and the standard error of the mean is $0.00701$. Therefore, the “true” $0.999$ quantile is taken to be 3.88 inches.

6.5 Application Results

Suppose the design budget is $n=30$. The number of points for the initial design is $n_{0}=20$ and $10$ additional points are chosen sequentially by optimizing a search criterion. The whole process is repeated $10$ times with different initial designs. The MC set contains 12,800 points simulated as described in Section 6.3. The finite candidate set is the same as the MC set. For the initial design with only $n_{0}=20$ points, however, we take a random uniform LHD to cover the input space globally, following the recommendations from the short column model in Section 5.

Results based on the discrepancy criterion are reported in Figure 10.

The initial estimates of the tail probability are typically only about 30% of the true value of 0.001 but vary tightly around the truth after 10 additional points, The quantile estimates similarly converge to the true quantile. Moreover, the extreme probability and quantile are estimated well using just the carefully designed 12,800 points in the MC and candidate sets and adjusted stratum weights. This computational efficiency is a result of the preliminary analysis that established a negative association between the MOEs of the joists and the associated deflection.

7 Diagnostics

In practice, the true tail probability or quantile is unknown, yet a user still wants guidance about whether the algorithm has converged. Furthermore, the question of convergence relates to the specific single experiment actually executed. Here we describe a diagnostic for convergence of a sequential search based on the discrepancy criterion.

From the formulation in (17), the negative of the absolute discrepancy is $-|\hat{m}_{{\boldsymbol{\psi}}}({\bf x})-y_{f}|/\sqrt{v_{{\boldsymbol{\psi}},\sigma^{2}}({\bf x})}$. Ideally, it will become more negative for all points in the MC set as the training set increases. Thus, we compute this criterion point-wise for the MC set at each step and draw a boxplot against the sample size. We would expect a declining trend with sample size.

Figures 11 and 12 show diagnostic plots for one repeat of the experiments with the short column function and the floor system computer model, respectively. Both experiments start from a uniform initial design.

Consistent with our expectation, both plots exhibit a decreasing trend. After adding 10 points to the initial design for the short column function, most of the discrepancies in the MC set are less than $-10$, i.e., a huge standardized distance between the predictive mean and the contour of interest. The pattern carries over to the computer model for the floor system, which again indicates the algorithm is converging.

8 Concluding Remarks

This paper illustrates the usefulness of a sequential strategy to estimate an extreme probability or its associated quantile. Throughout the paper, we have addressed some practical issues an engineer faces when using a sequential design. From the analyses, we have the following recommendations.

• •

Use the discrepancy criterion, which is more straightforward and converges faster than EI in our study.

• •

Use a uniform initial design, to over-sample the tails of the input distributions (for the floor model only one tail is over-sampled).

• •

It may be more efficient to choose the MC set and the candidate set according to a stratified weighting scheme, to enrich the sets with points in the failure region. That is more efficient than generating huge sets from the input distributions.

• •

It can be important to do a preliminary analysis to gain some insights on the relationship between inputs and output.

• •

The search criteria considered in this paper use the standard error of prediction to guide local/global search. As Bayesian methods can provide more realistic uncertainty estimates (Chen et al., 2017), their use is recommended for sequential search.

About the authors

\censor

Hao Chen was a PhD candidate in the Department of Statistics at the time of this research. He is currently a senior research scientist at R&D, LoyaltyOne Inc.

\censor

William J. Welch is a Professor in the Department of Statistics, University of British Columbia.

Funding

\censor

Research supported by the Natural Sciences and Engineering Research Council, Canada and Compute Canada/WestGrid.