Global Fitting of the Response Surface via Estimating Multiple Contours of a Simulator

TL;DR

This paper introduces two innovative sequential design methods utilizing contour estimation to improve response surface fitting for computationally expensive simulators, enhancing prediction accuracy for risk assessment.

Contribution

The paper presents novel sequential design strategies based on contour estimation to optimize input selection for better response surface approximation.

Findings

Proposed methods outperform traditional designs in accuracy

Numerical examples demonstrate improved prediction quality

Approaches are effective for complex physical system simulations

Abstract

Computer simulators are nowadays widely used to understand complex physical systems in many areas such as aerospace, renewable energy, climate modeling, and manufacturing. One fundamental issue in the study of computer simulators is known as experimental design, that is, how to select the input settings where the computer simulator is run and the corresponding response is collected. Extra care should be taken in the selection process because computer simulators can be computationally expensive to run. The selection shall acknowledge and achieve the goal of the analysis. This article focuses on the goal of producing more accurate prediction which is important for risk assessment and decision making. We propose two new methods of design approaches that sequentially select input settings to achieve this goal. The approaches make novel applications of simultaneous and sequential contour estimations. Numerical examples are employed to demonstrate the effectiveness of the proposed approaches.

Equations

$$y(\textbf{x}_{i})=\textbf{f}\mbox{\boldmath${\beta}$}+Z(\textbf{x}_{i}),i=1,2,...,n,$$

$$R({\text{x}}_i,{\text{x}}_j)$$

$$\hat{\mu }=({\text{1}}_n^{\text{mbox}T}{\text{R}}^{-1}{\text{1}}_n{)}^{-1}{\text{1}}_n^{\text{mbox}T}{\text{R}}^{-1}\text{y}$$

$${\hat{\sigma }}^2=\frac{(\text{y}-{\text{1}}_n\hat{\mu }{)}^T{\text{R}}^{-1}(\text{y}-{\text{1}}_n\hat{\mu })}{n}.$$

$$\hat{y}({\text{x}}_0)=E[y({\text{x}}_0)|\text{y}]=\mu +{\text{r}}^T({\text{x}}_0){\text{R}}^{-1}(\text{y}-\mu {\text{1}}_n),$$

$$s^2(x_0)=Var(y({\text{x}}_0)|\text{y})={\sigma }^2\left(1-{\text{r}}^T({\text{x}}_0){\text{R}}^{-1}\text{r}({\text{x}}_0)\right).$$

$$E(I(\text{x}))=(\hat{y}(\text{x})-y({\text{x}}_{j^{\ast }}){)}^2+Var(\hat{y}(\text{x}))=(\hat{y}(\text{x})-y({\text{x}}_{j^{\ast }}){)}^2+s^2(\text{x})$$

$$I(\text{x})=(\hat{y}(\text{x})-y({\text{x}}_{j^{\ast }}){)}^2$$

$${GE}_p={\left\{\sum _{i=1}^{n-1}\sum _{j=i+1}^n{\left(\frac{q({\text{x}}_i)q({\text{x}}_j)}{d({\text{x}}_i,{\text{x}}_j)}\right)}^p\right\}}^{1/p},$$

$${\text{x}}_{n+1}=\underset{{\text{x}}_0\in \chi }{\mathrm{argmin}}\sum _{i=1}^n{\left(\frac{\hat{q}({\text{x}}_i)\hat{q}({\text{x}}_0)}{d({\text{x}}_i,{\text{x}}_0)}\right)}^p.$$

$$S(a)=\{\text{x}\in \chi :y(\text{x})=a\}.$$

$$I(\textbf{x})=\epsilon^{2}(\textbf{x})-\min\big{\{}(y(\textbf{x})-a)^{2},\epsilon^{2}(\textbf{x})\big{\}},$$

$$E[I(\text{x})]$$

$$I(\text{x})={ϵ}^2(\text{x})-\min \{(y(\text{x})-a_1{)}^2,\dots ,(y(\text{x})-a_k{)}^2,{ϵ}^2(\text{x})\},$$

$$I(\text{x})$$

$$I(\textbf{x})=\left\{\begin{array}[]{ll}\epsilon^{2}(\textbf{x})-(y(\textbf{x})-a_{1})^{2},&\ \ \ a_{1}-\epsilon(\textbf{x})\leq y(\textbf{x})\leq\min\{a_{1}+\epsilon(\textbf{x}),(a_{1}+a_{2})/2\}\\ \ldots&\\ &\ \ \ \max\{a_{j}-\epsilon(\textbf{x}),(a_{j-1}+a_{j})/2\}\leq y(\textbf{x}),\\ \epsilon^{2}(\textbf{x})-(y(\textbf{x})-a_{j})^{2},&\ \ \ y(\textbf{x})\leq\min\{a_{j}+\epsilon(\textbf{x}),(a_{j}+a_{j+1})/2\},\\ &\ \ \ 2\leq j\leq k-1;\\ \ldots&\\ \epsilon^{2}(\textbf{x})-(y(\textbf{x})-a_{k})^{2},&\ \ \ \max\{a_{k}-\epsilon(\textbf{x}),(a_{k-1}+a_{k})/2\}\leq y(\textbf{x})\leq a_{k}+\epsilon(\textbf{x}),\\ 0,&\ \ \ \hbox{otherwise}.\end{array}\right.$$

$$v_{j1}(\textbf{x})=\left\{\begin{array}[]{ll}a_{1}-\epsilon(\textbf{x}),&j=1;\\ \max\{a_{j}-\epsilon(\textbf{x}),(a_{j-1}+a_{j})/2\},&2\leq j\leq k,\\ \end{array}\right.$$

$$v_{j2}(\textbf{x})=\left\{\begin{array}[]{ll}\min\{a_{j}+\epsilon(\textbf{x}),(a_{j}+a_{j+1})/2\},&1\leq j\leq k-1;\\ a_{k}+\epsilon(\textbf{x}),&j=k.\\ \end{array}\right.$$

$$E[I(\text{x})]$$

$$y=\frac{sin(10\pi x)}{2x}+(x-1{)}^4,0.5\le x\le \mathrm{2.5.}$$

$$y(\text{x})=[1+(4x_1+4x_2+1)](3+192x_1{x}_2),0\le x_1\le 1,0\le x_2\le 1.$$

$${\text{x}}_{opt}^{\ast }=\underset{1\le i\le m}{\mathrm{argmax}}{s}^2({\text{x}}_i^{\ast }).$$

$$RMSPE=\sqrt{\frac{1}{|{\mathcal{X}}_{pred}|}\sum _{\text{x}\in {\mathcal{X}}_{pred}}(\hat{y}(\text{x})-y(\text{x}){)}^2},$$

$$Maximum\; error=\underset{\text{x}\in {\mathcal{X}}_{pred}}{\max }|\hat{y}(\text{x})-y(\text{x})|.$$

$$y$$

$$y=\prod _{i=1}^3{x}_i,0\le x_i\le i,for,i=1,\dots ,3.$$

$$y=x_1{x}_2+x_3^2{x}_4^2,-1\le x_i\le 1,for,i=1,\dots ,4.$$

Full text

Global Fitting of the Response Surface via Estimating Multiple Contours of a Simulator

F. Yang†, *C. Devon Lin‡*111Address for correspondence: C. Devon Lin, Associate Professor, Department of Mathematics and Statistics, Queen’s University, Kingston, ON Canada (E-mail: devon.lin@queensu.ca)., and P. Ranjan§**

*†*College of Mathematics, Sichuan University, China

*‡*Department of Mathematics and Statistics, Queen’s University, Canada

§OM&QT Area, Indian Institute of Management Indore, India

Abstract Computer simulators are nowadays widely used to understand complex physical systems in many areas such as aerospace, renewable energy, climate modeling, and manufacturing. One fundamental issue in the study of computer simulators is known as experimental design, that is, how to select the input settings where the computer simulator is run and the corresponding response is collected. Extra care should be taken in the selection process because computer simulators can be computationally expensive to run. The selection shall acknowledge and achieve the goal of the analysis. This article focuses on the goal of producing more accurate prediction which is important for risk assessment and decision making. We propose two new methods of design approaches that sequentially select input settings to achieve this goal. The approaches make novel applications of simultaneous and sequential contour estimations. Numerical examples are employed to demonstrate the effectiveness of the proposed approaches.

Keywords: Computer experiment; Contour estimation; Gaussian process; Latin hypercube; Maximin design; Sequential design; Space-filling.

1 Introduction

Computer models or simulators are increasingly becoming popular for gaining insights of the physical processes and phenomena that are too expensive or infeasible to observe. For example, Greenberg (1979) developed a finite volume community ocean model (FVCOM) for simulating the flow of water in the Bay of Fundy; Bower et al. (2006) discussed the formation of galaxies using a simulator called GALFORM; and Bayarri et al. (2009) used a simulator called TITAN2D for modelling the maximum volcanic eruption flow height. Realistic computer simulators of complex processes can also be computationally expensive to run, and thus statistical surrogates trained on a handful of simulator runs are often used for the deeper understanding of the underlying phenomena. Sacks et al. (1989) proposed using a realization of the Gaussian process (GP) model as a surrogate for such processes.

The popular objectives of such computer experiments include global fitting, variable screening, and estimation of process features like the maximum, a pre-specified contour or a tail quantile region. Assuming the simulator under consideration is expensive to run, the number of simulator runs would be limited and thus one must be careful in choosing the inputs. Over the last two decades several innovative methodologies and algorithms have been developed to address some of the concerns. See Santner, Williams and Notz (2003), Fang, Li and Sudjianto (2005) and Rasmussen and Williams (2006) for details.

We focus on efficient designs for global fitting. In computer experiments literature, a popular technique is to use Latin hypercube designs (McKay et al., 1979) with some space-filling properties like maximin interpoint distance (Johnson et al., 1990, Morris and Mitchell, 1995), minimum pairwise coordinate correlation (Iman and Conover, 1982, Joseph and Hung, 2008), orthogonal array-based structure (Owen, 1992, Tang, 1993), projection property (Joseph, Gul and Ba, 2015), etc. Such designs aim at filling the input space as evenly as possible, but do not consider the complexity of the response surface. On the other hand, D-optimal designs (Johnson et al., 1990), integrated mean squared prediction error (IMSPE)-optimal designs (Sacks, Schiller and Welch, 1989) and maximum mean squared prediction error (MMSPE)-optimal designs (Sacks and Schiller, 1988) use the process response information in finding a design for global fitting.

Most of these designs follow one-shot approach, i.e., all design points are obtained at the same time. However, over the past decade, a few sequential designs have also been proposed for global fitting of the response surface that have higher prediction accuracy. For instance, the D-optimal design (Gramacy and Lee, 2009), expected improvement (EI) criterion based design (Lam and Notz, 2008) and minimum potential energy based design (Joseph et al., 2015). In this paper, we propose two new sequential design approaches for global fitting.

The main idea behind our proposed approaches comes from the fact that the estimation of a response surface can be approximated by the estimation of a large number of contours over the range of the responses. This further motivated us to generalize the contour estimation idea (Ranjan et al., 2008) for our objective. In this paper, we propose two generalizations. First, we recommend splitting the range of simulator outputs into $k$ equi-spaced contours and then develop a new EI criterion for the simultaneous estimation of these pre-specified multiple contours. Second, we propose a new adaptive approach of choosing contour levels for selecting the follow-up trial by maximizing the EI criterion for contour estimation. The performance of the proposed approaches have been compared with several state of the art designs for global fitting.

The remainder of the article is organized as follows. Section 2 presents a quick review of the Gaussian process (GP) model for building a surrogate of the computer model output, popular sequential design approaches for global fitting (Lam and Notz, 2008, Joseph et al., 2015) and the EI criterion for contour estimation (Ranjan et al., 2008). Section 3 presents the new multiple contours estimation-based EI method for constructing designs for global fitting of the response surface. In Section 4, we propose the new adaptive method of estimating the contour levels for choosing follow-up design points in the sequential framework. The performance comparison of the proposed methods and the existing approaches are discussed in Section 5. Finally, Section 6 summarizes the key findings and concluding remarks.

2 Background Review

This section reviews the necessary background and the existing relevant work for later development. More specifically, we provide a brief account of reviews on Gaussian process models used throughout, the existing sequential design approaches for global fitting as well as the contour estimation in Ranjan et al. (2008). Although these topics can easily be accessed in the literature, we include them here so that this article is a standalone document.

2.1 Gaussian Process Models

Gaussian process models are most widely used in computer experiments to emulate outputs from computer codes (e.g., Sacks et al., 1989). Its popularity is due to its simplicity, flexibility and the ability of providing the predictive uncertainty. Here we cover the key concepts of GP models and refer the reader to Santner, Williams and Notz (2003) and Rasmussen and Williams (2006) for details. For a training data of size $n$, let the $i$th input and output of a computer code be a $d$-dimensional vector $\textbf{x}_{i}=(x_{i1},\ldots,x_{id})$ and a scalar $y_{i}=y(\textbf{x}_{i})$, for $i=1,\ldots,n$. Typically, without the loss of generality, the design domain is assumed to be a unit hypercube, $\chi=(0,1)^{d}$. A GP model assumes

[TABLE]

where f is a vector of regression functions, ${\beta}$ is the vector of regression parameters, $Z(\textbf{x})$ is a stationary stochastic process with mean zero, constant variance $\sigma^{2}$, and the correlation between two outputs $y(\textbf{x}_{i})$ and $y(\textbf{x}_{j})$ being denoted by $R(\textbf{x}_{i},\textbf{x}_{j})=\hbox{corr}(\textbf{x}_{i},\textbf{x}_{j})$. In this article, we focus on the Gaussian process models with a constant mean, that is, \textbf{f}\mbox{\boldmath{\beta}}=\mu. Let $\textbf{y}=(y_{1},\ldots,y_{n})^{T}$ be the vector of responses for the training data and R be an $n\times n$ spatial correlation matrix with the $(i,j)$th element $R(\textbf{x}_{i},\textbf{x}_{j})$. A GP model in (1) is equivalent to assume that y follows a multivariate normal distribution with mean vector $\mu\textbf{1}_{n}$ and the covariance matrix $\sigma^{2}\textbf{R}$ with $\textbf{R}=(\textbf{R}_{ij})$, where $\textbf{1}_{n}$ is an $n$-dimensional column vector of all 1’s. Notationally, we denote $\textbf{y}\sim GP(\mu\textbf{1}_{n},\sigma^{2}\textbf{R})$. There are many choices of valid correlation functions. One popular choice is the Gaussian correlation function,

[TABLE]

where $\theta_{k}$ is the correlation parameter for the $k$th input variable. The unknown parameters in the model include the mean $\mu$, the variance $\sigma^{2}$, and $d$ correlation parameters $\theta_{1},\ldots,\theta_{d}$. They can be estimated via the maximum likelihood approach or Bayesian approach such as Markov Chain Monte Carlo (MCMC) (Santner et al., 2003, Fang et al., 2005, Currin et al., 1988, Linkletter et al., 2006). For the maximum likelihood approach, if the correlation parameters are known, the estimates of $\mu$ and $\sigma^{2}$ in (1) are

[TABLE]

and

[TABLE]

The best linear unbiased predictor (BLUP) at an input $\textbf{x}_{0}$ is given by

[TABLE]

where $\textbf{r}(\textbf{x}_{0})=(R(\textbf{x}_{0},\textbf{x}_{1}),\ldots,R(\textbf{x}_{0},\textbf{x}_{n}))^{\mbox{\tiny T}}$. Moreover, the predictive variance of $y(\textbf{x})$ is

[TABLE]

In practice, the unknown correlation parameters in (3) and (4) are replaced with the estimates. Thus, $\mu$, $\sigma^{2}$, R and $\textbf{r}(\textbf{x}_{0})$ in (5) and (6) are replaced by $\hat{\mu}$, $\hat{\sigma}^{2}$, $\hat{\textbf{R}}$ and $\hat{\textbf{r}}(\textbf{x}_{0})$, respectively. There are a number of R packages that can provide the GP model fitting. They include mlegp, GPfit, DiceKriging, tgp, RobustGaSP and SAVE (Dancik, 2018, MacDonald, Ranjan and Chipman, 2015, Roustant, Ginsbourger and Deville, 2018, Gramacy and Taddy, 2016, Gu, Palomo and Berger, 2018, Palomo, Paulo and Garcia-Donato, 2015). These R packages are different in terms of computational efficiency and stability. In general they shall provide similar results. For the reason of stability, we use the R package GPfit (MacDonald et al., 2015) in this article.

2.2 Existing Sequential Design Approaches for Global Fitting

The general setup of a sequential design approach starts with an initial design and adds one point or a batch of points at-a-time sequentially. We focus on the sequential approaches of adding one point at-a-time. The next follow-up point shall be chosen based on the information gathered from the existing data and shall be most informative among the candidate points. The process of adding points is repeated until a tolerance based stopping criterion is met or a pre-specified budget is exhausted. The step-by-step process of the sequential design approach is as follows.

• Step 1.

Choose an initial design of run size $n_{0}$. Let $n=n_{0}$.

• Step 2.

Build a statistical surrogate model using the available data $\{(\textbf{x}_{i},y_{i}),i=1,\ldots,n\}$.

• Step 3.

Choose the next design point $\textbf{x}_{n+1}$ based on a criterion. Run the computer code at the new input $\textbf{x}_{n+1}$ and obtain the corresponding response $y_{n+1}$.

• Step 4.

Let $n=n+1$ and repeat Steps 2 and 3 until it reaches the run size budget or satisfies the stopping criterion.

A few remarks are in order. First, the initial design typically comes with some space-filling property like maximin interpoint distance, minimum pairwise coordinate correlation, etc. If the initial run size $n_{0}$ is too small, the resulting surrogate model could be wildly inaccurate and mislead the follow-up design choice. On the other hand, if the $n_{0}$ is relatively large, it may not fully take the advantage of sequential design criterion in Step 3. The notion of expected improvement (EI) criterion has become extremely popular for choosing follow-up design points after Jones et al. (1998) developed an EI criterion for finding the global minimum of a computer simulator response. One recommendation, given by Ranjan et al. (2008), for the value of $n_{0}$ is $25-35\%$ of the ultimate run size budget. Such a recommendation is based on their sequential design approach for contour estimation. Second, the run size budget certainly depends on the computer code of interest. Loeppky et al. (2009) provided a rule of thumb for selecting a sample size, that is, 10 times the number of input variables. In our illustrative examples, the total runsize is at least $10d$. Third, in principle, any modelling methods such as GP, treed GP (TGP), or Bayesian additive regression trees (BART) (Gramacy and Lee, 2008, Chipman et al., 2012) can be used as a surrogate in Step 2. We focus on GP modelling in the examples.

2.2.1 Expected improvement criterion by Lam and Notz (2008)

Lam and Notz (2008) introduced a sequential design approach based on an expected improvement for global fit (EIGF) criterion which chooses the next input point that maximizes the expected improvement

[TABLE]

where the improvement function $I(\textbf{x})$ is defined as

[TABLE]

with $y(\textbf{x}_{j^{*}})$ being the observed output at the sampled point, $\textbf{x}_{j^{*}}$, that is closest in distance to the candidate point x. They use the Euclidean distance to determine this nearest sampled design point. The expectation in (7) is taken with respect to the predictive distribution of $y(\textbf{x})$ under the GP model, i.e., $y(\textbf{x})\sim N(\hat{y}(\textbf{x}),s^{2}(\textbf{x}))$. The EIGF criterion in (7) balances the local search and global search of the next potential design input that guides the search for the “informative” regions with significant variation in the response values.

2.2.2 Sequential minimum energy designs by Joseph et al. (2015)

Motivated by the fact in physics that the charged particles in a box repel and try to remain away from each other as much as possible, Joseph et al. (2015) viewed a space-filling design in the experimental region as the positions occupied by the charged particles in a box. The charge of each particle represents the experimental response. A minimum energy design is obtained by minimizing the potential energy. Let $q(\textbf{x})$ be the charge of the particle at the design input x and $d(\textbf{x}_{i},\textbf{x}_{j})$ denote the Euclidean distance between the $i$th and the $j$th input. Joseph et al. (2015) defined the potential energy of a design $\textbf{D}_{n}=\{\textbf{x}_{1},\ldots,\textbf{x}_{n}\}$ as

[TABLE]

where $p$ is in the range of $[1,\infty)$. They further proposed a sequential minimum energy design approach which works as follows. Let $\hat{q}(\textbf{x})=\{\hat{y}(\textbf{x})\}^{-1/(2d)}$, where $d$ is the dimensionality of the input x. Then the proposed one point at-a-time greedy algorithm finds the next follow-up design point given by

[TABLE]

The design generated by this algorithm is called sequential minimum energy design (SMED).

2.3 Contour Estimation via EI Criterion

The contour at level “$a$” of a simulator response surface consists of all the inputs x that yield the same response $a$, that is,

[TABLE]

Ranjan et al. (2008) developed an expected improvement criterion under the sequential design methodology for estimating a contour from an expensive to evaluate computer simulator with scalar responses. The proposed improvement function is,

[TABLE]

where $y(\textbf{x})$ has a normal predictive distribution, i.e., $y(\textbf{x})\sim N(\hat{y}(\textbf{x}),s^{2}(\textbf{x}))$, and $\epsilon(\textbf{x})=\alpha s(\textbf{x})$ for a positive constant $\alpha$. A suggested value for $\alpha$ is 1.96 for the reason that this value defines a region of interest around $S(a)$ to be 95% confidence interval under the normality assumption of the responses. Letting $v_{1}(\textbf{x})=a-\epsilon(\textbf{x})$ and $v_{2}(\textbf{x})=a+\epsilon(\textbf{x})$, the closed form of the expectation of the improvement function $I(\textbf{x})$ with respect to the predictive distribution of $y(\textbf{x})$ is given by,

[TABLE]

where $u_{1}=[v_{1}(\textbf{x})-\hat{y}(\textbf{x})]/s(\textbf{x})$, $u_{2}=[v_{2}(\textbf{x})-\hat{y}(\textbf{x})]/s(\textbf{x})$, and $\phi(\cdot)$ and $\Phi(\cdot)$ are the probability density function and the cumulative distribution function of a standard normal random variable, respectively. See Ranjan et al. (2008) and the associated Errata for the derivation of (12). The first term in (12) suggests an input with a large $s(\textbf{x})$ in the neighbourhood of the predicted contour, while the last term assigns the weights to points that are far away from the predicted contour with large uncertainties. The second term is often dominated by the other two terms in (12). Maximizing the EI criterion in (12) results in the inputs with high uncertainty near the predicted contour as well as those far away, achieving both aims of local search and global exploration.

3 Global Fitting by Estimating Multiple Contours

This section proposes a new method for constructing a sequential design for achieving higher prediction accuracy of the overall global fit. The basic sequential framework would remain the same as in Section 2.2, that is, start with a good initial design (e.g., maximin Latin hypercube) of size $n_{0}\ll n$ and then sequentially add the remaining $n-n_{0}$ points using some method that feeds on the objective of global fitting. Instead of the conventional approach of trying to evenly fill the input space, the proposed idea is to slice the response surface into multiple contours and then use the sequential design approach to simultaneously estimate those contours. Next, we generalize the EI criterion for contour estimation (Ranjan et al., 2008) for simultaneous estimation of multiple contours.

For a given integer $k>0$ and the set of scalar values $a_{1},\ldots,a_{k}\in[y_{min},y_{max}]$, suppose that we are interested in estimating $k$ contours $S(a_{1}),\ldots,S(a_{k})$, where $[y_{min},y_{max}]$ represents the range of the true simulator response, and $S(\cdot)$ is defined in (10). Without loss of generality, assume $a_{1}<a_{2}<\cdots<a_{k}$. For choosing the follow-up trial, we propose the improvement function at input x as

[TABLE]

where $y(\textbf{x})\sim N(\hat{y}(\textbf{x}),s^{2}(\textbf{x}))$ and $\epsilon(\textbf{x})=\alpha s(\textbf{x})$ for some positive constant $\alpha$. This improvement function will be non-zero only if $(y(\textbf{x})-a_{j})^{2}<\epsilon^{2}(\textbf{x})$ for some $j$. Therefore, the improvement function can be re-written as:

[TABLE]

Since $a_{1}<a_{2}<\cdots<a_{k}$, the improvement function can be further simplified as

[TABLE]

The term $\epsilon(\textbf{x})$ defines an uncertainty band around each contour that is a function of the predictive standard deviation $s(\textbf{x})$. For the design points already chosen, the radius of the band is exactly zero. In addition, the criterion will tend to be large for the samples from the sets $(\{\textbf{x}:y(\textbf{x})=a_{1}\},\{\textbf{x}:y(\textbf{x})=a_{2}\},\ldots,\{\textbf{x}:y(\textbf{x})=a_{k}\})$, where $s(\textbf{x})$ is large.

Similar to other sequential design approaches, we suggest choosing follow-up design points by maximizing the corresponding expected improvement, where the expectation is taken with respect to the predictive distribution, $y(\textbf{x})\sim N(\hat{y}(\textbf{x}),s^{2}(\textbf{x}))$. For $j=1,\ldots,k$, let $v_{j1}(\textbf{x})$’s and $v_{j2}(\textbf{x})$’s be defined as follows,

[TABLE]

and

[TABLE]

Then, the expectation of the improvement function in (13) is simply the sum of the individual contour estimation EI criterion of Ranjan et al. (2008) over $k$ cases, i.e.,

[TABLE]

where $u_{j1}=(v_{j1}(\textbf{x})-\hat{y}(\textbf{x}))/s(\textbf{x})$ and $u_{j2}=(v_{j2}(\textbf{x})-\hat{y}(\textbf{x}))/s(\textbf{x})$, $\phi(\cdot)$ and $\Phi(\cdot)$ are the probability density function and the cumulative distribution function of a standard normal random variable, respectively. The formulation in (3) reduces to (12) when the number of contour levels is $k=1$.

Note that the maximization of $E[I(\textbf{x})]$ over $I(\textbf{x})$ has two advantages. First, the true value of $y(\textbf{x})$ (and hence $I(\textbf{x})$) is unknown for any unsampled design point. Second, some regions of the design space may not have been sufficiently explored yet and the predictive variance of $\hat{y}(\textbf{x})$ is relatively high. For such an unsampled design point, even though the predicted response is not within the $\epsilon(\textbf{x})$-band of one of those $k$ contours, the true contours may lie in the unexplored region. As a result, the EI approach facilitates a balance between the local exploitation versus global exploration.

A naive way to choose $k$ - the number of contours, and contour levels, $a_{1},a_{2},\ldots,a_{k}$, for global fitting is to use equi-spaced $k$ contours in the simulator output range $[y_{min},y_{max}]$, and finding their optimal values appear to be a challenging task. We now present two illustrations of the proposed multiple contour estimation EI criterion (referred to as MC criterion) for global fitting with different values of $k$.

Example 1. Consider the computer model (Gramacy and Lee, 2012) that relates the one-dimensional input $x$ and the output $y$ as,

[TABLE]

The true relationship between the input $x$ and the output $y$ is displayed in the blue solid curve in Figure 1. Five initial design points are shown by black empty circles. We then sequentially add 15 design points using the MC criterion in (3). The numerical labels represent the order of the newly added design points. Figures 1(a), 1(b), and 1(c) illustrate the sequential design scheme with the MC criterion for 5, 20, 50 equally spaced contour levels within the ranges of the fitted surface.

Clearly Figure 1 reveals that the proposed MC criterion-based sequential design approach can choose the inputs that are around the areas where the function changes the direction and locate most of the points in the areas where the computer model is more complex. The results show that the final design is in general space-filling but with some nearby points, for example, points 10 and 15, in Figure 1(c).

Example 2. Consider a computer model with two-dimensional input variables $\textbf{x}=(x_{1},x_{2})$, and the output given by

[TABLE]

Suppose a maximin Latin hypercube design of 10 points is generated and the corresponding responses are collected from the computer model. First, we consider searching for the next follow-up design point for estimating only one contour at the level $a=300$. Figure 2 shows these 10 design points, the inputs with $I(\textbf{x})>0$, and the maximizer of the EI criterion for contour estimation in (12) for the candidate set on a regular $100\times 100$ rectangular grid.

Next, we consider the simultaneous estimation of three contours at levels $a_{1}=150$, $a_{2}=300$ and $a_{3}=600$ using the MC criterion in (3). Figure 3 shows the inputs from the same $100\times 100$ grid candidate set that achieve non-zero improvement (13), and the point that maximizes the MC criterion. This point in red in Figure 3 is in fact from the set of the points that yield non-zero improvement around the contour level $a_{1}=150$.

Figure 4 illustrates the complete sequential design scheme with $10$ initial design points and $30$ follow-up design points for simultaneously estimating three contours at levels $a_{1}=150$, $a_{2}=300$ and $a_{3}=600$. The red points are the new follow-up points and the label corresponds to the order the point is added. The last panel displays the squared distance between the estimated contour and the true contour at each stage. It can be observed from Figure 4 that the estimated uncertainty bands around the three contours become narrower and more accurate. It can also be seen that more points are added to estimate the contour level $a_{1}=150$ than to estimate the other two contours. Some points such as the second and third points are away from the contour bands.

It is clear from the two examples that the resulting designs do not have the conventional space-filling property. This is desirable as the objective is an overall good fit of the response surface, and not to explore the input space. However, as illustrated in Figure 4, a significant fraction of design points tend to line up on the pre-specified contours, which could lead to biased designs if $a_{1},...,a_{k}$ are not chosen appropriately. Next we propose an efficient method of selecting contour levels.

4 Sequential Estimation of Contours for Global Fitting

In this section we propose a new approach for choosing the follow-up design points. Different from the previous section where the simultaneous estimation of multiple contours was used for global fitting, we adopt the EI criterion for estimating only one contour level at each stage. More importantly, we propose using different contour levels at different stages. That is, we sequentially find the design points that maximize the criterion (12) with a different value of $a$ at each stage. The important issue is how to choose the contour level at each stage. We propose the following way to choose such contour level in an automatic way.

Suppose at stage $j$, the training data are $\{(\textbf{x}_{i},y_{i}),i=1,\ldots,n)\}$ and the corresponding emulator gives the predictive distribution as $y(\textbf{x})\sim N(\hat{y}(\textbf{x}),s^{2}(\textbf{x}))$ for any input x. Let the candidate set for the next follow-up point be $\textbf{x}^{*}_{1},\ldots,\textbf{x}^{*}_{m}$ and

[TABLE]

Then, we choose the contour level at stage $j$ as $a_{j}=\hat{y}(\textbf{x}^{*}_{opt})$. In other words, at each stage, we set the contour level to be the fitted response that has maximum predictive variance. This is to encourage exploring the area with maximum uncertainty.

Example 2 (contd.) Consider finding a design for global fitting of the computer simulator in Example 2. The procedure starts with an initial design of size $n_{0}=10$ obtained via maximin Latin hypercube sampling and $n-n_{0}=30$ follow-up points are chosen as per the proposed sequential strategy. Figure 5 displays the follow-up design points found by the proposed method, i.e., the sequential contour estimation-based EI criterion.

Note that the resulting design is more randomized as compared to a systematic layout of points on the contour lines shown in Figure 4. Again, the design is not completely space-filling and it has some pairs of close-by points.

5 Simulated Examples

In this section, we conduct a simulation study to demonstrate the effectiveness of the proposed sequential design approaches. Specifically, we compare the proposed approaches with the following methods:

• (a)

a one-shot maximin Latin hypercube design;

• (b)

the sequential $D$-optimal design in the R package tgp;

• (c)

the sequential approach by Lam and Notz (2008);

• (d)

the sequential minimum energy design in Joseph et al. (2015);

• (e)

the proposed multiple contours estimation-based criterion in Section 3;

• (f)

the proposed sequential contour estimation-based criterion in Section 4.

For approach (e), we use the 10 contour levels, that is, $k$ is set to be 10. These methods are denoted by ‘maximinLHD’, ‘tgp’,‘EIGF’, ‘SMED’, ‘MC$\_$10’, and ‘SC$\_$var’.

Several criteria can be used to evaluate the performance of different design approaches in comparison. We adopt the root mean square prediction error (RMSPE) given by

[TABLE]

where $\hat{y}(\textbf{x})$ and $y(\textbf{x})$ are the predicted response and the true response at the new input x in the hold-out set $\mathcal{X}_{pred}$. Another criterion we use is the maximum error provided by

[TABLE]

For each example below, the initial design for sequential designs is a maximin Latin hypercube design of $n_{0}$ runs generated using the R package SLHD (Ba, 2015). The model fitting is implemented using the default setting of the function GP_fit in the R package GPfit. The test data is a random Latin hypercube design of $1000d$ points where $d$ is the number of input variables. The parameter $\alpha$ in ‘MC$\_$10’ and ‘SC$\_$var’ is set to be 2.

Example 3. We consider computer model with two input variables $x_{1}$ and $x_{2}$,

[TABLE]

This model is known as Branin function (Dixon and Szego, 1978). We use $n_{0}=10$ initial design points. The total run size budget is 30. Figure 6 displays the boxplots of RMSPEs and maximum errors of the different design approaches over 50 simulations. The results show that in this example the one-shot approach ‘maximinLHD’ is the worst while the approaches ‘MC$\_$10’ and ‘SC$\_$var’ are comparably better than the others.

Example 4. We consider computer model with three input variables $x_{1}$, $x_{2}$ and $x_{3}$,

[TABLE]

We use $n_{0}=20$ initial design points. The total run size budget is 60. Figure 7 displays the boxplots of RMSPEs and maximum errors of the different design approaches over 50 simulations. Again, here the one-shot approach ‘maximinLHD’ is the worst. The approaches ‘tgp’, ‘MC$\_$10’ and ‘SC$\_$var’ are comparably better than the others in terms of RMSPE and the proposed approaches ‘MC$\_$10’ and ‘SC$\_$var’ are significantly better than the others in terms of maximum errors.

Example 5. We consider computer model with four input variables $x_{1}$, $x_{2}$, $x_{3}$ and $x_{4}$,

[TABLE]

We use $n_{0}=27$ initial design points. The total run size budget is 80. Figure 8 displays the boxplots of RMSPEs and maximum errors of the different design approaches over 50 simulations. Here, the approach ‘EIGF’ is the worst followed by the approach ‘maximinLHD’ in terms of both criteria. The performance of the other four approaches are similar based on RMSPE. However, based on the maximum error, the proposed approaches give more accurate predictions.

6 Concluding Remark

In this article we have developed two sequential design approaches for accurately predicting a complex computer code. The approaches are based on the expected improvement criteria for simultaneously or sequentially estimating contours. We used a Gaussian process (GP) model as a surrogate for the computer simulator, which is an integral component of the proposed criteria for identifying the follow-up trials. Numerical examples are given to demonstrate that the proposed approaches can significantly outperform the existing approaches.

Note that if some other surrogate is used instead of GP model, then also the key ideas like formulation of improvement function and sequential estimation of contour levels can be retained. Of course, the resultant expected improvement criteria would change, and in fact, one may not even end up with a closed form expression of the final design criterion for selecting follow-up points. Future work also include the application of the proposed contour estimation-based sequential design approaches for global fitting for computer experiments with both qualitative and quantitative factors (Deng et al., 2017) and dynamic computer experiments (Zhang, Lin and Ranjan, 2018).