%2multibyte Version: 5.50.0.2960 CodePage: 1252
%\documentclass{article}%
\documentclass[ijoc,nonblindrev]{informs3} % current default for manuscript submission 
\OneAndAHalfSpacedXII % current default line spacing


% Natbib setup for author-year style
\usepackage{natbib}
 \bibpunct[, ]{(}{)}{,}{a}{}{,}%
 \def\bibfont{\small}%
 \def\bibsep{\smallskipamount}%
 \def\bibhang{24pt}%
 \def\newblock{\ }%
 \def\BIBand{and}%



\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{xcolor}
\usepackage{graphicx}%
\setcounter{MaxMatrixCols}{30}
%TCIDATA{OutputFilter=latex2.dll}
%TCIDATA{Version=5.50.0.2960}
%TCIDATA{Codepage=1252}
%TCIDATA{CSTFile=40 LaTeX article.cst}
%TCIDATA{Created=Friday, March 16, 2018 12:23:05}
%TCIDATA{LastRevised=Friday, February 28, 2020 12:46:46}
%TCIDATA{<META NAME="GraphicsSave" CONTENT="32">}
%TCIDATA{<META NAME="SaveForMode" CONTENT="1">}
%TCIDATA{BibliographyScheme=Manual}
%TCIDATA{<META NAME="DocumentShell" CONTENT="Standard LaTeX\Blank - Standard LaTeX Article">}
%TCIDATA{Language=American English}
%BeginMSIPreambleData
\providecommand{\U}[1]{\protect\rule{.1in}{.1in}}

\usepackage{graphicx}

%EndMSIPreambleData
\newtheorem{theorem}{Theorem}
\newtheorem{acknowledgement}[theorem]{Acknowledgement}
\newtheorem{algorithm}[theorem]{Algorithm}
\newtheorem{axiom}[theorem]{Axiom}
\newtheorem{case}[theorem]{Case}
\newtheorem{claim}[theorem]{Claim}
\newtheorem{conclusion}[theorem]{Conclusion}
\newtheorem{condition}[theorem]{Condition}
\newtheorem{conjecture}[theorem]{Conjecture}
\newtheorem{corollary}[theorem]{Corollary}
\newtheorem{criterion}[theorem]{Criterion}
\newtheorem{definition}[theorem]{Definition}
\newtheorem{example}[theorem]{Example}
\newtheorem{exercise}[theorem]{Exercise}
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{notation}[theorem]{Notation}
\newtheorem{problem}[theorem]{Problem}
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{remark}[theorem]{Remark}
\newtheorem{solution}[theorem]{Solution}
\newtheorem{summary}[theorem]{Summary}
\newenvironment{proof}[1][Proof]{\noindent\textbf{#1.} }{\ \rule{0.5em}{0.5em}}

\MANUSCRIPTNO{JOC-2019-05-OA-114} 
\begin{document}


% Author's names for the running heads
% Sample depending on the number of authors;
% \RUNAUTHOR{Jones}
\RUNAUTHOR{Kleijnen and Van Beers}
% \RUNAUTHOR{Jones, Miller, and Wilson}
% \RUNAUTHOR{Jones et al.} % for four or more authors
% Enter authors following the given pattern:
%\RUNAUTHOR{}

% Title or shortened title suitable for running heads. Sample:
% \RUNTITLE{Bundling Information Goods of Decreasing Value}
% Enter the (shortened) title:
\RUNTITLE{LOO-CV Supplement}

\KEYWORDS{validation, cross-validation, Kriging, Gaussian process, Monte Carlo}

% Full title. Sample:
% \TITLE{Bundling Information Goods of Decreasing Value}
% Enter the full title:


\TITLE{Online supplement of \textquotedblleft Statistical tests for cross-validation of Kriging models\textquotedblright}
%\TITLE{Online supplement of "Statistical tests for cross-validation of %Kriging models"}
\maketitle

% Block of authors and their affiliations starts here:
% NOTE: Authors with same affiliation, if the order of authors allows, 
%   should be entered in ONE field, separated by a comma. 
%   \EMAIL field can be repeated if more than one author
\ARTICLEAUTHOR{%
\AUTHOR{Jack P.C. Kleijnen and Wim C.M. van Beers}
\AFF{Department of Management \\Tilburg School of Economics and Management (TiSEM) \\Tilburg University (TiU) \\Postbox 90153, 5000 LE Tilburg, Netherlands \\E-mail: kleijnen@tilburguniversity.edu\\ Orcid iD: https://orcid.org/0000-0001-8413-2366}

%%%\AUTHOR{Author2}
%%%\AFF{Author2 affiliation, \EMAIL{}, \URL{}}
% Enter all authors
} % end of the block



\date{}

\begin{abstract}


\end{abstract}


\textbf{Appendix 1: Basics of ordinary Kriging}\newline

We focus on \textit{ordinary Kriging }(OK) instead of \textit{universal
Kriging} (UK). We use the notation in Kleijnen (2020), except that we focus on
\textit{leave-one-out cross-validation }\textbf{(}LOO-CV) so we denote the
\textit{new }input combination to be predicted by $\mathbf{x}_{i}$ instead of
the usual $\mathbf{x}_{0}$, and we use the $n-1$ \textit{old} or
\textit{training} combinations $\mathbf{X}_{-i}$ instead of the $n$ old
combinations $\mathbf{X}$. OK\textit{\ }assumes the following (meta)model:%
\begin{equation}
y(\mathbf{x})=\mu+M(\mathbf{x}) \label{GP}%
\end{equation}
where $\mathbf{x}$ is a combination of the $d$ simulation inputs $x_{j}$ ($j$
= 1, ..., $d$), $\mu$ is the constant mean E$[y(\mathbf{x)]}$, and
$M(\mathbf{x})\mathbf{\ }$is a zero-mean stationary GP. We let $\mathbf{X}%
_{-i}$ denote the $(n-1)\times d$ matrix with the rows $\mathbf{x}_{i^{\prime
}}$ = $(x_{i^{\prime};1}$, ..., $x_{i^{\prime};d})$ with $i$ $\neq$
$i^{\prime} $. Furthermore, we use the symbol $\boldsymbol{\Sigma}_{M}$ =
$(\sigma_{i^{\prime};i^{\prime\prime}})$ = $\mathbf{(}$Cov$[y(\mathbf{x}%
_{i^{\prime}}\mathbf{)},y(\mathbf{x}_{i^{\prime\prime}}\mathbf{)}])$ (with
$i^{\prime},i^{^{\prime\prime}}$ $\neq$ $i$) to denote the $(n-1)\times(n-1)$
matrix with the covariances between the metamodel's $(n-1) $ old outputs.
Analogously, we let $\mathbf{\sigma}_{M}\mathbf{(x}_{i})$\textbf{\ = }%
$(\sigma_{i;i^{\prime}})$ = $\mathbf{(}$Cov$[y(\mathbf{x}_{i}\mathbf{)}%
,y(\mathbf{x}_{i^{\prime}}\mathbf{)}])$ denote the $(n-1)$-dimensional vector
with the covariances between the new output $y(\mathbf{x}_{i})$ and the
$(n-1)$ old outputs $y(\mathbf{x}_{i^{\prime}}\mathbf{)}$. We let
$\mathbf{1}_{n-1}$ denote the $(n-1)$-dimensional vector with all its elements
equal to $1$. Finally, we let $\mathbf{w}_{-i}$ denote the $(n-1)$-dimensional
vector with the observed simulation outputs for the $(n-1)$ non-deleted combinations.

If and only if we assume that OK is a \textit{valid} metamodel of the
underlying simulation model, then we may write $y$ = $w$. We can then derive
the\textit{\ best linear unbiased predictor }(BLUP) $\widehat{y}%
(\mathbf{x}_{i})$, using the $n-1$ non-deleted I/O combinations $(\mathbf{X}%
_{-i},\mathbf{w}_{-i})$. This BLUP is a weighted average of these
$\mathbf{w}_{-i}$:%
\begin{equation}
\widehat{y}(\mathbf{x}_{i})=\boldsymbol{\lambda}_{-i}^{\prime}\mathbf{w}_{-i}.
\label{Kriging estimate}%
\end{equation}
\bigskip If these weights\ $\boldsymbol{\lambda}_{-i}$ are selected optimally,
then\ (by definition) the resulting BLUP has minimum variance and is unbiased.
We can prove that this BLUP is%
\begin{equation}
\widehat{y}(\mathbf{x}_{i})=\mu+\mathbf{\sigma}_{M}\mathbf{(x}_{i})^{\prime
}\mathbf{\Sigma}_{M}^{-1}(\mathbf{w}_{-i}\mathbf{-}\mu\mathbf{1}_{n-1}).
\label{Krige predictor}%
\end{equation}
We note that this equation implies that if $\mathbf{w}_{-i}$
%TCIMACRO{\TEXTsymbol{>}}%
%BeginExpansion
$>$%
%EndExpansion
\textbf{\ }$\mu\mathbf{1}_{n-1}$ (or $(\mathbf{w}_{-i}\mathbf{-}\mu
\mathbf{1}_{n-1})$ $>$ $0$), then $\widehat{y}(\mathbf{x}_{i})$
%TCIMACRO{\TEXTsymbol{>} }%
%BeginExpansion
$>$
%EndExpansion
$\mu$. Furthermore, if we denote Var$(y)$ by $\tau^{2}$, then we can prove%
\begin{equation}
\text{MSPE}[\widehat{y}(\mathbf{x}_{i})]\text{ }\text{=}\text{ }\tau
^{2}-\mathbf{\sigma}_{M}(\mathbf{x}_{i})^{\prime}\mathbf{\Sigma}_{M}%
^{-1}\mathbf{\sigma}_{M}(\mathbf{x}_{i})+\frac{[1-\mathbf{1}_{n-1}^{\prime
}\boldsymbol{\Sigma}_{M}^{-1}\mathbf{\sigma}_{M}(\mathbf{x}_{i})]^{2}%
}{\mathbf{1}_{n-1}^{\prime}\mathbf{\Sigma}_{M}^{-1}\mathbf{1}_{n-1}}.
\label{MPSE2}%
\end{equation}
This MSPE\textbf{\ }equals Var$[\widehat{y}(\mathbf{x}_{i})]$ if the metamodel
is valid so $\widehat{y}(\mathbf{x}_{i})$ is unbiased.

Sometimes, it is convenient to switch from covariances to
\textit{correlations}; i.e., sometimes we switch to the correlation matrix
$\mathbf{R}$ = $(\rho_{i^{\prime};i^{\prime^{\prime}}})$, which equals
$\tau^{-2}\mathbf{\Sigma}_{M}$; analogously, $\boldsymbol{\rho}(\mathbf{x}%
_{i})$ = $\tau^{-2}\boldsymbol{\sigma}_{M}(\mathbf{x}_{i})$. There are several
types of correlation functions; see (e.g.) Rasmussen and Williams (2006, pp.
80--104). In simulation, however, the most popular function is the
\textit{separable}, \textit{anisotropic}\textbf{\ }(instead of isotropic)
\textit{Gaussian correlation function}%
\begin{equation}
\rho(\mathbf{h,}\boldsymbol{\theta})\text{ = }%
%TCIMACRO{\tprod \limits_{j=1}^{d}}%
%BeginExpansion
{\textstyle\prod\limits_{j=1}^{d}}
%EndExpansion
\exp\left(  -\theta_{j}h_{j}^{2}\right)  =\exp(-\sum_{j=1}^{d}\theta_{j}%
h_{j}^{2})\text{ with }\theta_{j}\geq0. \label{Krige cor funct}%
\end{equation}
with the \textit{correlation factors }$\boldsymbol{\theta}$ = $(\theta_{j}) $
and the \textit{distance} vector $\mathbf{h}$ = $(h_{j})$ where $h_{j}$ =
$\left\vert x_{g;j}-x_{g^{\prime};j}\right\vert $ and $g,g^{\prime}$ = 1, ...,
$n$. So, $\mathbf{R}$ = $\mathbf{R(h,}\boldsymbol{\theta}\mathbf{)}$. If
$\theta_{j}$ $\downarrow$ $0$, then $\exp$ $\left(  -\theta_{j}h_{j}%
^{2}\right)  $ $\uparrow$ 1 for any distance $h_{j}$; i.e., input $j $ has a
strong effect for any $h_{j}$. If $\theta_{j}$ $\uparrow$ $\infty$, then
$\exp$ $\left(  -\theta_{j}h_{j}^{2}\right)  $ $\downarrow$ 0, so input $j$
has virtually\textbf{\ }no effect. We use the symbol $\boldsymbol{\psi}$ to
denote the vector with the $(2+d)$ \textit{Kriging (hyper)parameters}
$(\mu,\tau^{2},\theta_{1},...,\theta_{d})^{\prime}$.

In practice, we must estimate these (nuisance) OK parameters $\boldsymbol{\psi
}$. The most popular estimator is the \textit{maximum likelihood estimator
}(MLE), which we denote by $\widehat{\boldsymbol{\psi}}$ =
$\widehat{\boldsymbol{\psi}}(\mathbf{X},\mathbf{w})$. This MLE is the solution
of%
\begin{equation}
\min_{\mathbf{\psi}}\ln[\left\vert \tau^{2}\mathbf{R}\right\vert
]+(\mathbf{\mathbf{w}}_{-i}\mathbf{-}\mu\mathbf{\mathbf{1}}_{n-1}%
\mathbf{)}^{\prime}(\tau^{2}\mathbf{R)}^{-1}(\mathbf{w}_{-i}\mathbf{-}%
\mu\mathbf{\mathbf{1}}_{n-1}\mathbf{)}\text{ with }\mathbf{\theta}%
\geq\mathbf{0} \label{MLE}%
\end{equation}
\bigskip where $\mathbf{R}$ = $\mathbf{R(h,}\boldsymbol{\theta}\mathbf{)}$ and
$\left\vert \mathbf{R}\right\vert $ denotes the determinant of $\mathbf{R}$.
This MLE implies the following explicit formula for $\widehat{\mu}$:%
\begin{equation}
\widehat{\mu}=(\mathbf{1}_{n-1}^{\prime}\mathbf{R}^{-1}\mathbf{1}_{n-1}%
)^{-1}\mathbf{1}_{n-1}^{\prime}\mathbf{R}^{-1}\mathbf{\mathbf{w}_{-i},}
\label{mhu^}%
\end{equation}
which is a weighted mean of the simulation outputs $\mathbf{w}_{-i}$.
Furthermore, this MLE implies
\begin{equation}
\widehat{\tau}^{2}=\frac{1}{n}(\mathbf{w}_{-i}-\widehat{\mu}\mathbf{1}%
_{n-1})^{\prime}\mathbf{R}^{-1}(\mathbf{w}_{-i}-\widehat{\mu}\mathbf{1}%
_{n-1}), \label{tau^}%
\end{equation}
so $\widehat{\tau}^{2}$ = $\widehat{\tau}^{2}(\widehat{\mu},\mathbf{R})$ is a
weighted mean of the squared \textit{residuals} $(\mathbf{w}_{-i}%
-\widehat{\mu}\mathbf{1}_{n-1})$. However, (\ref{MLE}) does not give an
explicit formula for $\widehat{\boldsymbol{\theta}}$, but requires an
iterative search. Actually, solving (\ref{MLE}) is a mathematical challenge;
e.g., different $\widehat{\boldsymbol{\psi}}$ may result from different
software packages or from different starting values for the same package; see
Erickson et al. (2018) and our MC experiments.

To obtain $\widehat{y}(\mathbf{x}_{i},\widehat{\boldsymbol{\psi}}_{-i})$ and
$s^{2}[\widehat{y}(\mathbf{x}_{i},\widehat{\boldsymbol{\psi}}_{-i})]$, we
\textit{plug} the estimator $\widehat{\boldsymbol{\psi}}_{-i}$ into
(\ref{Krige predictor}) and (\ref{MPSE2}). Most frequentist publications
ignore the consequences of this plugging-in; i.e., they ignore the fact that
$\widehat{y}(\mathbf{x}_{i},\widehat{\boldsymbol{\psi}}_{-i})\mathbf{\ }%
$becomes nonlinear and $s^{2}[\widehat{y}(\mathbf{x}_{i}%
,\widehat{\boldsymbol{\psi}}_{-i})]$ underestimates the true Kriging
variance. Our bootstrap algorithm (in the main text) gives an
unbiased variance estimator.

If we did not re-estimate $\boldsymbol{\psi}$ (resulting in
$\widehat{\boldsymbol{\psi}}_{-i}$), then we could apply a
computational shortcut to compute\textbf{ }$\widehat{y}(\mathbf{x}%
_{0},\widehat{\boldsymbol{\psi}})$ and\textbf{ }$s^{2}[\widehat{y}%
(\mathbf{x}_{0},\widehat{\boldsymbol{\psi}})]$; see Ginsbourger
and\ Sch\"{a}rer (2021) and Sch\"{o}bi et al. (2015).

Note: Bayesian approaches combine the prior distribution on $\boldsymbol{\psi
}$\ with the simulation I/O data that give $\widehat{\boldsymbol{\psi}} $, to
obtain the posterior distribution using Bayes theorem. To estimate this
posterior distribution, Bayesians apply \textit{Markov chain MC }(MCMC)
algorithms; see Efron (2015). Our bootstrap algorithm is also a type of
MC.\newline

\textbf{Appendix 2: Literature review of cross-validation in Kriging}\newline

We focus on \textit{recent} publications on \textit{cross-validation}
(CV), so we do not discuss an old milestone publication such as Stone (1974).
\textit{Leave-one-out CV }(LOO-CV) is discussed in\ Gramacy (2020) and Santner
et al. (2018); however, those two textbooks do not cover the test
that we derive in this paper. Bartz-Beielstein (2016) discusses $k$-fold CV
(if $k$ = 1, then this equals LOO-CV), to select an \textquotedblleft
ensemble\textquotedblleft\ (combination) of different types of metamodels
including Kriging metamodels. Kamath and Fan (2018) uses LOO-CV and $k$-fold
CV to compare several types of metamodels including Kriging models. Afendras
and Markatou (2019), Antunes et al. (2019), Fuhg et al. (2020), and
Vi\'{u}dez-Moreiras (2019) also discuss $k$-fold CV. Xiao et al. (2018) uses
$k$-fold CV (including LOO-CV) in sequential sampling for reliability
problems, focusing on a so-called learning function to select the next point
to be simulated; this function is not a measure of the metamodel accuracy (see
Xiao et al., p. 409). Shu et al. (2017) also uses LOO-CV to select the next
point to be simulated; that publication measures metamodel accuracy through
either the \textit{root mean squared prediction error} (RMSPE) or the
\textit{maximum absolute error} (MAE) instead of our LOO-CV statistic; Garbo
and German (2019) also uses LOO-CV with the MSPE (which is the square of the
RMSPE) to select the next point. Rasmussen and Williams (2006) discusses
LOO-CV within a Bayesian framework. Bhosekar and Ierapetritou (2018) surveys
$k$-fold CV including LOO-CV for Kriging and other types of metamodels in
simulation. Both Strano et al. (2018) and Yin et al. (2018) apply
LOO-CV to select the best type of metamodel. Van Steenkiste et al. (2018) uses
$k$-fold CV for sequential \textit{sensitivity analysis }(SA) through
metamodels including Kriging metamodels, and the \textquotedblleft root
relative squared error (CV$_{\text{RRSE}}$)\textquotedblleft---which is
related to the coefficient of determination $R^{2}$---or the \textquotedblleft
Bayesian estimation error quotient\textquotedblleft; also see Gorissen et al.
(2010). Bacchi et al. (2018) discusses $k$-fold CV for Kriging metamodels,
including MSPE, $R^{2}$, and \textquotedblleft residual
analysis\textquotedblright\ where residuals are defined as $w_{i}%
-\widehat{y}_{i}$ ($i$ = 1, ..., $n$). Kleijnen (2015, pp. 114--121) discusses
LOO-CV, focussing on linear regression metamodels---but also presenting
references and websites for Kriging metamodels. We note that\textbf{ }Law
(2015) is the most popular textbook in discrete-event simulation, and covers
Kriging and linear regression metamodels---but does not discuss CV.

Concerning the \textit{analysis} of the CV results, most publications use
visual inspection of a scatterplot with the simulated outputs $w$ versus the
predicted outputs $\widehat{y}$; see Da Costa et al. (2018), Kamath and Fan
(2018), Kyprioti et al. (2020), Lupera Calahorrano et al. (2016), Quirante et
al. (2018), and Roustant et al. (2012). Moreover, Roustant et al. (2012, p.
35) uses the \textquotedblleft standardized residuals\textquotedblright%
---which are identical to our \textquotedblleft prediction error standardized
(PES)\textquotedblright---analyzed through the \textquotedblleft Normal
QQ-plot of standardized residuals\textquotedblright---instead of Bonferroni's
inequality---without re-estimating the Kriging parameters\textbf{.}
Parnianifard et al. (2018, p. 4) \textquotedblleft follow[s] Kleijnen
(2015)\textquotedblleft, eyeballs the scatterplot, computes the
\textquotedblleft standardized residuals\textquotedblright, and claims that
these residuals should fall in the range $[-3,3]$; however, that
standardization uses the RMSPE instead of our measure (Shu et al. (2017) also
uses the RMSPE, as we have already mentioned). Williams et al. (2019) uses the
Studentized prediction errors, but uses \textquotedblleft per
comparison\textquotedblright\ error rates instead of an \textquotedblleft
experimentwise\textquotedblright\ error rate. Many publications use the RMSPE,
which should be as small as possible; e.g., Forrester and Keane (2009)
considers a Kriging model with a RMSPE smaller than 2\% as a reasonably good
model. Ji et al. (2018) uses $R^{2}$, and seems to find $R^{2}$
%TCIMACRO{\TEXTsymbol{>} }%
%BeginExpansion
$>$
%EndExpansion
0.80 \ acceptable. Zhang et al. (2017) uses the \textquotedblleft correlation
coefficient\textquotedblright\ (but leaves that coefficient undefined; we
conjecture it is $R^{2}$) and a threshold (that\ seems subjective). De
Carvalho et al. (2017) applies LOC-CV to Kriging and radial basis function
(RBF) metamodels, quantifying prediction errors through $R^{2}$, the
\textquotedblleft relative average absolute error\textquotedblright\ (RAAE),
and the \textquotedblleft relative maximum absolute error\textquotedblright%
\ (RMAE); that article concludes: \textquotedblleft Most of the prediction
residuals were lower than 3\%\textquotedblleft. The review in Bhosekar and
Ierapetritou (2018) includes a table (namely, Table 4) with seven metrics for
CV; however, these metrics do not include our statistic. Fei et al. (2019)
lists six CV metrics, and gives qualitative thresholds (such as
\textquotedblleft close to 0\textquotedblright\ and \textquotedblleft close to
\'{1}\textquotedblright). Bastos and O'Hagan (2009) considers validation of
Kriging metamodels in a Bayesian framework, but uses a training set and a
validation set instead of CV; Jin and Jung (2016) uses Bastos and O'Hagan
(2009) and heuristic thresholds for $R^{2}$ and RMSPE. Both Luminari et al.
(2018, pp. 73--74) and Zhang et al. (2018)\textbf{ }compute the quadratic
difference between the Kriging predictor that uses all $n$ combinations and
the Kriging predictor that uses only $n-k$ combinations (where $k$ $\geq$ 1);
next Luminari et al. (2018) accepts a value below 6\%.

To analyze $k$-fold CV, Antunes et al. (2019) uses five well-known metrics;
namely, root mean squared error, mean absolute error (MAE), root relative
squared error, and Pearson's correlation coefficient. Demory et al. (2019)
also uses $R^{2}$, RMSPE, and the maximum absolute Studentized prediction
error---but not our test based on Bonferroni's inequality. Ginsbourger
and\ Sch\"{a}rer (2021) uses a Q-Q plot to analyze either\textbf{
}$w(\mathbf{x}_{i})-\widehat{y}(\mathbf{x}_{i},\widehat{\boldsymbol{\psi}%
})/s[\widehat{y}(\mathbf{x}_{i},\widehat{\boldsymbol{\psi}})]$ (using\textbf{
}$\widehat{\boldsymbol{\psi}}$, not $\widehat{\boldsymbol{\psi}}_{-i}$) or a
transformation of $w(\mathbf{x}_{i})-\widehat{y}(\mathbf{x}_{i}%
,\widehat{\boldsymbol{\psi}})$ that accounts for correlations among the $n$
prediction errors\textbf{ }$w(\mathbf{x}_{i})-\widehat{y}(\mathbf{x}%
_{i},\widehat{\boldsymbol{\psi}})$ ($i$ = 1, ..., $n$). Nahvi et al. (2019)
combines LOO-CV with a probabilistic divergence measure---namely, the
Kullback--Leibler divergence---to select the next simulation input combination
in multi-objective simulation-optimization---whereas we focus on prediction.
Mubarak et al. (2015) compares five Kriging variants, and uses the RMSPE in
$k$-fold CV to select a variant.

We note that the preceding literature review focuses on LOO-CV
in Kriging, but LOO-CV is also applied in many other modelling areas;
e.g., Lin et al. (2019) uses LOO-CV in simulation via $k$ nearest
neighbors. Our references demonstrate that CV and Kriging are applied in many
scientific disciplines.

Computational time for CV in Kriging may be a serious problem in case of large
data sets\ or \textquotedblleft big\textquotedblright\ $n$; see Dubrule
(1983).\ Kleijnen and Van Beers (2020) discusses big data in Kriging, but not
CV. Wang and Haaland (2019) also gives a (mathematically advanced) discussion
of numeric errors in Kriging.\newline

\textbf{Appendix 3: Latin hypercube sampling in Kriging}\newline

\textit{Latin hypercube sampling}\textbf{\ }(LHS) samples\textit{\ without
replacement}. If LHS samples midpoints, then LHS samples each midpoint only
once in the sample of size $n$. Altogether, $\mathbf{x}_{j}$ = $(x_{i;j})$
($i$ = 1, ..., $n$; $j$ = 1, ..., $d$) is a permutation of the $n$ midpoints;
i.e., $\mathbf{X}_{n\times d}$ is a stratified sample of $n$ points from a $d
$-dimensional grid.

Note: We claim that LHS with \textit{sampling} of points within subintervals
(instead of midpoints) has the disadvantage that it may give two values---in
two neighboring subintervals---that are very close together. So, the two
resulting outputs are close together (assuming a smooth I/O function, as
Kriging does) and give little new information (because these outputs have a
high positive correlation).\textbf{\ }However, sampling of midpoints ensures
that we obtain $n$ realizations of $x_{j}$ that are \textquotedblleft wide
apart\textquotedblright---which we conjecture gives better estimates of the
correlation parameters $\theta_{j}$ in the correlation function
(\ref{Krige cor funct}). Sampling midpoints never gives $x$ $\downarrow$ 0 or
$x$ $\uparrow$ 1 ($x$ = 0 or $x$ = 1 is impossible because $x$ is continuous;
$x$ = $\epsilon$ or $x$ = 1 - $\epsilon$ is also undesirable); we wish to
avoid these two extreme values because Kriging can use the output near $x$ = 0
only to predict the output at $x$
%TCIMACRO{\TEXTsymbol{>} }%
%BeginExpansion
$>$
%EndExpansion
0 (not at $x$
%TCIMACRO{\TEXTsymbol{<} }%
%BeginExpansion
$<$
%EndExpansion
0), so Kriging can use $x$ = 0 only in one direction; a similar argument holds
for $x$ $\uparrow$ 1.

LHS implies that projection of the $n$ points in $\mathbf{X}_{n\times d}$ onto
the $d$ individual axes gives $n$ \textit{non-collapsing} values per axis.
Consequently, it is relatively easy to estimate the $d$ individual separable
correlation functions in (\ref{Krige cor funct}).%

%TCIMACRO{\TeXButton{B}{\begin{table}[tbp] \centering}}%
%BeginExpansion
\begin{table}[tbp] \centering
%EndExpansion%
\begin{tabular}
[c]{|ccccccccccc|}\hline
$i$ & \multicolumn{1}{c|}{1} & \multicolumn{1}{c|}{2} & \multicolumn{1}{c|}{3}
& \multicolumn{1}{c|}{4} & \multicolumn{1}{c|}{5} & \multicolumn{1}{c|}{6} &
\multicolumn{1}{c|}{7} & \multicolumn{1}{c|}{8} & \multicolumn{1}{c|}{9} &
10\\\hline\hline
$x_{1}$ & \multicolumn{1}{||r}{-0.5} & \multicolumn{1}{|r}{-0.7} &
\multicolumn{1}{|r}{1.5} & \multicolumn{1}{|r}{1.7} &
\multicolumn{1}{|r}{-1.5} & \multicolumn{1}{|r}{-0.1} &
\multicolumn{1}{|r}{-0.3} & \multicolumn{1}{|r}{0.5} &
\multicolumn{1}{|r}{1.1} & \multicolumn{1}{|r|}{1.3}\\\hline
$x_{2}$ & \multicolumn{1}{||r}{0.9} & \multicolumn{1}{|r}{-0.7} &
\multicolumn{1}{|r}{0.1} & \multicolumn{1}{|r}{0.7} & \multicolumn{1}{|r}{1.5}
& \multicolumn{1}{|r}{1.9} & \multicolumn{1}{|r}{0.3} &
\multicolumn{1}{|r}{-1.9} & \multicolumn{1}{|r}{-1.7} &
\multicolumn{1}{|r|}{-0.5}\\\hline
\multicolumn{11}{c}{}\\\hline
$i$ & \multicolumn{1}{c|}{11} & \multicolumn{1}{c|}{12} &
\multicolumn{1}{c|}{13} & \multicolumn{1}{c|}{14} & \multicolumn{1}{c|}{15} &
\multicolumn{1}{c|}{16} & \multicolumn{1}{c|}{17} & \multicolumn{1}{c|}{18} &
\multicolumn{1}{c|}{19} & 20\\\hline\hline
$x_{1}$ & \multicolumn{1}{||r}{-1.1} & \multicolumn{1}{|r}{0.9} &
\multicolumn{1}{|r}{0.7} & \multicolumn{1}{|r}{-1.7} &
\multicolumn{1}{|r}{-1.9} & \multicolumn{1}{|r}{0.1} &
\multicolumn{1}{|r}{1.9} & \multicolumn{1}{|r}{-0.9} &
\multicolumn{1}{|r}{0.3} & \multicolumn{1}{|r|}{-1.3}\\\hline
$x_{2}$ & \multicolumn{1}{||r}{1.3} & \multicolumn{1}{|r}{-0.9} &
\multicolumn{1}{|r}{1.1} & \multicolumn{1}{|r}{-1.3} &
\multicolumn{1}{|r}{1.7} & \multicolumn{1}{|r}{-0.3} &
\multicolumn{1}{|r}{-1.5} & \multicolumn{1}{|r}{-1.1} &
\multicolumn{1}{|r}{0.5} & \multicolumn{1}{|r|}{-0.1}\\\hline
\end{tabular}
\caption{A LHS design for twenty combinations of two inputs}\label{LHS n=20 hori}%
%TCIMACRO{\TeXButton{E}{\end{table}}}%
%BeginExpansion
\end{table}%
%EndExpansion


LHS does not impose a strict mathematical relationship between $n$ and $d$
(whereas a grid with 10 values per input has $10^{d}$ points). Nevertheless,
if LHS uses a \textquotedblleft small\textquotedblright\ $n$ and a
\textquotedblleft large\textquotedblright\ $d$, then LHS covers $[0,1]^{d}$
\textit{sparsely}. For example, if $d$ = 8 (as in the borehole example) and
$n$ = 80, then we sample only 80 gridpoints from the total of $80^{8}$ =
$1.\,7\times10^{15}$ gridpoints \ This sparsity implies that there are only a
few old points close to the new point, so $\widehat{y}$ may have a high
$\sigma(\widehat{y})$. For LHS in Kriging aimed at SA, Loeppky et al. (2009)
gives the\textit{\ rule-of-thumb} $n$ = 10$d$, which implies $n$ $\geq$ 10 if
$d$ $\geq$ 1. (If $n$ $\geq$ 10, then we expect to estimate $\theta_{j}$ quite
accurately.) An example of LHS with $d$ = 2 and $\ n$ = 20 is given
in Table \ref{LHS n=20 hori}; these $n$ = 20 observations cover $[0,1]^{2}$
relatively \textit{sparsely}. In general, the volume of $[l_{d},h_{d}]^{d}$
increases exponentially with $d$, whereas $n$ = 10$d$ increases only linearly
with $d$. This sparsity implies that there are only a few points among the
$n-1$ points in $(\mathbf{X}_{-i},\mathbf{w}_{-i})$) that are
\textit{relatively close} to $\mathbf{x}_{i}$ (to be predicted). The Euclidean
distance between $\mathbf{x}_{i}$ and $\mathbf{x}_{i^{\prime}}$ is
$[\sum_{j=1}^{d}(x_{i;j}-x_{i^{\prime};j})^{2}]^{1/2}$, which indeed increases
as $d$ increases.

There is much \textit{software} for LHS. Microsoft's Excel spreadsheet
software has add-ins that include LHS; also see Oracle's Crystal Ball,
Palisade's @Risk, and Frontline Systems' Risk Solver. LHS is also available in
the MATLAB Statistics toolbox, the R package, the Open TURNS software, and
Sandia's DAKOTA software. Various LHS algorithms are referenced in Kleijnen
(2015, p. 200); recent algorithms are detailed in Dong and Nakayama (2017),
and Le Guiban et al. (2018). (Panagiotopoulos et al. (2018) uses LHS variants
for RBFs---instead of Kriging---metamodels for optimization through genetic
algorithms.)\newline

\textbf{Appendix 4: Boxplots for estimated OK parameters in MC experiment}%
\newline%

%TCIMACRO{\FRAME{ftbpFU}{4.4581in}{3.3313in}{0pt}{\Qcb{Boxplots for
%$\protect\widehat{\mu}_{r}$, $\protect\widehat{\tau}_{r}^{2}$,
%$\protect\widehat{\theta}_{1;r}$ and $\protect\widehat{\theta}_{2;r}$ with $r$
%= 1, ..., 100 in MC experiment with fixed $\QTR{bf}{X}_{20\times2}$ and
%$\QTR{bs}{\psi}$; horizontal dotted lines give this $\QTR{bs}{\psi}$%
%-}}{\Qlb{boxplotGramacy}}{boxplotgramv2.tif}%
%{\special{ language "Scientific Word";  type "GRAPHIC";
%maintain-aspect-ratio TRUE;  display "USEDEF";  valid_file "F";
%width 4.4581in;  height 3.3313in;  depth 0pt;  original-width 10.4374in;
%original-height 7.7712in;  cropleft "0";  croptop "1";  cropright "1";
%cropbottom "0";  filename 'BoxplotGramV2.tif';file-properties "XNPEU";}} }%
%BeginExpansion
\begin{figure}
  \centering
  \includegraphics[scale=0.65]{1BoxplotGramV2.jpg}
  \caption{Boxplots for $\protect\widehat{\mu}_{r}$, $\protect\widehat{\tau}%
_{r}^{2}$, $\protect\widehat{\theta}_{1;r}$ and $\protect\widehat{\theta
}_{2;r}$ with $r$ = 1, ..., 100 in MC experiment with fixed $\mathbf{X}%
_{20\times2}$ and $\boldsymbol{\psi}$; horizontal dotted lines give this
$\boldsymbol{\psi}$}
  \label{boxplotGramacy}
\end{figure}
%EndExpansion
%\bigskip

In general, boxplots use robust estimators for location and scale; namely, the
median and the interquartile (see the textbook on robust statistic, Huber and
Ronchetti (2009)). More specifically, we use\textbf{\ }\textit{boxplots with
whiskers} to summarize the marginal distributions of the components of
$\widehat{\boldsymbol{\psi}}$. After some preliminary experimentation, we let
DACE search for $\widehat{\theta}_{1;r}$ and $\widehat{\theta}_{2;r}$ ($r$ =
1, ..., $m$ = 100), in the range $[0.01,50.00])$. To construct
boxplots,\textbf{\ }we use the MATLAB function \textit{boxplot}. The resulting
Fig. \ref{boxplotGramacy} displays boxplots that include the central mark
denoting the median, the bottom and top edges denoting the 25th and 75th
percentiles, the whiskers denoting the most extreme data points not considered
outliers, and the dots denoting outliers. This Figure is based on
$\mathbf{X}_{20\times2}$ generated by LHS; an example is given in
Table\textbf{ \ref{LHS n=20 hori}.} The $\widehat{\boldsymbol{\psi}}%
$-components should not deviate\ \textquotedblleft too much\textquotedblright%
\ from the true values $\mu$ = $-0.1142$, $\tau^{2}$ = 0.0589, $\theta_{1}$ =
3.8555, and $\theta_{2}$ = 1.1970; these true values are indicated by
horizontal dotted lines in the Figure.\textbf{\ }Visually inspection shows
that the sample median of $\widehat{\mu}_{r}$ lies close to the true mean.
However, the sample median of $\widehat{\tau}_{r}^{2}$ is approximately 0.2,
whereas $\tau^{2}$ = 0.0589. The sample median of $\widehat{\theta}_{1;r}$ is
approximately 2 ($\theta_{1}$ is 3.8555) and the sample median of
$\widehat{\theta}_{2;r}$ is approximately 0.5 ($\theta_{2}$ is 1.1970).
Altogether, this inspection confirms that it is difficult to estimate
$\boldsymbol{\psi}$---as Erickson et al. (2018) also observes (when comparing
software for computing $\widehat{y}$ and $s(\widehat{y})$, which depend on
$\widehat{\boldsymbol{\psi}}$).%

%TCIMACRO{\FRAME{ftbpFU}{5.8894in}{4.1537in}{0pt}{\Qcb{Boxplots for OK
%parameters in MC experiment with 100 replications, 20 combinations of two
%inputs and slope g = 0, 1, 5, and 25 of linear trend; horizontal lines give
%true values of OK parameters}}{\Qlb{boxpsitrend_4g}}{boxpsitrend_4g.tif}%
%{\special{ language "Scientific Word";  type "GRAPHIC";
%maintain-aspect-ratio TRUE;  display "USEDEF";  valid_file "F";
%width 5.8894in;  height 4.1537in;  depth 0pt;  original-width 11.0316in;
%original-height 7.7496in;  cropleft "0";  croptop "1";  cropright "1";
%cropbottom "0";  filename 'BoxPsiTrend_4g.tif';file-properties "XNPEU";}} }%
%BeginExpansion
\begin{figure}
  \centering
  \includegraphics[scale=0.7]{2BoxPsiTrend_4gV3.jpg}
  \caption{ Boxplots for OK parameters in MC experiment with 100 replications, 20
combinations of two inputs and slope g = 0, 1, 5, and 25 of linear trend;
horizontal lines give true values of OK parameters}
  \label{boxpsitrend_4g}
\end{figure}
%EndExpansion

Fig. \ref{boxpsitrend_4g} displays boxplots for the Kriging parameters
$\widehat{\mu}_{r}$, $\widehat{\tau}_{r}^{2}$, $\widehat{\theta}_{1;r}$ and
$\widehat{\theta}_{2;r}$ with $r$ = 1, ..., 100 in the MC experiment with
input $\mathbf{X}_{20\times2}$ and a \textit{linear trend} with slope $\beta$
= $g\tau$ with $g$ = 0, 1, 5, and 25, respectively. We observe that $g$ = 25
gives a box for $\widehat{\theta}_{1;r}$ with zero length; the first, second,
and third quantiles have the same value (namely, 0.4819). Furthermore,
$\widehat{\theta}_{2;r}$ varies between 0.0602 and 0.8839, but 65 (out of 100)
estimates coincide with the sample median.\newline

\textbf{Appendix 5: Details of Gramacy (2016)'s example}\newline%

%TCIMACRO{\FRAME{ftbpFU}{2.5443in}{2.1551in}{0pt}{\Qcb{I/O data $(z_{1;i}%
%,z_{2;i},v_{i})$ with $i$ = 1, ..., 20 of example in Gramacy (2016)}%
%}{\Qlb{I/O Gramacy}}{stemgram.tif}{\special{ language "Scientific Word";
%type "GRAPHIC";  maintain-aspect-ratio TRUE;  display "USEDEF";
%valid_file "F";  width 2.5443in;  height 2.1551in;  depth 0pt;
%original-width 8.4267in;  original-height 7.1252in;  cropleft "0";
%croptop "1";  cropright "1";  cropbottom "0";
%filename 'StemGram.tif';file-properties "XNPEU";}} }%
%BeginExpansion
\begin{figure}
  \centering
  \includegraphics[scale=0.6]{3StemGram.jpg}
  \caption{I/O data $(z_{1;i},z_{2;i},v_{i})$ with $i$ = 1, ..., 20 of example
in Gramacy (2016)}
  \label{I/O Gramacy}
\end{figure}
%EndExpansion
%

%TCIMACRO{\FRAME{ftbpFU}{4.9095in}{2.156in}{0pt}{\Qcb{Scatterplots
%$(z_{1;i},v_{i})$ and $(z_{2;i},v_{i})$ using the I/O data of Fig. 3}%
%}{\Qlb{Scatterplots Gramacy}}{graphgramm_z1v_z2v.tif}%
%{\special{ language "Scientific Word";  type "GRAPHIC";
%maintain-aspect-ratio TRUE;  display "USEDEF";  valid_file "F";
%width 4.9095in;  height 2.156in;  depth 0pt;  original-width 13.4997in;
%original-height 5.8851in;  cropleft "0";  croptop "1";  cropright "1";
%cropbottom "0";  filename 'GraphGramm_z1v_z2v.tif';file-properties "XNPEU";}}
%}%
%BeginExpansion
\begin{figure}
  \centering
  \includegraphics[scale=0.65]{4GraphGramm_z1v_z2v.jpg}
  \caption{Scatterplots $(z_{1;i},v_{i})$
and $(z_{2;i},v_{i})$ using the I/O data of Fig. 3}
 \label{Scatterplots Gramacy}
\end{figure}
%EndExpansion


Fig. \ref{I/O Gramacy} is a 3D plot with $n$ = 20 vertical lines that start at
$(z_{1;i},z_{2;i},0)$ with $i$\ = 1, ..., $n$\ = 20 (empty circles) and end at
$(z_{1;i},z_{2;i},v_{i})$ (solid circles; each $v_{i}$ is
negative). Fig. \ref{Scatterplots Gramacy} displays the two corresponding
scatterplots $(z_{j;i},v_{i})$ with $j$ = 1, 2.%

%TCIMACRO{\FRAME{ftbpFU}{2.5443in}{2.1551in}{0pt}{\Qcb{Gramacy (2016)' s
%example with I/O data at $N$ = 40,401 points}}{\Qlb{Gramacy fig}%
%}{gramfunction.tif}{\special{ language "Scientific Word";  type "GRAPHIC";
%maintain-aspect-ratio TRUE;  display "USEDEF";  valid_file "F";
%width 2.5443in;  height 2.1551in;  depth 0pt;  original-width 8.4267in;
%original-height 7.1252in;  cropleft "0";  croptop "1";  cropright "1";
%cropbottom "0";  filename 'GramFunction.tif';file-properties "XNPEU";}} }%
%BeginExpansion
\begin{figure}
  \centering
  \includegraphics[scale=0.2]{5GramFunction.jpg}
  \caption{Gramacy (2016)' s example with I/O data at $N$ = 40,401 points}
\label{Gramacy fig}
\end{figure}
%EndExpansion


Actually, $(z_{1;i},z_{2;i},v_{i})$ is determined by the following
\textit{white box} (explicit mathematical) function:%
\begin{align}
v(z_{1},z_{2})  &  =-f(z_{1})f(z_{2})\text{ with }-2\leq z_{j}\leq2\text{
and}\nonumber\\
f(z_{j})  &  =e^{-(z_{j}-1)^{2}}+e^{-0.8(z_{j}+1)^{2}}-0.05\sin[8(z_{j}%
+0.1)]\text{.} \label{Gramacy example}%
\end{align}
Gramacy (2016) gives\textbf{ }the 3D plot of this function using as many as
$N$ = 40,401 combinations defined by the $201\times201$ grid in $[-2,2]^{2}%
$; see Fig\textbf{.} \ref{Gramacy fig}. This plot shows that the
function is very \textquotedblleft wiggly\textquotedblright, which agrees with
our interpretation of the plots with only $n$ = 20 I/O data. Both
Fig. \ref{Scatterplots Gramacy} and Fig. \ref{Gramacy fig} show that $v$ is
extreme if $z_{j}$ is extreme, whereas OK assumes a constant mean---which
seems a reasonable assumption (the borehole example,\textbf{\ }however, has a
monotonically increasing I/O function; see Appendix 6).

Note: Gramacy (2016) considers $N$ = 40,401 I/O combinations, resulting in
numerical problems when computing $\widehat{y}$ and $s(\widehat{y})$.
Therefore, that publication selects only $n$ = 50 combinations---among these
$N$ combinations---that lie close to $x_{0}$, using a rather complicated
selection procedure. Kleijnen and Van Beers (2020), however, simply selects
the $n$ nearest neighbors of $\mathbf{x}_{0}$.%

%TCIMACRO{\FRAME{ftbpFU}{4.3111in}{3.0441in}{0pt}{\Qcb{Blue circles denote
%$(\protect\widehat{\theta}_{1;r},\protect\widehat{\theta}_{2;r})$ with $r$ =
%1, ..., 100 and 0 $\leq$ $\protect\widehat{\theta}_{j;r}$ $\leq$ 50 ($j$ = 1,
%2) in 100 LHS designs for Gramacy (2016)'s example; red line has slope $1$ and
%intercept $0$; black line denotes first-order polynomial fitted through least
%squares}}{\Qlb{scatter thetas}}{scatterthetas.tif}%
%{\special{ language "Scientific Word";  type "GRAPHIC";
%maintain-aspect-ratio TRUE;  display "USEDEF";  valid_file "F";
%width 4.3111in;  height 3.0441in;  depth 0pt;  original-width 8.0834in;
%original-height 5.6879in;  cropleft "0";  croptop "1";  cropright "1";
%cropbottom "0";  filename 'ScatterThetas.tif';file-properties "XNPEU";}} }%
%BeginExpansion
\begin{figure}
  \centering
  \includegraphics[scale=0.6]{6ScatterThetas.jpg}
  \caption{Blue circles denote $(\protect\widehat{\theta}_{1;r}%
,\protect\widehat{\theta}_{2;r})$ with $r$ = 1, ..., 100 and 0 $\leq$
$\protect\widehat{\theta}_{j;r}$ $\leq$ 50 ($j$ = 1, 2) in 100 LHS designs for
Gramacy (2016)'s example; red line has slope $1$ and intercept $0$; black line
denotes first-order polynomial fitted through least squares}
\label{scatter thetas}
\end{figure}
%EndExpansion


Using\textbf{ }$(\mathbf{Z},v)$, DACE gives $\widehat{\theta}_{1}$ = $3.8555$
and $\widehat{\theta}_{2}$ = $1.1970$. We find it hard to explain\textbf{\ }%
these specific values (also see our discussion of the prior distribution for
$\theta_{j}$ in the Bayesian approach to GP, in Section 1). We may compare the
\textit{relative }values of $\widehat{\theta}_{j}$ in this example, because we
(not OK) know the white-box simulation model (\ref{Gramacy example}) so
$\theta_{1}/\theta_{2}$ = 1. However, DACE gives $\widehat{\theta}%
_{1}/\widehat{\theta}_{2}$ = $3.2$. To further investigate $\widehat{\theta
}_{1}/\widehat{\theta}_{2}$, we \ sample $\mathbf{(Z}_{n;r},v_{n;r})$ with $n$
= 20 and $r$ = 1, ..., 100. DACE's search for $\widehat{\theta}_{1}$ and
$\widehat{\theta}_{2}$ uses the same range $[0;50]$ in all 100
macroreplications. These macroreplications give Fig. \ref{scatter thetas},
where the blue dots denote $(\widehat{\theta}_{1;r},\widehat{\theta}_{2;r})$.
The pair\textbf{ }$(0.63,10)$ occurs four times, $(0.078,25)$ occurs two
times, and\textbf{ }$(0.04,20)$ and $(0.16,20$) are very close. Two
macroreplications give $\widehat{\theta}_{2}$ = 50, which is the upper limit
in DACE's search for the MLEs (in practice, we would relax this limit, and
continue DACE's search). Furthermore, 53 of these 100 pairs give
$\widehat{\theta}_{2}$
%TCIMACRO{\TEXTsymbol{>} }%
%BeginExpansion
$>$
%EndExpansion
$\widehat{\theta}_{1}$; i.e., 53 pairs lie above the line with slope $1$ and
intercept $0$ (see the solid, non-dashed red line). This result does not
contradict our preceding\textbf{ }statement $\theta_{1}/\theta_{2}$ = 1. The
first-order polynomial fitted through least squares (LS), has intercept
3.89\ and slope 0.94; see the dashed black line. This intercept is close to
the \textquotedblleft true\textquotedblright\ value $\theta_{1}$ = $\theta
_{2}$ $\thickapprox$ 3 where the value 3 is explained in the next Note. Most
of the 100 (estimated) pairs cluster around the (true) combination $(3,3)$;
i.e., DACE's search does give \textquotedblleft reasonable\textquotedblright%
\ estimates. We let $\rho_{\boldsymbol{\theta}}$ denote Pearson's correlation
coefficient for $(\widehat{\theta}_{1},\widehat{\theta}_{2})$; $\rho$ is
further discussed in Kleijnen (2015, p. 113). Then the fitted line gives
$\widehat{\rho}_{\boldsymbol{\theta}}$ = 0.49. The LS criterion gives results
that are sensitive to outliers, so we might remove outliers such as the two
combinations with $\widehat{\theta}_{2}$ = 50. However, we do not further
refine our analysis of the fitted line, because this line does not contradict
our statement $\theta_{1}$ = $\theta_{2}$. The rather low value of
$\widehat{\rho}_{\boldsymbol{\theta}}$ (namely, 0.49) is explained by the
rather high standard deviations of $\widehat{\theta}_{1}$ and $\widehat{\theta
}_{2}$; namely, 6.98 and 13.47, while the sample means are 5.08 and 8.68 and
the sample medians are 2.50 and 3.54.

Note: We do obtain $\widehat{\theta}_{1}$ $\approx\widehat{\theta}_{2}$ if we
use a grid instead of LHS; e.g., grids with 4, 5, 10, and 25
observations\textbf{ }per input $z_{j}$ give the following pairs
$(\widehat{\theta}_{1},\widehat{\theta}_{2})$: $(0.3125,0.2210)$,
$(2.5000,2.4148)$, $(2.9730,2.9730)$, and $(3.2421,3.3856)$. Because\textbf{
}$\boldsymbol{\psi}$ $\boldsymbol{\neq}$ $\widehat{\boldsymbol{\psi}}$, the
predictor $\widehat{y}$ is not the BLUP. To further investigate Gramacy's
example, we replace the OK metamodel by a simpler metamodel; namely,
a \textit{second-degree polynomial}. We expect that in this
polynomial the first-order effects of $z_{1}$ and $z_{2}$ are the same, and so
are their purely quadratic effects. Indeed, when we fit this polynomial to the
I/O data of our LHS design with $n$ = 20, then we obtain \ estimated
first-order effects\ that are nearly the same (namely, 0.03), and purely
quadratic effects that are also nearly the same (namely, 0.09) (the two-factor
interaction is 0.00, and the intercept is -0.98). The $4\times4$ grid gives
estimated first-order effects\ that are exactly the same (namely, 0.005), and
purely quadratic effects that are also exactly the same (namely, 0.104) (the
two-factor interaction is -0.000, and the intercept is -0.963). Altogether,
the LHS and the grid designs confirm our expectation; namely, $z_{1}$ and
$z_{2}$ have the same first-order and purely quadratic effects in a
second-degree polynomial metamodel.%
%TCIMACRO{\FRAME{ftbpFU}{3.1557in}{2.3497in}{0pt}{\Qcb{Scatterplot for
%$(v_{i},\protect\widehat{y}_{-i})$ augmented with $\protect\widehat{y}_{-i}\pm
%z_{\alpha_{\text{E}}/(2n)}s(\protect\widehat{y}_{-i})$ with $n$ = 20, and
%$\protect\widehat{y}_{-i}$ = $v_{i}$, in Gramcy (2016)'s example}%
%}{\Qlb{Scatter Gramacy}}{fig7a_markers_star.tif}%
%{\special{ language "Scientific Word";  type "GRAPHIC";
%maintain-aspect-ratio TRUE;  display "USEDEF";  valid_file "F";
%width 3.1557in;  height 2.3497in;  depth 0pt;  original-width 10.4374in;
%original-height 7.7496in;  cropleft "0";  croptop "1";  cropright "1";
%cropbottom "0";  filename 'Fig7A_Markers_Star.tif';file-properties "XNPEU";}}
%}%
%BeginExpansion
\begin{figure}
  \centering
  \includegraphics[scale=0.6]{7Fig7A_Markers_Star.jpg}
  \caption{Scatterplot for $(v_{i},\protect\widehat{y}_{-i})$ augmented with
$\protect\widehat{y}_{-i}\pm z_{\alpha_{\text{E}}/(2n)}s(\protect\widehat{y}%
_{-i})$ with $n$ = 20, and $\protect\widehat{y}_{-i}$ = $v_{i}$, in Gramcy
(2016)'s example}
\label{Scatter Gramacy}
\end{figure}
%EndExpansion


Fig. \ref{Scatter Gramacy} displays the augmented scatterplot; all $n$ = 20
CIs cover $\widehat{y}_{-i}$ = $v_{i}$, so our test does not reject the OK
metamodel. The corresponding $\max|$PES$_{i}|$ is 2.33, which is not
significant (because $z_{1-[\alpha_{\text{E}}/(2n)]}$ = $z_{1-\alpha
_{\text{E}}/40}$ so $\alpha_{\text{E}}$ = 0.20 gives $z_{0.995}$ = 2.58;
obviously, $\alpha_{\text{E}}$
%TCIMACRO{\TEXTsymbol{<} }%
%BeginExpansion
$<$
%EndExpansion
0.20 gives $z_{1-\alpha_{\text{E}}/40}$
%TCIMACRO{\TEXTsymbol{>} }%
%BeginExpansion
$>$
%EndExpansion
2.58). \newline

\textbf{Appendix 6: Details of the borehole example}\newline%

%TCIMACRO{\FRAME{ftbpFU}{3.864in}{3.2197in}{0pt}{\Qcb{Scatterplots
%$(z_{j;i},v_{i})$ with $j$ = 1, ..., 8 and $i$ = 1, ...., 80 using I/O data of
%borehole model}}{\Qlb{2 scatters borehole}}{borehole8scatters.tif}%
%{\special{ language "Scientific Word";  type "GRAPHIC";  display "USEDEF";
%valid_file "F";  width 3.864in;  height 3.2197in;  depth 0pt;
%original-width 8.0834in;  original-height 9.9791in;  cropleft "0";
%croptop "1";  cropright "1";  cropbottom "0";
%filename 'Borehole8Scatters.tif';file-properties "XNPEU";}} }%

%BeginExpansion
\begin{figure}
  \centering
  \includegraphics[scale=0.65]{8Borehole8Scatters.jpg}
  \caption{Scatterplots $(z_{j;i},v_{i})$ with $j$ = 1, ..., 8 and $i$ = 1,
...., 80 using I/O data of borehole model}
\label{2 scatters borehole}
\end{figure}
%EndExpansion


Fig. \ref{2 scatters borehole} displays $d$ = 8 scatterplots $(z_{j;i},v_{i})$
with $j$ = 1, ..., 8 and $i$ = 1, ..., 80. These plots suggest that $z_{1}$ =
$r_{w}$ has the strongest (positive) effect, whereas the other seven inputs
have no clear (main or first-order) effects. So a better metamodel might use
UK with a first-order polynomial in $r_{w}$; however, we limit our
investigation to OK. This borehole simulation model uses the
following white-box function:%
\begin{equation}
v=\frac{2\pi T_{u}(H_{u}-H_{l})}{\ln(r/r_{w})(1+\frac{2LT_{u}}{\ln
(r/r_{w})r_{w}^{2}d_{w}}+\frac{T_{u}}{T_{l}})}.\label{borehole eq}%
\end{equation}
For this model,\textbf{ }there are computer codes in MATLAB and R available on

https://www.sfu.ca/\symbol{126}ssurjano/borehole.html.%

%TCIMACRO{\FRAME{ftbpFU}{4.1226in}{1.4451in}{0pt}{\Qcb{Borehole example: I/O
%function with $N$ = $21\times21$ input combinations $(r_{w},L)$ (left pane)
%and $(r,T_{u})$ (right pane) versus output $w$}}{\Qlb{borehole plots}%
%}{graphboreh_rw-l_r-tu.tif}{\special{ language "Scientific Word";
%type "GRAPHIC";  maintain-aspect-ratio TRUE;  display "USEDEF";
%valid_file "F";  width 4.1226in;  height 1.4451in;  depth 0pt;
%original-width 12.5519in;  original-height 4.3647in;  cropleft "0";
%croptop "1";  cropright "1";  cropbottom "0";
%filename 'GraphBoreH_rw-l_r-Tu.tif';file-properties "XNPEU";}} }%

%BeginExpansion
\begin{figure}
  \centering
  \includegraphics[scale=0.5]{9GraphBoreH_rw-l_r-Tu.jpg}
  \caption{Borehole example: I/O function with $N$ = $21\times21$ input
combinations $(r_{w},L)$ (left pane) and $(r,T_{u})$ (right pane) versus
output $w$}
\label{borehole plots}
\end{figure}
%EndExpansion



The white box defined by (\ref{borehole eq}) is nonlinear and has non-additive
effects. Our \textit{interpretation} of the simulation I/O data uses the
original data, $(\mathbf{Z,v)}$. We make 3D plots (as we did in Fig.
\ref{Gramacy fig} for Gramacy's example); i.e., we plot the output $v_{i}$
($i$ = 1, ...., $n$) versus two original inputs while keeping the other six
original inputs constant at the midpoints of their ranges. Santner et al.
(2018) estimates that $r_{w}$ is the most important input and that the three
inputs $L$, $H_{l}$, and $H_{u}$ have approximately equally important effects.
In the left-hand pane of Fig. \ref{borehole plots} we use a grid of $N$ =
$21\times21$ input combinations $(r_{w},L)$ while keeping the other six inputs
constant; the right-hand pane gives a similar plot for $(r,T_{u})$. We point
out that the range of $w$ is relatively small in the right-hand\textbf{\ }%
panel, so the effect of $r$ is relatively small (the effect of $T_{u}$ is
unimportant). We also make such 3D plots for all other combinations of two
inputs; these plots are very similar, so we do not display them. Moreover, we
make contourplots. Together these\textbf{\ }3D plots and contourplots confirm
Santner et al.'s conclusion that $r_{w}$ has the most important main effect,
and has important interactions. Moreover, we conclude that the I/O function is
\textit{monotonic}; i.e., the borehole model looks complicated, but seems to
give a simpler I/O function than Gramacy's example in Appendix 5.

DACE with $(\mathbf{X}_{80\times8},\mathbf{v}_{80})$ gives
$\widehat{\boldsymbol{\psi}}$. This $\widehat{\boldsymbol{\psi}}$ includes
$\widehat{\mu}$ = 76.74, which agrees with the average output in the left-hand
pane of Fig. \ref{borehole plots}. Furthermore, $\widehat{\tau}$ = 40.10,
which also agrees with the spread in the outputs in this pane. Finally, DACE
searches for $\widehat{\theta}_{j}$ in the interval $[0.001,5]$, which is
based on some preliminary tests with DACE. This search gives the following
$\widehat{\theta}_{j}$-values, where we display the first three decimal units
because the lower limit of this search (namely, 0.001) also uses three
decimals: 0.310, 0.067, 0.015, 0.017, 0.001, 0.007, 0.011, 0.002. So
$\widehat{\theta}_{5}$ = 0.001 equals the lower limit; i.e., a better estimate
of $\theta_{5}$ might have a lower value than 0.001.

Originally we conjectured that $\widehat{\theta}_{1}$ would be the smallest of
the $\widehat{\theta}_{j}$-values, because $z_{1}$ = $r_{w}$ seems the most
important input in Fig. \ref{borehole plots}. The magnitudes of
$1/\widehat{\theta}_{j}$ (not $\widehat{\theta}_{j}$) are also shown in the
boxplot in Sun et al. (2018, p. 12), which confirms our results. More
precisely, our definition of the Gaussian correlation function equals the
definition in Lophaven et al. (2002, p. 6), which equals the definition in Sun
et al. (2018, p. 6)---provided we replace our $\theta_{j}h_{j}^{2}$ by
$h_{j}^{2}/\theta_{j}$. Moreover, different software may give different
$\widehat{\theta}$, as Erickson et al. (2018) demonstrates.

Now, however, we revisit our original conjecture; i.e., now we notice that the
Gaussian correlation function is determined not only by $\widehat{\theta}_{j}%
$, but also by $h_{j}^{2}$ = $(x_{g;j}-x_{g^{\prime};j})^{2}$ with
$g,g^{\prime}$ = 1, ..., $n$. The original inputs $z_{i;j}$ give (say)
$\theta_{z;j}$, and $h_{z;j}$, which depend on the range $[l_{j},u_{j}]$. In
$\mathbf{Z}_{n\times d}$ (determined by $\mathbf{X}_{n\times d}$ sampled by
LHS) the common length of the $n$ subranges of $z_{j}$ is $(u_{j}-l_{j})/n$.
So the observed distances between two outputs at the midpoints of these
subranges are a multiple of this length; e.g., there are $n-1$ observations on
outputs with inputs at the distance $(u_{j}-l_{j})/n$ and there is a single
observation on the two outputs with the longest distance $(u_{j}%
-l_{j})(n-1)/n$. In general, there are $i$ observations on outputs with inputs
that are $(u_{j}-l_{j})(n-i)/n$ apart (with $i$ = 1, ..., $n$). The range of
$z_{1}$ = $r_{w}$ is very small compared with the ranges of the other seven
original inputs (Gramacy (2016)'s example has inputs with the same range;
namely, $[-2,2]$).%

%TCIMACRO{\FRAME{ftbpFU}{3.5708in}{3.2041in}{0pt}{\Qcb{Scatterplot for
%$(v_{i},\protect\widehat{y}_{-i})$ augmented with $\protect\widehat{y}_{-i}\pm
%z_{\alpha_{\text{E}}/(2n)}s(\protect\widehat{y}_{-i})$ with $\alpha_{\text{E}%
%}$ = 0.20 and $n$ = 80, and $\protect\widehat{y}_{-i}$ = $v_{i}$, in borehole
%example}}{\Qlb{Scatter borehole}}{augmentedscatterboreholerevised.tif}%
%{\special{ language "Scientific Word";  type "GRAPHIC";
%maintain-aspect-ratio TRUE;  display "USEDEF";  valid_file "F";
%width 3.5708in;  height 3.2041in;  depth 0pt;  original-width 8.1042in;
%original-height 7.2601in;  cropleft "0";  croptop "1";  cropright "1";
%cropbottom "0";
%filename 'AugmentedScatterBoreholeRevised.tif';file-properties "XNPEU";}} }%

%BeginExpansion
\begin{figure}
  \centering
  \includegraphics[scale=0.65]{10AugmentedScatterBoreholeRevised.jpg}
  \caption{Scatterplot for $(v_{i},\protect\widehat{y}_{-i})$ augmented with
$\protect\widehat{y}_{-i}\pm z_{\alpha_{\text{E}}/(2n)}s(\protect\widehat{y}%
_{-i})$ with $\alpha_{\text{E}}$ = 0.20 and $n$ = 80, and $\protect\widehat{y}%
_{-i}$ = $v_{i}$, in borehole example}
\label{Scatter borehole}
\end{figure}
%EndExpansion
\newline

We summarize the results of our LOO-CV test, in the
augmented scatterplot with 80 CIs and \ $\alpha_{\text{E}}$ = 0.20,
\textbf{in }Fig. \ref{Scatter borehole}.\newline

\textbf{References}\newline

Afendras, G. and M. Markatou (2019), Optimality of training/test size and
resampling effectiveness in cross-validation. \textit{Journal of Statistical
Planning and Inference}, 199, pp. 86--301

Antunes, F., M. Amorim, B. Ribeiro, and F. Pereira (2019), Active learning
metamodeling for policy analysis: application to an emergency medical service
simulator. \textit{Simulation Practice and Theory}, \ accepted

Bacchi, V., E. Antoshchenkova, H. Jomard, L. Bardet, C-M. Duluc, O. Scotti,
and H. Hebert (2018), Development of a methodological framework for the
assessment of seismic induced tsunami hazard through uncertainty
quantification: application to the Azores-Gibraltar Fracture Zone.
\textit{Natural Hazards and Earth System Sciences} (under review)

Bartz-Beielstein, T. (2016), Stacked generalization of surrogate models; a
practical approach. TH K\"{o}ln

Bastos,L.S. and A. O'Hagan (2009), Diagnostics for Gaussian process emulators.
\textit{Technometrics}, 51, no. 4, pp. 425--438

Bhosekar, A. and M. Ierapetritou (2018), Advances in surrogate based modeling,
feasibility analysis, and optimization: a review. \textit{Computers and
Chemical Engineering}, 108, pp. 250--267

Da Costa, J.J., F. Chainet, B. Celse,M. Lacoue-Negre, C. Ruckebusch, N.
Caillol, and D. Espinat (2018), Kriging modeling to predict viscosity index of
base oils.\textit{\ Energy Fuels}, 32, pp. 2588--2597

De Carvalho, R.N., G.B. Machado, and M.J. Cola\c{c}o (2017), Estimating
gasoline performance in internal combustion engines with simulation
metamodels. \textit{Fuel}, 93, pp. 230--240

Demory, B., M. Henner, A. Lagnoux, and T.M.N. Nguyen (2019). Expected
improvement applied to an industrial context: prediction of new geometries
increasing the efficiency of fans. ffhal-02044258f

Dong, H. and M.K. Nakayama (2017), Quantile estimation with Latin hypercube
sampling. \textit{Operations Research}, 65, no. 6, pp. 1678--1695

Dubrule, O. (1983), Cross validation of kriging in a unique neighborhood.
\textit{Mathematical Geology} 15, pp. 687--699

Efron, B. (2015), Frequentist accuracy of Bayesian estimates. Royal
Statistical Society, Series B, 77, no. 3, pp. 617--646

Erickson, C.B., B.E. Ankenman, S.M. Sanchez (2018), Comparison of Gaussian
process modeling software.\textit{\ European Journal of Operational Research},
266, pp. 179--192

Fei, L., J.Y. Dantan, C. Baudouin, and S. Du (2019), Calibration and
uncertainty estimation of non-contact coordinate measurement systems based on
Kriging models. \textit{Precision Engineering}, 57, pp. 16--29

Forrester, A. and A. Keane (2009), Recent advances in surrogate-based
optimization. \textit{Progress in Aerospace Sciences}, 45, no. 1--3, pp. 50--79

Fuhg, J.N., A. Fau, and U. Nackenhorst (2020), State-of-the-art and
comparative review of adaptive sampling methods for Kriging. \textit{Archives
of Computational Methods in Engineering}, in press

Garbo, A. and B.J. German (2019), Performance assessment of a cross-validation
sampling strategy with active surrogate model selection. \textit{Structural
and Multidisciplinary Optimization}, 59, no. 6, pp. 2257--2272

Ginsbourger, D. and C. Sch\"{a}rer (2021), Fast calculation of Gaussian
process multiple-fold cross-validation residuals and their covariances.
arXiv:2101.03108v1 [stat.ME] 8 Jan 2021

Gorissen, D., I. Couckuyt, E. Laermans, and T. Dhaene (2010), Multiobjective
global surrogate modeling, dealing with the 5-percent problem.
\textit{Engineering with Computer}s, 26, pp. 81--98

Gramacy, R.B. (2016), laGP: large-scale spatial modeling via local approximate
Gaussian processes in R. \textit{Journal of Statistical Software}, 72, no. 1,
pp. 1--46

Gramacy, R.B. (2020), \textit{Surrogates: Gaussian process modeling, design,
and optimization for the applied sciences}. CRC Press, Taylor \& Francis
Group, Boca Raton

Huber, P.J. and E.M. Ronchetti (2009), \textit{Robust statistics}. Wiley, New York

Ji, D., W. Dong, T. Hong, T. Dai, Z. Zheng, S. Yang, and X. Zhu (2018),
Assessing parameter importance of the weather research and forecasting model
based on global sensitivity analysis methods.\textit{\ Journal of Geophysical
Research: Atmospheres}, 123, in press

Jin, S-S and H-J Jung (2016), Self-adaptive sampling for sequential surrogate
modeling of time-consuming finite element analysis. \textit{Smart Structures
and Systems}, 17, no. 4, pp. 611--629

Kamath, C. and Y.J. Fan (2018, Regression with small data sets: a case study
using code surrogates in additive manufacturing. \textit{Knowledge and
Information Systems}, 57, pp. 475--493

Kleijnen, J.P.C. (2015), \textit{Design and analysis of simulation
experiments; second edition}. Springer

Kleijnen, J. P. C. (2020), Kriging: methods and applications. \textit{Handbook
on Model Order Reduction. Volume I: Methods \& Algorithms}, edited by P.
Benner, G. Rozza, S. Grivet-Talocia, W.H.A. Schilders, A. Quarteroni, and L.M.
Silveira, Walter De Gruyter, Berlin (a preprint is available on

https://sites.google.com/site/kleijnenjackpc/home/news)

Kleijnen, J.P.C. and W. van Beers (2020), Prediction for big data through
Kriging: small sequential and one-shot designs. \textit{American Journal of
Mathematical and Management Sciences}, Open Access

https://www.tandfonline.com/doi/full/10.1080/01966324.2020.1716281

Kyprioti, A.P., A.A. Tafanidis1, N.C. Nadal Caraballo, and M.O. Campbel
(2020), Incorporation of sea level rise in storm surge surrogate modeling.
\textit{Natural Hazards}, https://doi.org/10.1007/s11069-020-04322-z

Law, A.M. (2015),\textit{\ Simulation modeling and analysis; fifth edition}.
McGraw-Hill, Boston

Le Guiban, K., A. Rimmel, M-A. Weisser, and J. Tomasik (2018), The first
approximation algorithm for the maximin Latin hypercube design problem.
\textit{Operations Research}, 66, no. 1, pp. 253--266

Lin, Y., B.L. Nelson, and L. Pei (2019), Virtual statistics in simulation via
$k$ nearest neighbors. \textit{INFORMS Journal on Computing}, 31, no. 3, pp. 576--592

Luminari, N., C. Airiau, and A. Bottaro (2018), Effects of porosity and
inertia on the apparent permeability tensor in fibrous media.
\textit{International Journal of Multiphase Flow}, 106, pp. 60--74

Loeppky, J.L., J. Sacks, and W. Welch (2009) Choosing the sample size of a
computer experiment: a practical guide. \textit{Technometrics}, 51, no. 4, pp. 366--376

Lophaven, S.N., H.B. Nielsen, and J. Sondergaard (2002), \textit{DACE: a
Matlab Kriging toolbox, version 2.0}. IMM Technical University of Denmark,
Kongens Lyngby

Lupera Calahorrano, G.J., A. Shokry, G. Campanya, and A. Espu\~{n}a (2016),
Application of the meta-multiparametric methodology to the control of
emissions in the industry under continuous and discrete uncertain parameters.
\textit{Chemical Engineering Research and Design}, 115, pp. 365--373

Mubarak, N., I. Hussain, M. Faisal, T. Hussain, M.Y. Shad, N.M. AbdEl-Salam,
and J. Shabbir (2015), Spatial distribution of sulfate concentration in
groundwater of South-Punjab, Pakistan.\textit{ Water Quality Exposure and
Health,} DOI: 10.1007/s12403-015-0165-7

Nahvi, A., M.K. Sadoughi, A. Arabzadeh, A. Sassani, C. Hu, H. Ceylan, and S.
Kim (2019), Multi-objective Bayesian optimization of super hydrophobic
coatings on asphalt concrete surfaces. \textit{Journal of Computational Design
and Engineering}, 6, pp. 693--704

Panagiotopoulos, D., O. Iqbal, Z. Mourelatos, and D. Papadimitriou (2018),
Optimal water jacket flow distribution using a new group-based space-filling
design of experiments algorithm. SAE Technical Paper 2018-01-1017, 2018, doi:10.4271/2018-01-1017

Parnianifard, A., A.S. Azfanizam, M.K.A. Ariffin, and M.I.S. Ismail (2018),
Kriging-assisted robust black-box simulation optimization in direct speed
control of DC motor under uncertainty.\textit{\ IEEE Transactions on
Magnetics}, pp. 1--10

Quirante, N., J. Javaloyes-Ant\'{o}n, and J.A. Caballero (2018), Hybrid
simulation-equation based synthesis of chemical processes.\textit{\ Chemical
Engineering Research and Design}, 132, pp. 766--784

Rasmussen, C.E. and C.K.I. Williams (2006), \textit{Gaussian processes for
machine learning}. The MIT Press,

http://www.gaussianprocess.org/gpml/chapters/RW.pdf

Roustant, O., D. Ginsbourger, and Y. Deville (2012). DiceKriging, DiceOptim:
two R packages for the analysis of computer experiments by Kriging-based
metamodeling and optimization. \textit{Journal of Statistical Software}, 51,
no. 1, pp. 1---55

Santner, T.J., B.J. Williams, and W.I. Notz (2018), \textit{The design and
analysis of computer experiments; second revised edition}. Springer, New York

Sch\"{o}bi,R., B. Sudret, and J. Wiart (2015), Polynomial-chaos-based
Kriging.\textit{ International Journal for Uncertainty Quantification}, 5, no.
2, pp. 171--193

Shu, L., P. Jiang, L. Wan, Q. Zhou, X. Shao, and Y. Zhang, (2017),
Metamodel-based design optimization employing a novel sequential sampling
strategy\textit{.Engineering Computations}, 34, no. 8, pp.2547--2564

Stone, M. (1974), Cross-validatory choice and assessment of statistical
predictions. \textit{Journal Royal Statistical} \textit{Society, Series B},
36, no. 2, pp. 111--147

Strano, M., Q. Semeraro, L. Iorio, and R. Sofia (2018), Hierarchical
metamodeling of the air bending process. \textit{Journal of Manufacturing
Science and Engineering}, 140, pp. 1--19

Sun, F., R.B. Gramacy, B. Haaland, E. Lawrence , and A. Walker (2018),
Evaluating satellite drag from large simulation experiments.
arXiv:1712.00182v2 [stat.CO] 21 Aug 2018

Van Steenkiste, T., J. van der Herten, I. Couckuyt, and T. Dhaene (2018),
Sequential sensitivity analysis of expensive black-box simulators with
metamodelling. \textit{Applied Mathematical Modelling, }61, pp. 66 8--6 81

Vi\'{u}dez-Moreiras, D. (2019), Performance influences on metamodelling for
aerodynamic surrogate-based optimization of an aerofoil.\textit{\ Engineering
Optimization}, 51:3, pp. 427--446

Wang, W. and B. Haaland (2019), Controlling sources of inaccuracy in
stochastic kriging. \textit{Technometrics}, 61, no. 3, pp. 309--321

Williams, C.A., K.C.A. Wedgwood, H. Mohammadi, K. Prouse, O.W. Tomlinson, and
K. Tsaneva-Atanasova (2019), Cardiopulmonary responses to maximal aerobic
exercise in patients with cystic fibrosis. \textit{PloS One}, 14, no. 2, https://doi.org/10.1371/journal.pone.0211219

Xiao, N-C. M.J. Zuo,and W. Guo (2018), Efficient reliability analysis based on
adaptive sequential sampling design and cross-validation. \textit{Applied
Mathematical Modelling}, 58, pp. 404--420

Yin, S., Z. Wang, Z. Zhu, X. Zou, and W. Wang (2018), Using Kriging with a
heterogeneous measurement error to improve the accuracy of extreme
precipitation return level estimation. \textit{Journal of Hydrology}, 562, pp. 518--529

Zhang, J., A.A. Taflanidis, and J. C. Medina (2017), Sequential approximate
optimization for design under uncertainty problems utilizing Kriging
metamodeling in augmented input space. \textit{Computer Methods in Applied
Mechanics and Engineering}, 315, pp. 369--395

Zhang, J., A.A. Taflanidis, N.C. Nadal-Caraballo, J.A. Melby, and F. Diop
(2018), Advances in surrogate modeling for storm surge prediction: storm
selection and addressing characteristics related to climate change.
\textit{Natural Hazards,} 94, pp. 1225--1253


\end{document}