# Statistical Uncertainty Analysis for Stochastic Simulation

Wei Xie1, Barry L. Nelson2, Russell R. Barton3 Corresponding author: Wei Xie
(e-mail: w.xie@northeastern.edu). 1Department of Mechanical and Industrial
Engineering, Northeastern University 2Department of Industrial Engineering and
Management Sciences, Northwestern University 3The Mary Jean and Frank P. Smeal
College of Business Administration, Pennsylvania State University

###### Abstract

When we use simulation to evaluate the performance of a stochastic system, the
simulation often contains input distributions estimated from real-world data;
therefore, there is both simulation and input uncertainty in the performance
estimates. Ignoring either source of uncertainty underestimates the overall
statistical error. Simulation uncertainty can be reduced by additional
computation (e.g., more replications). Input uncertainty can be reduced by
collecting more real-world data, when feasible. This paper proposes an
approach to quantify overall statistical uncertainty when the simulation is
driven by independent parametric input distributions; specifically, we produce
a confidence interval that accounts for both simulation and input uncertainty
by using a metamodel-assisted bootstrapping approach. The input uncertainty is
measured via bootstrapping, an equation-based stochastic kriging metamodel
propagates the input uncertainty to the output mean, and both simulation and
metamodel uncertainty are derived using properties of the metamodel. A
variance decomposition is proposed to estimate the relative contribution of
input to overall uncertainty; this information indicates whether the overall
uncertainty can be significantly reduced through additional simulation alone.
Asymptotic analysis provides theoretical support for our approach, while an
empirical study demonstrates that it has good finite-sample performance.

###### keywords:

Input uncertainty, confidence intervals, bootstrap, stochastic kriging,
simulation output analysis, metamodel

##  1 Introduction

Stochastic simulation is used to characterize the behavior of complex systems
that are driven by random input processes. By “input process” we mean a
sequence of one or more random variables with a fully specified joint
distribution. In this paper we consider independent and identically
distributed (i.i.d.) input processes that are themselves mutually independent,
which means that the input processes can be full specified by their individual
marginal distributions.

The distributions of these input processes are often estimated from real-world
data. Thus, a complete statistical characterization of stochastic system
performance requires quantifying both simulation and input estimation
errors.111As with any mathematical or computer model, simulations are also
subject to a host of non-statistical errors in abstraction; these are not
considered in this paper. There are robust methods that are adequate for
quantifying simulation error for many practical problems. However, the impact
of input estimation error (which we call “input uncertainty”) can overwhelm
the simulation error, as demonstrated in Barton Schruben (2001); ignoring it
may lead to unfounded confidence in the simulation assessment of system
performance, which could be the basis for critical and expensive decisions.
What we call input and simulation uncertainty are also known as epistemic and
aleatory uncertainty, respectively; see for instance Kleijnen (2008). Our
method accounts for both sources of error.

In this paper we address problems with parametric input distributions that are
mutually independent with parameters that are estimated from a finite sample
of real-world data. Of course, there exist practical problems in which the
input processes are not independent, and there may also be significant
uncertainty about the correct parametric distribution as well as its parameter
values. Nevertheless, the case of i.i.d. input processes represented by a
parametric distribution is prevalent in practice, making our contribution a
useful step. We build on Barton . (2014), which proposed a metamodel-assisted
bootstrapping approach that forms a confidence interval (CI) to account for
the impact of input uncertainty when estimating the system’s mean performance.
In that paper, bootstrap resampling of the real-world data was used to
approximate the input uncertainty, while a metamodel predicted the
simulation’s mean response at different parameter settings corresponding to
bootstrap resampled data sets.

In a metamodel-assisted bootstrapping framework there are three types of
error: the input estimation error, the simulation estimation error and the
error in the metamodel itself. The latter two types of error are not easily
separated, so we call their combined effect “metamodel uncertainty.”

Barton . (2014) showed that metamodel uncertainty can be ignored when the
simulation budget is not tight and an appropriate type of metamodel and
experiment design are used; they provided a follow-up test to insure that the
metamodel fit is close enough. In this setting their method yields a valid CI.
However, if the true mean response surface is complex, especially for high-
dimensional problems (i.e., many input distributions), and the computational
budget is tight, then the impact of metamodel uncertainty can no longer be
ignored without underestimating the error, which manifests itself in a CI that
is too short. Computationally intensive stochastic simulations are the norm
for a number of fields: Spatial stochastic simulations, e.g., of oil
reservoirs, can take hours for a single run, and depend on many stochastic
parameters (Bangerth . 2006, Wang . 2012). Simulations of semiconductor
manufacturing (Fowler Rose 2004) and biological systems (Ghosh . 2011, Kastner
. 2002) can be similarly intensive.

This paper is a significant enhancement of Barton . (2014). Here we propose an
approach to form an interval estimate that accounts for both input and
metamodel uncertainty in estimating a stochastic system’s mean performance.
When there is little metamodel uncertainty the new method performs like Barton
. (2014), but it does not experience a degradation in coverage when metamodel
uncertainty is significant.

When the statistical error measured by our CI is too large for the estimate to
be useful, then the decision maker may want to know how the error can be
reduced. Our approach leads naturally to a measure of the relative
contribution of input to overall uncertainty that indicates whether the error
can be reduced by an additional computational investment.

The next section describes other approaches to attack the input uncertainty
problem and contrasts them with our method. This is followed in Section 3 by a
formal description of the problem of interest and a brief review of the
metamodel-assisted bootstrapping approach in Section 4. In Section 5 we
provide an algorithm to build an interval estimator accounting for both input
and metamodel uncertainty, and give a method to estimate their relative
contributions. We then report results from an empirical study of a difficult
problem in Section 6 and conclude the paper in Section 7. All proofs are in
the Appendix.

##  2 Background

Various approaches to account for input uncertainty have been proposed. The
Bayesian methods use the posterior distributions of the inputs given the real-
world data to quantify the input distribution uncertainty, and the impact on
the system mean is estimated by drawing samples from these posterior
distributions and running simulations at each sample point (Chick 2001; Chick
Ng 2002; Zouaoui Wilson 2003, 2004). This could be computationally expensive
when the time for each simulation run is significant because simulations need
to be run at a large number of posterior sample points. In addition, for each
input prior distribution we need to derive a corresponding posterior
distribution which might be nonstandard and complex.

A second approach is based on direct bootstrapping; it quantifies the impact
of input uncertainty using bootstrap resampling of the input data, and runs
simulations at each bootstrap resample point to estimate the impact on the
system mean (Barton Schruben 1993, 2001; Barton 2007; Cheng Holland 1997).
Compared with the Bayesian approach, the direct bootstrap can be adapted to
any input process without additional analysis and it is suitable for complex
and nonstandard input distributions. However, similar to the Bayesian
approach, this method also runs simulations at each resample point. Since the
number of bootstrap resample points to construct a CI is recommended to be a
few thousand, the direct bootstrapping method is also computationally
expensive. More subtly, since the statistic that is bootstrapped is the random
output of a simulation it is not a smooth function of the input data; this
violates the asymptotic validity of the bootstrap.

Differing from the previous two approaches that estimate the system mean
response at each sample point by running simulations, a third approach
introduces an equation-based metamodel of the mean response (see Cheng Holland
(2004) and references therein). Specifically, it assumes that the parametric
families of the inputs are known, uses maximum likelihood estimators (MLEs) of
the unknown parameters, and represents input-parameter uncertainty by the
large-sample normal distribution of the MLEs. This uncertainty is propagated
to the output mean by a linear function of the parameters that is based on a
Taylor series approximation. Since the metamodel can be constructed using
simulation results from a small number of runs, this method does not need
substantial computational effort. However, a metamodel based on a locally
linear approximation is only appropriate when there is a large quantity of
real-world data so that the MLEs locate in a small neighborhood of the true
parameters with high probability; it is not suitable when the underlying
response surface is highly non-linear and only a modest quantity of real-world
data are available. In addition, the asymptotic normal approximation for the
input distribution parameters can be poor with sample sizes encountered in
some applications.

The metamodel-assisted bootstrapping approach introduced by Barton . (2014)
addresses some of the shortcomings in the prior work. Compared with Cheng
Holland (2004), the bootstrap provides a more accurate approximation of the
input uncertainty than the asymptotic normal distribution of the parameter
estimates in many situations (Horowitz, 2001). Further, the use of a general-
form metamodel provides higher fidelity than a locally linear approximation.
Compared with Bayesian and direct bootstrap approaches, the use of a metamodel
reduces the impact of simulation error on the accuracy of CIs and reduces the
computational effort because it does not run simulations at a large number of
sampled or resampled points; instead, an equation-based metamodel is
constructed based on a designed experiment at a smaller number of parameter
settings. In addition, employing a metamodel makes the bootstrap statistic a
smooth function of the input data so that the asymptotic validity concerns
faced by the direct bootstrap method disappear. The numerical results in
Barton . (2014) provide evidence that metamodel-assisted bootstrapping is
effective and superior to competitors when there is little metamodel
uncertainty, motivating its extension in this paper to more general and
complex input-uncertainty problems in which the impact of metamodel
uncertainty can no longer be ignored. The end result is a robust method for
quantifying statistical uncertainty.

##  3 Problem Description

To make the description of the input uncertainty problem and our solution to
it easier to follow we will use the queueing network in Figure 1 as an example
and return to it in our empirical study in Section 6. Consider estimating the
steady state expected number of customers in this network. The interarrival
times follow a gamma distribution, $A\sim\mbox{gamma}(\alpha_{A},\beta_{A})$,
and the service times at the $i$th station also follow a gamma distribution,
$S_{i}\sim\mbox{gamma}(\alpha_{S_{i}},\beta_{S_{i}})$. Customers finishing
service at stations $1,2,3$ must make decisions about their next station.
These routing decisions follow Bernoulli distributions
$P_{i}\sim\mbox{Ber}(p_{i}),i=1,2,3$. The parameters of the input
distributions, $\alpha_{A},\beta_{A}$,
$\\{(\alpha_{S_{i}},\beta_{S_{i}}),i=1,2,3,4\\}$ and $\\{p_{i},i=1,2,3\\}$ are
all unknown and estimated from real-world data. Notice that the inputs include
both continuous and discrete distributions. Our goal is to build a CI that
covers the steady-state expected number of customers in the system when the
input parameters assume their true but unknown values. We assume that at these
“true values” the system is in fact stable, and if we have enough real-world
data (which we may not) then the simulation with estimated parameters will
also be stable.

Figure 1: Queueing network example.

More generally, the stochastic simulation output is a function of random
numbers and $L$ independent input distributions
$F\equiv\\{F_{1},F_{2},\ldots,F_{L}\\}$. For notation simplification, we do
not explicitly include the random numbers. The output from the $j$th
replication of a simulation with input distribution $F$ can be written as

$Y_{j}(F)=\mu(F)+\epsilon_{j}(F)$ (1)

where $\mu(F)={\rm E}[Y_{j}(F)]$ denotes the unknown output mean and
$\epsilon_{j}(F)$ represents the simulation error with mean zero. Notice that
the simulation output depends on the choice of input distributions. The true
“correct” input distributions, denoted by
$F^{c}\equiv\\{F^{c}_{1},F^{c}_{2},\ldots,F^{c}_{L}\\}$, are unknown and are
estimated from a finite sample of real-world data. Our goal is to quantify the
impact of the statistical error by finding a $(1-\alpha)100\%$ CI
$[Q_{L},Q_{U}]$ such that

$\mbox{Pr}\\{\mu(F^{c})\in[Q_{L},Q_{U}]\\}=1-\alpha.$ (2)

This is a relatively general statement of the problem which could encompass
multivariate input distributions (i.e., if one or more of $F_{\ell}^{c}$ are
multivariate) and also unknown distribution families.

However, in this paper we assume that the input distributions are univariate,
the families are known, but the parameter values are not. Therefore, the input
distributions $F$ can be completely characterized by the collection of
parameters denoted by ${\bf x}$ with dimension $d$. We will define what we
mean by “parameters” below. By abusing the notation, we rewrite Equation (1)
as

$Y_{j}(\mathbf{x})=\mu(\mathbf{x})+\epsilon_{j}(\mathbf{x})$ (3)

where $\mathbf{x}\in\Psi$ and $\Psi\equiv\\{{\bf x}\in\Re^{d}{:}\mbox{ the
random variable }Y({\bf x})\mbox{ is defined and $\mu({\bf x})$ is finite}\\}$
denotes the region of interest. The collection of true parameters is denoted
by ${\bf x}_{c}$ and is assumed to be in the interior of $\Psi$. We also
assume $\mu({\bf x})$ is continuous for ${\bf x}\in\Psi$. Thus, our goal can
be restated as finding a $(1-\alpha)100\%$ CI such that

$\mbox{Pr}\\{\mu({\bf x}_{c})\in[Q_{L},Q_{U}]\\}=1-\alpha.$ (4)

Let $m_{\ell}$ denote the number of i.i.d. real-world observations available
from the $\ell$th input process
$\mathbf{Z}_{\ell,m_{\ell}}\equiv\left\\{Z_{\ell,1},Z_{\ell,2},\ldots,Z_{\ell,m_{\ell}}\right\\}$
with $Z_{\ell,i}\stackrel{{\scriptstyle i.i.d}}{{\sim}}F^{c}_{\ell}$,
$i=1,2,\ldots,m_{\ell}$. Let
$\mathbf{Z}_{\mathbf{m}}=\\{\mathbf{Z}_{\ell,m_{\ell}},\ell=1,2,\ldots,L\\}$
be the collection of samples from all $L$ input distributions in $F^{c}$,
where $\mathbf{m}=(m_{1},m_{2},\ldots,m_{L})$. The real-world data are a
particular realization of $\mathbf{Z}_{\mathbf{m}}$, say
$\mathbf{z}_{\mathbf{m}}^{(0)}$. Since the unknown input distributions are
estimated from $\mathbf{z}_{\mathbf{m}}^{(0)}$, we assume the parameters are
functions of $\mathbf{Z}_{\mathbf{m}}$ denoted by
$\mathbf{X}_{\mathbf{m}}=\mathbf{X}(\mathbf{Z}_{\mathbf{m}})$. Therefore,
input uncertainty is fully captured by the sampling distribution of
$\mu(\mathbf{X}_{\mathbf{m}})$.

Since the underlying response surface $\mu(\cdot)$ is unknown, we approximate
it by a metamodel fitted to a set of stochastic simulation outputs. Let
$\widehat{\mu}(\cdot)$ denote the metamodel. Fitting the metamodel introduces
a source of uncertainty in propagating the sampling distribution of
$\mathbf{X}_{\mathbf{m}}$ to the output mean: metamodel uncertainty. The
contributions of this paper are to construct an interval estimator that
accounts for both input and metamodel uncertainty, and to quantify the
relative impact of each contributor.

##  4 Metamodel-Assisted Bootstrapping Framework

Barton . (2014) introduced the metamodel-assisted bootstrapping approach to
quantify input uncertainty. We review it here.

The representation of the $L$ input distributions plays an important role in
the implementation of metamodel-assisted bootstrapping. Since this paper
focuses on problems with independent parametric distributions having unknown
parameters, $F$ can be uniquely characterized by the corresponding collection
of each distribution’s parameters or, in many cases, its moments. The $\ell$th
input distribution includes $h_{\ell}$ unknown parameters. Suppose that this
$h_{\ell}$-parameter distribution is uniquely specified by its first (finite)
$h_{\ell}$ moments, which is true for the distributions that are most often
used in stochastic simulation. The moments are chosen as the independent
variables for the metamodel because when they are close, the corresponding
distributions will be similar and therefore generate similar outputs. An
extended argument for choosing moments instead of the natural distribution
parameters as independent variables can be found in Barton . (2014). This
characterization will not work for all distributions, including some so-called
“heavy-tailed” distributions.

Let $\mathbf{x}_{[\ell]}$ denote an $h_{\ell}\times 1$ vector of the first
$h_{\ell}$ moments for the $\ell$th input distribution and
$d=\sum_{\ell=1}^{L}h_{\ell}$. By stacking ${\bf x}_{[\ell]}$ with
$\ell=1,2,\ldots,L$ together, we have a $d\times 1$ dimensional independent
variable ${\bf x}$. Notice that $F$ is completely characterized by the
collection of moments ${\bf x}$, and from here on $\mu(\cdot)$ will be a
function of distribution moments organized in this way. Denote the true
moments by ${\bf x}_{c}$.

For the queueing network in Figure 1, there are $L=8$ input distributions
including arrivals, services at stations $1,2,3,4$ and the three routing
processes. Thus, $h_{1}=h_{2}=h_{3}=h_{4}=h_{5}=2$ and $h_{6}=h_{7}=h_{8}=1$.
The distributions for the three routing processes can be completely specified
by their corresponding means. The distributions for the arrival and service
processes can be uniquely characterized by the corresponding first two
standardized moments: mean and standard deviation. To avoid a scaling effect,
the standard deviation is used instead of the second raw moment. For the
$\ell$th process, let $\tau_{\ell}$ denote the first moment (mean) and
$\varsigma_{\ell}$ denote the standard deviation. Then ${\bf
x}=(\tau_{1},\varsigma_{1},\tau_{2},\varsigma_{2},\tau_{3},\varsigma_{3},\tau_{4},\varsigma_{4},\tau_{5},\varsigma_{5},\tau_{6},\tau_{7},\tau_{8})^{\top}$
with $d=13$.

The true moments ${\bf x}_{c}$ are unknown and estimated based on a finite
sample $\mathbf{Z}_{\mathbf{m}}$ from $F^{c}$. As noted above,
$\mathbf{X}_{\mathbf{m}}$ is a $d\times 1$ dimensional moment estimator that
is a function of $\mathbf{Z}_{\mathbf{m}}$ written as
$\mathbf{X}_{\mathbf{m}}=\mathbf{X}(\mathbf{Z}_{\mathbf{m}})$. Specifically,
$\mathbf{X}_{\ell,m_{\ell}}=\mathbf{X}_{\ell}(\mathbf{Z}_{\ell,m_{\ell}})$ and
$\mathbf{X}_{\mathbf{m}}^{T}=(\mathbf{X}_{1,m_{1}}^{T},\mathbf{X}_{2,m_{2}}^{T},\ldots,\mathbf{X}_{L,m_{L}}^{T})$.
Let $F_{\mathbf{X}_{\mathbf{m}}}^{c}$ represent the true, unknown distribution
of $\mathbf{X}_{\mathbf{m}}$. Then “input uncertainty” refers to the
distribution of $\mu(\mathbf{X}_{\mathbf{m}})$ with
$\mathbf{X}_{\mathbf{m}}\sim F_{\mathbf{X}_{\mathbf{m}}}^{c}$. Given a finite
sample of real-world data $\mathbf{z}_{\mathbf{m}}^{(0)}$, Barton . (2014)
used bootstrap resampling to approximate $F_{\mathbf{X}_{\mathbf{m}}}^{c}$ and
a metamodel to represent $\mu({\bf x})$.

###  4.1 Bootstrap Resampling

We use distribution-free bootstrapping, meaning that we resample from the
empirical distribution of the data rather than a fitted parametric
distribution. Under some regularity conditions, the bootstrap can provide an
asymptotically consistent approximation for the sampling distribution of a
moment estimator (Shao Tu, 1995). For the general performance of the bootstrap
in representing the sampling distribution of an estimator, see Horowitz
(2001).

Implementation of the bootstrap in metamodel-assisted bootstrapping is as
follows.

  1. 1.

Draw $m_{\ell}$ samples with replacement from
$\mathbf{z}_{\ell,m_{\ell}}^{(0)}$, denoted by
$\mathbf{Z}_{\ell,m_{\ell}}^{(1)}$, and calculate the corresponding
$h_{\ell}\times 1$ vector of bootstrap moment estimates denoted by
$\widehat{{\mathbf{X}}}_{\ell,m_{\ell}}^{(1)}={{\bf
X}}_{\ell}(\mathbf{Z}_{\ell,m_{\ell}}^{(1)})$ for $\ell=1,2,\ldots,L$. Then
stack the results for all $L$ processes to obtain a $d\times 1$ vector
$\widehat{{\mathbf{X}}}_{\mathbf{m}}^{(1)}$.

  2. 2.

Repeat the previous step $B$ times to generate
$\widehat{{\mathbf{X}}}_{\mathbf{m}}^{(b)},b=1,2,\ldots,B$.

The bootstrap resampled moments are drawn from the bootstrap distribution
denoted by
$\widehat{F}_{\mathbf{X}_{\mathbf{m}}}(\cdot|\mathbf{z}_{\mathbf{m}}^{(0)})$.

The most straightforward approach to propagate input uncertainty to the output
mean is direct bootstrapping. Given a total computational budget of $N$
simulation replications, the system mean response at each of $B$ bootstrap
resamples is estimated by the sample mean of $n=\lfloor N/B\rfloor$
replications, denoted by $\bar{Y}(\widehat{{\bf
X}}_{\mathbf{m}}^{(b)})=\sum_{j=1}^{n}Y_{j}(\widehat{{\bf X}}_{\bf
m}^{(b)}),b=1,2,\ldots,B$. From these mean estimates a CI is formed. Direct
bootstrapping consumes the total simulation budget to estimate the simulation
mean responses and to form the CI. Thus, for direct bootstrapping the impact
of Monte Carlo estimation error on the CI and the impact of simulation
estimation error on the mean responses are not separable.

Barton . (2014) assume that there is enough computational budget available so
that the difference between $\widehat{\mu}({\bf x})$ and $\mu({\bf x})$ can be
ignored; i.e. $\widehat{\mu}({\bf x})\approx\mu({\bf x})$. Let
$\widehat{\mu}_{b}\equiv\widehat{\mu}(\widehat{\mathbf{X}}_{\mathbf{m}}^{(b)})$
for $b=1,2,\ldots,B$. Barton . (2014) quantified input uncertainty by a
$(1-\alpha)100\%$ bootstrap percentile interval induced by the sorted
responses
$[Q^{*}_{L},Q^{*}_{U}]=[\widehat{\mu}_{(\lceil\frac{\alpha}{2}B\rceil)},\widehat{\mu}_{(\lceil(1-\frac{\alpha}{2})B\rceil)}]$,
where $\widehat{\mu}_{(i)}$ is the $i$th smallest value of
$\widehat{\mu}_{1},\widehat{\mu}_{2},\ldots,\widehat{\mu}_{B}$. Here, the
superscript “*” indicates that the input distribution is approximated with the
bootstrap. This interval provides asymptotically correct coverage when all
input distributions meet certain moment conditions and the metamodel is
continuously differentiable with nonzero gradient in a neighborhood of
$\mathbf{x}_{c}$ (Barton ., 2014). Specifically, they establish the iterated
limit

$\lim_{m\rightarrow\infty}\lim_{B\rightarrow\infty}\mbox{Pr}\\{\mu({\bf
x}_{c})\in[Q^{*}_{L},Q^{*}_{U}]\\}=1-\alpha$

where as $m\rightarrow\infty$ we have $m_{\ell}/m\rightarrow 1$, for
$\ell=1,2,\ldots,L$.

However, with a tight computational budget for building the metamodel we
cannot guarantee that $\widehat{\mu}({\bf x})\approx\mu({\bf x})$ holds for
arbitrarily complex systems especially for problems with many input
distributions. Thus, we desire an interval estimator that accounts for both
input and metamodel uncertainty. Stochastic kriging (SK), introduced by
Ankenman . (2010), facilitates this. SK is flexible. Unlike the locally linear
approximation in Cheng Holland (2004), it does not require any strong
assumption about the form of the underlying true response surface
$\mu(\cdot)$. Based on our previous study (Xie ., 2010), a SK metamodel fit to
a small number of simulation runs can provide good global predictions and also
a characterization of metamodel uncertainty for a wide variety of examples.
The characterization of metamodel uncertainty is a key contribution of the new
method presented here.

Compared with direct bootstrapping, metamodel-assisted bootstrapping separates
the choice of $B$ from the budget $N$ and reduces the influence of simulation
estimation error. Instead doing simulations at $B$ samples from the bootstrap,
we run simulations at well-chosen design points and build an equation-based
metamodel $\widehat{\mu}({\bf x})$ to predict the mean response at different
input distributions represented by bootstrap resampled moments. Notice that
once we have the metamodel, we can use any $B$ we want to control the Monte
Carlo estimation error of the interval estimator, even $B>N$. Further, the
metamodel can efficiently use the computational budget to reduce the
uncertainty introduced when propagating the input uncertainty to the output
mean.

###  4.2 Stochastic Kriging Metamodel

Kriging is a widely used interpolation method. Since the outputs from
stochastic simulations include simulation variability that often changes
significantly across the design space, SK was introduced to distinguish the
uncertainty about the response surface from the simulation error inherent in
stochastic simulation output. This section provides a brief review of SK.

Suppose that the underlying true (but unknown) response surface can be thought
of as a realization of a stationary Gaussian Process (GP). This fiction has
been shown to provide a very useful framework for quantifying uncertainty
about the unknown surface implied by a deterministic computer experiment and
has been successfully employed in a wide variety of applications (see, for
instance, Santner . (2003)). SK extends this framework to include the
variability of the simulation output $Y$ using the model

$Y_{j}(\mathbf{x})=\mathbf{\beta}_{0}+W(\mathbf{x})+\epsilon_{j}(\mathbf{x}).$
(5)

The independent variable ${\bf x}$ is interpreted as a location in space. In
this paper, ${\bf x}$ denotes a $d\times 1$ vector of moments that uniquely
characterize the input distributions. The variation in the simulation output
over the $\mathbf{x}$ space is divided into extrinsic (response-surface)
uncertainty $W(\mathbf{x})$ and intrinsic (simulation output) uncertainty
$\epsilon_{j}(\mathbf{x})$. The term “intrinsic” uncertainty refers to the
variability inherent in the sampling that generates stochastic simulation
outputs and “extrinsic” uncertainty refers to our lack of knowledge about the
response surface.

SK uses a mean-zero, second-order stationary GP $W({\bf x})$ to account for
the spatial dependence of the response surface. Thus, the uncertainty about
the true response surface $\mu({\bf x})$ is represented by a GP $M({\bf
x})\equiv\mathbf{\beta}_{0}+W(\mathbf{x})$ (note that $\beta_{0}$ can be
replaced by a more general trend term $\mathbf{f}({\bf
x})^{\top}\mbox{\boldmath$\beta$}$ without affecting our method). For many,
but not all, simulation settings the output is an average of a large number of
more basic outputs, so a normal approximation can be applied: $\epsilon({\bf
x})\sim\mbox{N}(0,\sigma^{2}_{\epsilon}({\bf x}))$. Of course, normality of
the simulation output will not always hold, but could be empirically tested if
there is a concern.

In SK, the covariance between $W(\mathbf{x})$ and $W(\mathbf{x}^{\prime})$
quantifies how knowledge of the surface at some design points affects the
prediction of the surface. A parametric form of the spatial covariance,
denoted by
$\Sigma(\mathbf{x},\mathbf{x}^{\prime})=\mbox{Cov}[W(\mathbf{x}),W(\mathbf{x}^{\prime})]=\tau^{2}r(\mathbf{x}-\mathbf{x}^{\prime})$,
is typically assumed where $\tau^{2}$ denotes the variance and $r(\cdot)$ is a
correlation function that depends only on the distance
$\mathbf{x}-\mathbf{x}^{\prime}$. Based on our previous study (Xie ., 2010),
we use the product-form Gaussian correlation function

$r(\mathbf{x}-\mathbf{x}^{\prime})=\exp\bigg{(}-\sum_{j=1}^{d}\theta_{j}(x_{j}-x^{\prime}_{j})^{2}\bigg{)}$

for the empirical evaluation in Section 6; however, our results do not require
it. Let $\boldsymbol{\theta}=(\theta_{1},\theta_{2},\ldots,\theta_{d})$
represent the correlation parameters; for different correlation functions the
dimension of $\boldsymbol{\theta}$ could change. In any event, $M({\bf x})$
can be represented by a Gaussian process $M({\bf
x})\sim\mathrm{GP}(\mathbf{\beta}_{0},\tau^{2}r(\mathbf{x}-\mathbf{x}^{\prime})).$

To reduce the uncertainty about $\mu({\bf x})$ we choose an experiment design
consisting of pairs $\mathcal{D}\equiv\\{({\bf
x}_{i},n_{i}),i=1,2,\ldots,k\\}$ at which to run simulations and collect
observations, where $({\bf x}_{i},n_{i})$ denotes the location and the number
of replications, respectively, at the $i$th design point. The design that we
recommend is described in more detail in the Appendix, but it is not the only
design that could be effective. The simulation outputs at $\mathcal{D}$ are
$\mathbf{Y}_{\mathcal{D}}\equiv\left\\{(Y_{1}(\mathbf{x}_{i}),Y_{2}(\mathbf{x}_{i}),\ldots,Y_{n_{i}}(\mathbf{x}_{i}));i=1,2,\ldots,k\right\\}$
and the sample mean at design point ${\bf x}_{i}$ is
$\bar{Y}(\mathbf{x}_{i})=\sum_{j=1}^{n_{i}}Y_{j}(\mathbf{x}_{i})/n_{i}$. Let
the sample means at all $k$ design points be
$\bar{\mathbf{Y}}_{\mathcal{D}}=(\bar{Y}(\mathbf{x}_{1}),\bar{Y}(\mathbf{x}_{2}),\ldots,\bar{Y}(\mathbf{x}_{k}))^{T}$.
Since the use of common random numbers is detrimental to prediction (as
opposed to optimization; see Chen . (2012)), the simulations at different
design points are independent and the variance of
$\bar{\mathbf{Y}}_{\mathcal{D}}$ is represented by a $k\times k$ diagonal
matrix $C=\mbox{diag}\left\\{\sigma^{2}_{\epsilon}({\bf
x}_{1})/n_{1},\sigma^{2}_{\epsilon}({\bf
x}_{2})/n_{2},\ldots,\sigma^{2}_{\epsilon}({\bf x}_{k})/n_{k}\right\\}$.

Let $\Sigma$ be the $k\times k$ spatial covariance matrix of the design points
and let $\Sigma({\bf x},\cdot)$ be the $k\times 1$ spatial covariance vector
between each design point and a fixed prediction point ${\bf x}$. If the
parameters $(\tau^{2},\boldsymbol{\theta},C)$ are known, then the metamodel
uncertainty can be characterized by a refined GP $M_{p}({\bf x})$ that denotes
the conditional distribution of $M({\bf x})$ given all simulation outputs,

$M_{p}({\bf x})\sim\mathrm{GP}(m_{p}({\bf x}),\sigma^{2}_{p}({\bf x}))$ (6)

where $m_{p}(\cdot)$ is the minimum mean squared error (MSE) linear unbiased
predictor

$m_{p}({\bf x})=\widehat{\beta}_{0}+\Sigma({\bf
x},\cdot)^{\top}(\Sigma+C)^{-1}(\bar{\mathbf{Y}}_{\mathcal{D}}-\widehat{\beta}_{0}\cdot
1_{k\times 1}),$ (7)

and the corresponding variance is

$\sigma^{2}_{p}({\bf x})=\tau^{2}-\Sigma({\bf
x},\cdot)^{\top}(\Sigma+C)^{-1}\Sigma({\bf
x},\cdot)+\mathbf{\eta}^{\top}[1_{k\times 1}^{\top}(\Sigma+C)^{-1}1_{k\times
1}]^{-1}\mathbf{\eta}$ (8)

where $\widehat{\beta}_{0}=[1_{k\times 1}^{\top}(\Sigma+C)^{-1}1_{k\times
1}]^{-1}1_{k\times 1}^{\top}(\Sigma+C)^{-1}\bar{\mathbf{Y}}_{\mathcal{D}}$ and
$\mathbf{\eta}=1-1_{k\times 1}^{\top}(\Sigma+C)^{-1}\Sigma({\bf x},\cdot)$
(Ankenman ., 2010). With the parameters $(\tau^{2},\boldsymbol{\theta},{C})$
known, $M_{p}({\bf x})$ depends on the simulation outputs only through
$\bar{\mathbf{Y}}_{\mathcal{D}}$. Thus, $M_{p}({\bf x})$ is a random function
having the conditional distribution of $M({\bf x})$ given
$\bar{\mathbf{Y}}_{\mathcal{D}}$. Notice that $\sigma^{2}_{p}(\mathbf{x})$
reflects both metamodel and simulation error, including the constant term
$\widehat{\beta_{0}}$, with the intrinsic simulation sampling error affecting
$\sigma^{2}_{p}(\mathbf{x})$ through the matrix $C$.

Since in reality the spatial correlation parameters $\tau^{2}$ and
$\boldsymbol{\theta}$ are unknown, MLEs are typically used for prediction with
the log-likelihood function

$\ell(\beta_{0},\tau^{2},\boldsymbol{\theta})=-\ln[(2\pi)^{k/2}]-\frac{1}{2}\ln[|\Sigma+C|]-\frac{1}{2}(\bar{\mathbf{Y}}_{\mathcal{D}}-\beta_{0}\cdot
1_{k\times
1})^{\top}[\Sigma+C]^{-1}(\bar{\mathbf{Y}}_{\mathcal{D}}-\beta_{0}\cdot
1_{k\times 1})$

where $\Sigma$ is a function of $\tau^{2}$ and $\boldsymbol{\theta}$. The
sample variance is used as an estimate for the simulation variance at design
points $C$. By plugging
$(\widehat{\beta}_{0},\widehat{\tau}^{2},\widehat{\boldsymbol{\theta}},\widehat{C})$
into Equations (7) and (8) we can obtain the estimated mean
$\widehat{m}_{p}(\mathbf{x})$ and variance
$\widehat{\sigma}_{p}^{2}(\mathbf{x})$. Thus, the metamodel we use is
$\widehat{\mu}({\bf x})=\widehat{m}_{p}({\bf x})$ with marginal variance
estimated by $\widehat{\sigma}^{2}_{p}({\bf x})$.

Ankenman . (2010) demonstrate that $\widehat{m}_{p}({\bf x})$ is still an
unbiased predictor even with the plug-in estimator $\widehat{C}$, and further
that the variance inflation of $\sigma^{2}_{p}({\bf x})$ caused by using
$\widehat{C}$ is typically small. We performed an empirical study whose
results indicate that if we use an adequate experiment design, such as the
one-stage space-filling design used in this paper, then the performance of
metamodel-assisted bootstrapping is also not sensitive to the estimation error
in $\widehat{\tau}^{2}$ and $\widehat{\boldsymbol{\theta}}$; see the Appendix.
However, it is known that the estimator (8) with plug-in MLEs may sometimes
underestimate the prediction variance; see Den Hertog . (2006).

In the derivations that follow we will assume that the parameters
$(\tau^{2},\boldsymbol{\theta},C)$ are known. This is necessary (and common in
the kriging literature) because including the effect of parameter estimation
is mathematically intractable. To apply the methods in practice (including our
empirical study below), we form plug-in estimators by inserting
$\widehat{\tau}^{2},\widehat{\boldsymbol{\theta}},\widehat{C}$.

##  5 Confidence Interval and Variance Decomposition

Our approach is to use metamodel-assisted bootstrapping to provide a CI for
the true mean performance. To be robust the CI should account for both input
and metamodel uncertainty. Since $m_{p}({\bf x})$ is an unbiased predictor
under the Gaussian process assumption, $\sigma^{2}_{p}({\bf x})=0$ for all
${\bf x}$ would imply that there is no metamodel uncertainty due either to a
finite number of design points ${\bf x}_{i}$ or finite number of replications
$n_{i}$; that is, $m_{p}({\bf x})=\mu({\bf x})$. Unfortunately, with anything
short of complete information, there will always be some metamodel
uncertainty; and if the budget is tight relative to the complexity of the true
response surface, then the effect of metamodel uncertainty could be
substantial, resulting in significant undercoverage of the confidence interval
of Barton . (2014) as we show in Section 6. The new interval introduced here
does not suffer this degradation, and therefore is robust to the amount of
simulation effort that can be expended and can be recommended for general use.

The kriging literature is the foundation for our work; see for instance
Santner . (2003). Kriging provides inference about the value of an unknown
function $\mu(\cdot)$ at a fixed prediction point ${\bf x}_{0}$ where the
function has not been evaluated based on values of the function at a set of
design points. Kriging models uncertainty about the function as a GP
$M(\cdot)$ by assuming $\mu(\cdot)$ is a realization of $M(\cdot)$. An
interval constructed to cover the conditional distribution of $M({\bf x}_{0})$
given the values at the design points is often interpreted as a CI for
$\mu({\bf x}_{0})$ (e.g., Picheny . (2010)). The success of this paradigm is
not because the function of interest is actually random—it is not—but because
in many problems the conditional GP appears to be a robust characterization of
the remaining response-surface uncertainty.

We adopt the kriging paradigm but with two key differences: our prediction
point ${\bf x}_{c}$ is also unknown and must be estimated from real-world
data, and our function $\mu(\cdot)$ can only be evaluated in the presence of
stochastic simulation noise. Given the simulation outputs
$\bar{\mathbf{Y}}_{\mathcal{D}}$, the remaining uncertainty about $\mu(\cdot)$
is characterized by the conditional GP $M_{p}(\cdot)$. To account for the
impact from both input and metamodel uncertainty, we construct an interval
$[C_{L},C_{U}]$ covering $M_{p}({\bf x}_{c})$ with probability
$(1-\alpha)100\%$; that is,

$\mbox{Pr}\\{M_{p}({\bf x}_{c})\in[C_{L},C_{U}]\\}=1-\alpha.$ (9)

Since the conditional coverage is $1-\alpha$, the unconditional coverage of
$M({\bf x}_{c})$ is $1-\alpha$ as well. The revised objective (9) is connected
to our objective (4) through the assumption that the function $\mu(\cdot)$ is
a realization of the GP $M(\cdot)$. A procedure that delivers an interval
satisfying (9) will be a good approximation for a CI procedure satisfying (4)
if $M_{p}(\cdot)$ faithfully represents the remaining uncertainty about
$\mu(\cdot)$. This is clearly an approximation because in any real problem
$\mu(\cdot)$ is a fixed function, therefore we refer to $[C_{L},C_{U}]$ as an
approximation for the CI (ACI).

In a practical setting, what is the next step if the interval $[C_{L},C_{U}]$
is so wide that we are uncomfortable making decisions based on estimates with
that level of error? We suggest gaining some sense of the relative
contribution from each source of uncertainty as a guide toward either running
more simulations or collecting more real-world input data or both. For many
problems collecting additional input data is not feasible or we would have
done so already; in such cases knowing that input uncertainty is substantial
and cannot be reduced allows us to exercise caution in how we use the
simulation results.

In this section, we first present a procedure to build an ACI that satisfies
Equation (9) asymptotically. The asymptotic consistency of this interval is
proved under the assumption that the true response surface is a realization of
a GP with all parameters known except $\beta_{0}$. Next a variance
decomposition is proposed to measure the relative contribution of input
uncertainty to overall statistical uncertainty, and we study its asymptotic
properties as well. This is a measure of input uncertainty due to all input
distributions. A method for attributing the input uncertainty to the $L$
distributions is provided by Song Nelson (2013). Finally, we address problems
that can arise when the system’s mean performance fails to exist, or the
system is undefined, for some values of the sample moments, and explain why
metamodel-assisted bootstrapping tends to be tolerant of the former situation
and can be adjusted for the latter.

Assumptions that are needed for the asymptotic analysis are the following:

Assumptions:

  1. 1.

The $\ell$th input distribution is uniquely determined by its first $h_{\ell}$
moments and it has finite first $4h_{\ell}$ moments for $\ell=1,2,\ldots,L$.

  2. 2.

We have i.i.d observations
$Z_{\ell,1}^{(0)},Z_{\ell,2}^{(0)},\ldots,Z_{\ell,m_{\ell}}^{(0)}$ from the
$\ell$th distribution for $\ell=1,2,\ldots,L$. As $m\rightarrow\infty$, we
have $m_{\ell}/m\rightarrow c_{\ell}$, $\ell=1,2,\ldots,L$, for a constant
$c_{\ell}>0$.

  3. 3.

The $\epsilon_{j}(\mathbf{x})\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}\mbox{N}(0,\sigma^{2}_{\epsilon}(\mathbf{x}))$ for any ${\bf
x}$, and $M(\mathbf{x})$ is a stationary, separable GP with a continuous
correlation function satisfying

$1-r(\mathbf{x}-\mathbf{x}^{\prime})\leq\frac{c}{|\mbox{log}(\parallel\mathbf{x}-\mathbf{x}^{\prime}\parallel_{2})|^{1+\gamma}}\mbox{
for all }\parallel\mathbf{x}-\mathbf{x}^{\prime}\parallel_{2}\leq\delta$ (10)

for some $c>0$, $\gamma>0$ and $\delta<1$, where
$\parallel\mathbf{x}-\mathbf{x}^{\prime}\parallel_{2}=\sqrt{\sum_{j=1}^{d}({x}_{j}-{x}^{\prime}_{j})^{2}}.$

  4. 4.

The input processes $Z_{\ell j}^{(0)}$, simulation noise $\epsilon_{j}({\bf
x})$ and GP $M(\mathbf{x})$ are mutually independent and the bootstrap process
is independent of all of them.

Assumptions 1–2 give sufficient conditions for the almost sure (a.s.)
consistency of bootstrap moment estimators
$\widehat{\mathbf{X}}_{\mathbf{m}}\stackrel{{\scriptstyle
a.s.}}{{\rightarrow}}{\bf x}_{c}$ as $m\rightarrow\infty$ (see Lemma 1 in the
Appendix). Under Assumption 3, a GP $M(\cdot)$ with a correlation function
satisfying Condition (10) has continuous sample paths almost surely (Adler
2010, Theorem 3.4.1). Condition (10) is satisfied by many correlation
functions used in practice, and in particular any power exponential
correlation function $r({\bf x}-{\bf
x}^{\prime})=\exp\left(-\sum_{j=1}^{d}\theta_{j}|x_{j}-x_{j}^{\prime}|^{p}\right)$
with $0<p\leq 2$ and $\theta_{j}>0$ (Santner ., 2003). Assumption 4 indicates
that input data are collected independently of the simulation model, and that
our uncertainty about the mean response surface as represented by $M({\bf x})$
is independent of the stochastic simulation noise (although both can depend on
${\bf x}$).

###  5.1 ACI Procedure

Based on a hierarchical approach, we propose the following procedure to build
$(1-\alpha)100\%$ bootstrap percentile ACIs to achieve (9):

  1. 1.

Given real-world data $\mathbf{z}_{\mathbf{m}}^{(0)}$, choose experiment
design $\mathcal{D}=\\{({\bf x}_{i},n_{i}),i=1,2,\ldots,k\\}$ as described in
the Appendix.

  2. 2.

Run simulations at design points to obtain outputs $\mathbf{Y}_{\mathcal{D}}$.
Compute the sample average $\bar{Y}({\bf x}_{i})$ and sample variance
$S^{2}({\bf x}_{i})$ of the simulation outputs, $i=1,2,\ldots,k$. Fit the SK
metamodel parameters $(\beta_{0},\tau^{2},\boldsymbol{\theta},C)$ to obtain
$\widehat{m}_{p}({\bf x})$ and $\widehat{{\sigma}}_{p}^{2}(\mathbf{x})$ using
$\left(\bar{Y}({\bf x}_{i}),S^{2}({\bf x}_{i}),{\bf x}_{i}\right)$,
$i=1,2,\ldots,k$.

  3. 3.

For $b=1\mbox{ to }B$

    1. (a)

Generate bootstrap resample
$\mathbf{Z}^{(b)}_{\mathbf{m}}\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}\mathbf{z}^{(0)}_{\mathbf{m}}$ and compute sample moments
$\widehat{\mathbf{X}}_{\mathbf{m}}^{(b)}$.

    2. (b)

Let
$\widehat{\mu}_{b}\equiv\widehat{m}_{p}(\widehat{\mathbf{X}}_{\mathbf{m}}^{(b)})$.

    3. (c)

Draw
$\widehat{M}_{b}\sim\mbox{N}\left(\widehat{m}_{p}(\widehat{\mathbf{X}}_{\mathbf{m}}^{(b)}),\widehat{\sigma}_{p}^{2}(\widehat{\mathbf{X}}_{\mathbf{m}}^{(b)})\right)$.

  4. Next $b$

  5. 4.

Report estimated CI and ACI, respectively,

$\displaystyle\mbox{CI}_{0}$ $\displaystyle\equiv$
$\displaystyle\left[\widehat{\mu}_{(\lceil
B\frac{\alpha}{2}\rceil)},\widehat{\mu}_{(\lceil
B(1-\frac{\alpha}{2})\rceil)}\right]$ $\displaystyle\mbox{CI}_{+}$
$\displaystyle\equiv$ $\displaystyle\left[\widehat{M}_{(\lceil
B\frac{\alpha}{2}\rceil)},\widehat{M}_{(\lceil
B(1-\frac{\alpha}{2})\rceil)}\right]$

where
$\widehat{\mu}_{(1)}\leq\widehat{\mu}_{(2)}\leq\cdots\leq\widehat{\mu}_{(B)}$
and $\widehat{M}_{(1)}\leq\widehat{M}_{(2)}\leq\cdots\leq\widehat{M}_{(B)}$
are the sorted values.

In this procedure, Step 1 provides an experiment design to build a SK
metamodel, which is central to the metamodel-assisted bootstrapping approach.
Since the input uncertainty is quantified with bootstrap resampled moments, we
want the metamodel to correctly predict the responses at these points
$\widehat{\mathbf{X}}_{\mathbf{m}}\sim\widehat{F}_{\mathbf{X}_{\mathbf{m}}}(\cdot|\mathbf{z}_{\mathbf{m}}^{(0)})$.
Thus, the metamodel needs to be accurate and precise in a design space that
covers the “most likely” bootstrap moment estimates, which can be achieved by
the experiment design proposed by Barton . (2014). Their design is data-
driven; specifically, they first find the smallest ellipsoid denoted by $E$
that covers the most likely bootstrap resampled moments. They then generate a
space-filling design that covers $E$. This design methodology, which is
summarized in the Appendix, yielded accurate metamodels in the examples they
studied.

Based on the experiment design provided in Step 1, we run simulations and
construct a metamodel in Step 2 by fitting
$(\beta_{0},\tau^{2},\boldsymbol{\theta},C)$. Given the metamodel, we predict
the simulation’s mean responses at different input settings corresponding to
bootstrap resampled moments and construct interval estimators as shown in Step
3. Notice that Step 3(a) accounts for the input uncertainty and Step 3(c)
accounts for the input and metamodel uncertainty. Thus, this procedure
provides two types of intervals:

  * •

$\mbox{CI}_{0}$, proposed in Barton . (2014), returns an estimate of
$[Q_{L},Q_{U}]$ in Equation (4) by assuming $\widehat{m}_{p}({\bf x})=\mu({\bf
x})$; that is, it only accounts for input uncertainty and will be in error if
there is substantial metamodel uncertainty.

  * •

$\mbox{CI}_{+}$ returns an estimate of $[C_{L},C_{U}]$ in Equation (9). This
ACI accounts for both input and metamodel uncertainty.

As the metamodel uncertainty decreases, $\mbox{CI}_{0}$ and $\mbox{CI}_{+}$
become closer and closer to each other. Before evaluating the finite-sample
performance of $\mbox{CI}_{+}$ in Section 6, we establish its asymptotic
consistency for objective (9).

In Theorems 1–3 that follow, we replace $\widehat{\mu}_{b}$ and
$\widehat{M}_{b}$ in Steps 3(b)–(c) of the ACI procedure with

${\color[rgb]{0,0,0}\mu_{b}\equiv m_{p}(\widehat{{\bf X}}_{{\bf
m}}^{(b)})\mbox{ and
}}M_{b}\sim\mbox{N}\left(m_{p}(\widehat{\mathbf{X}}_{\mathbf{m}}^{(b)}),\sigma_{p}^{2}(\widehat{\mathbf{X}}_{\mathbf{m}}^{(b)})\right)$

recalling that $m_{p}(\cdot)$ and $\sigma_{p}^{2}(\cdot)$ are the marginal
mean and variance of the conditional GP $M_{p}(\cdot)$ when
$(\tau^{2},\boldsymbol{\theta},C)$ are known.

######  Theorem 5.1.

Suppose that Assumptions 1–4 hold. Then the interval $[M_{(\lceil
B\frac{\alpha}{2}\rceil)},M_{(\lceil B(1-\frac{\alpha}{2})\rceil)}]$ is
asymptotically consistent, meaning the iterated limit

$\lim_{m\rightarrow\infty}\lim_{B\rightarrow\infty}\Pr\\{M_{(\lceil
B\alpha/2\rceil)}\leq M_{p}(\mathbf{x}_{c})\leq M_{(\lceil
B(1-\alpha/2)\rceil)}\\}=1-\alpha.$ (11)

In brief, under the assumption that $\mu({\bf x})$ is a realization of a GP,
$M_{p}({\bf x})$ characterizes the remaining metamodel uncertainty after
observing $\bar{{\mathbf{Y}}}_{\mathcal{D}}$. And since the input uncertainty
is asymptotically correctly quantified by the bootstrap moment estimator
$\widehat{\mathbf{X}}_{\mathbf{m}}$, the distribution of
$M_{p}(\widehat{\mathbf{X}}_{\mathbf{m}})$ accounts for both input and
metamodel uncertainty. Theorem 5.1 shows that this interval satisfies
objective (9) asymptotically. We are particularly interested in situations
when the simulation effort is limited (Barton . (2014) addressed the ample
budget case), so the consistency result in Theorem 5.1 is only with respect to
the real-world data. The detailed proof is provided in the Appendix.

In practice, including our empirical evaluation in Section 6,
$(\tau^{2},\boldsymbol{\theta},C)$ must be estimated, and the impact of
parameter estimation (other than ${\bf x}_{c}$ and $\beta_{0}$) is not covered
by Theorem 1. We address sensitivity to parameter estimation in the Appendix.

###  5.2 Variance Decomposition

$\mbox{CI}_{+}$ accounts for input and metamodel uncertainty. When the width
of $\mbox{CI}_{+}$ is too large to be useful, it is important to know the
relative contribution from each source. Since the total output variability is
the convolution of the input uncertainty and simulation/metamodel uncertainty,
it is hard to separate the effects from these sources. To estimate the
relative contributions Zouaoui Wilson (2003), Ng Chick (2006), Ankenman Nelson
(2012) and Song Nelson (2013) assume that the simulation noise has a constant
variance. In this section, we propose a variance decomposition that does not
require the homogeneity assumption.

Suppose that the parameters $(\tau^{2},\boldsymbol{\theta},C)$ are known, the
metamodel uncertainty can be characterized by a GP and the simulation error
follows a normal distribution. Then the metamodel uncertainty, given the
simulation result $\bar{\mathbf{Y}}_{\mathcal{D}}$, is characterized by a GP
$M_{p}({\bf x})\sim\mbox{N}(m_{p}(\mathbf{x}),\sigma^{2}_{p}(\mathbf{x}))$.
Conditional on $\bar{\mathbf{Y}}_{\mathcal{D}}$, both $m_{p}(\mathbf{x})$ and
$\sigma^{2}_{p}(\mathbf{x})$ are fixed functions. For notation simplification,
all of following derivations are conditional on the simulation outputs
$\bar{\mathbf{Y}}_{\mathcal{D}}$, but we will suppress the
“$|\bar{\mathbf{Y}}_{\mathcal{D}}$”.

The random variable $M_{p}({\mathbf{X}}_{\mathbf{m}})$ accounts for input
uncertainty through the sampling distribution of $\mathbf{X}_{\mathbf{m}}$ and
the metamodel uncertainty through the random function $M_{p}(\cdot)$. To
quantify the relative contribution of input and metamodel uncertainty, we
decompose the total variance of $M_{p}({\mathbf{X}}_{\mathbf{m}})$ into two
parts:

$\displaystyle\sigma^{2}_{T}$ $\displaystyle\equiv$
$\displaystyle\mbox{Var}[M_{p}(\mathbf{X}_{\mathbf{m}})]$ (12)
$\displaystyle=$
$\displaystyle\mbox{E}\\{\mbox{Var}[M_{p}(\mathbf{X}_{\mathbf{m}})|\mathbf{X}_{\mathbf{m}}]\\}+\mbox{Var}\\{\mbox{E}[M_{p}(\mathbf{X}_{\mathbf{m}})|\mathbf{X}_{\mathbf{m}}]\\}$
$\displaystyle=$
$\displaystyle\mbox{E}[\sigma^{2}_{p}(\mathbf{X}_{\mathbf{m}})]+\mbox{Var}[m_{p}(\mathbf{X}_{\mathbf{m}})].$

The term
$\sigma^{2}_{M}\equiv\mbox{E}[\sigma^{2}_{p}(\mathbf{X}_{\mathbf{m}})]$ is a
measure of the metamodel uncertainty: the expected metamodel variance weighted
by the density of moment estimator $\mathbf{X}_{\mathbf{m}}$. This weighting
makes sense because the accuracy of the metamodel in regions with higher
density is more important for the estimation of system mean performance. The
term $\sigma^{2}_{I}\equiv\mbox{Var}[m_{p}(\mathbf{X}_{\mathbf{m}})]$ is a
measure of input uncertainty when we replace the unknown true response surface
$\mu(\cdot)$ with its best linear unbiased estimate $m_{p}(\cdot)$.

What is the contribution of each term to ACI coverage? If the metamodel
uncertainty disappears (i.e., $\sigma^{2}_{p}(\cdot)=0$), then
$\sigma^{2}_{M}=0$, $\mbox{CI}_{0}$ and $\mbox{CI}_{+}$ coincide and they
provide asymptotically consistent coverage (Barton ., 2014). Metamodel
uncertainty is reduced by simulation effort. On the other hand, as
$m\rightarrow\infty$ (more and more real-world input data),
$\mathbf{X}_{\mathbf{m}}\stackrel{{\scriptstyle a.s.}}{{\rightarrow}}{\bf
x}_{c}$ and since $m_{p}({\bf x})$ is continuous we have $\sigma^{2}_{I}=0$;
therefore, the width of $\mbox{CI}_{0}$ shrinks to zero as does coverage since
there is remaining metamodel uncertainty in general. However, because
$\mbox{CI}_{+}$ accounts for metamodel uncertainty it still provides
asymptotically consistent coverage. This effect is demonstrated by the
empirical study in Section 6.

Our decomposition allows us to express the total variance in Equation (12) as
the sum of two variances measuring input and metamodel uncertainty:
$\sigma^{2}_{T}=\sigma^{2}_{I}+\sigma^{2}_{M}.$ In the metamodel-assisted
bootstrapping framework, we can estimate each variance component as follows:

  * •

Total variance:
$\widehat{\sigma}^{2}_{T}=\sum_{b=1}^{B}(M_{b}-\bar{M})^{2}/(B-1)$, where
$\bar{M}=\sum_{b=1}^{B}M_{b}/B.$

  * •

Input variance:
$\widehat{\sigma}^{2}_{I}=\sum_{b=1}^{B}({\color[rgb]{0,0,0}\mu}_{b}-\bar{{\mu}})^{2}/(B-1)$,
where $\bar{{\mu}}=\sum_{b=1}^{B}{\color[rgb]{0,0,0}\mu}_{b}/B$.

  * •

Metamodel variance:
$\widehat{\sigma}^{2}_{M}=\sum_{b=1}^{B}{\color[rgb]{0,0,0}\sigma}_{p}^{2}(\widehat{\mathbf{X}}^{(b)}_{\mathbf{m}})/B$.

The ratio $\widehat{\sigma}_{I}/\widehat{\sigma}_{T}$ provides an estimate of
the relative contribution from input uncertainty on $\mbox{CI}_{+}$. If it is
close to 1, the contribution from metamodel uncertainty can be ignored. Thus,
this ratio can help a decision maker determine where to put more effort: If
the input variance dominates, then get more real-world data (if possible). If
the metamodel variance dominates, then it can be reduced by more simulation,
which can be a combination of additional design points and additional
replications at existing design points. If neither dominates, then both
activities are necessary to reduce $\mbox{CI}_{+}$ to a practically useful
size.

The asymptotic properties of these variance component estimators are shown in
the following theorems.

######  Theorem 5.2.

Suppose that Assumptions 1–4 hold. Then conditional on
$\bar{\mathbf{Y}}_{\mathcal{D}}$ the variance component estimators
$\widehat{\sigma}^{2}_{M},\widehat{\sigma}^{2}_{I},\widehat{\sigma}^{2}_{T}$
are consistent as $m,B\rightarrow\infty$, where as $m\rightarrow\infty$ we
have $m_{\ell}/m\rightarrow c_{\ell}$, $\ell=1,2,\ldots,L$, for a constant
$c_{\ell}>0$. Specifically,

  * •

As $m\rightarrow\infty$, the input uncertainty disappears:

$\lim_{m\rightarrow\infty}\sigma^{2}_{M}=\sigma^{2}_{p}({\bf
x}_{c}),\lim_{m\rightarrow\infty}\sigma^{2}_{I}=0\mbox{ and
}\lim_{m\rightarrow\infty}\sigma^{2}_{T}=\sigma^{2}_{p}({\bf x}_{c}).$

  * •

As $m\rightarrow\infty$ and $B\rightarrow\infty$ in an iterated limit, the
variance component estimators are consistent:

$\displaystyle\lim_{m\rightarrow\infty}\lim_{B\rightarrow\infty}\widehat{\sigma}_{M}^{2}$
$\displaystyle=$
$\displaystyle\lim_{m\rightarrow\infty}\sigma_{M}^{2}=\sigma_{p}^{2}({\bf
x}_{c}),$
$\displaystyle\lim_{m\rightarrow\infty}\lim_{B\rightarrow\infty}\widehat{\sigma}_{I}^{2}$
$\displaystyle=$ $\displaystyle\lim_{m\rightarrow\infty}\sigma_{I}^{2}=0,$
$\displaystyle\lim_{m\rightarrow\infty}\lim_{B\rightarrow\infty}\widehat{\sigma}_{T}^{2}$
$\displaystyle=$
$\displaystyle\lim_{m\rightarrow\infty}\sigma_{T}^{2}=\sigma^{2}_{p}({\bf
x}_{c}).$

Theorem 5.2 demonstrates that the variance components estimators
$\widehat{\sigma}^{2}_{I}$, $\widehat{\sigma}^{2}_{M}$ and
$\widehat{\sigma}^{2}_{T}$ are consistent. However, we can see that when
$m\rightarrow\infty$ the input uncertainty disappears. Since
$\lim_{m\rightarrow\infty}\sigma^{2}_{I}=\lim_{m\rightarrow\infty}\lim_{B\rightarrow\infty}\widehat{\sigma}^{2}_{I}=0$
is not interesting, we study the consistency of scaled versions of
$\sigma^{2}_{I}$ and $\widehat{\sigma}^{2}_{I}$ in Theorem 3, showing that
$m\sigma^{2}_{I}$ and $m\widehat{\sigma}_{I}^{2}$ converge to the same non-
zero constant.

######  Theorem 5.3.

Suppose that Assumptions 1–4 and the following additional assumptions hold:

  1. 5.

The first three derivatives of the correlation function of the GP $M({\bf x})$
exist and the third derivative is bounded; and

  2. 6.

$m_{\ell}/m\rightarrow 1$ for $\ell=1,2,\ldots,L$.

Then
$\lim_{m\rightarrow\infty}m\sigma_{I}^{2}=\lim_{m\rightarrow\infty}\lim_{B\rightarrow\infty}m\widehat{\sigma}^{2}_{I}=\sigma^{2}_{\mu}$
almost surely, where $\sigma^{2}_{\mu}$ is a positive constant.

Theorems 5.2–5.3 give the asymptotic properties of the variance component
estimators, guaranteeing $\widehat{\sigma}_{I}/\widehat{\sigma}_{T}$ is a
consistent estimator for the relative contribution of input to overall
uncertainty. We will empirically evaluate its finite-sample performance in
Section 6 where we form the variance component estimators by inserting
$(\widehat{\tau}^{2},\widehat{\boldsymbol{\theta}},\widehat{C})$ for the
unknown parameters $(\tau^{2},\boldsymbol{\theta},C)$.

###  5.3 Unstable and Undefined Moments

A fundamental assumption of simulation is that the expectation $\mu({\bf
x}_{c})$ exists. This assumption does not imply, however, that it exists for
all possible values of ${\bf x}$, ${\bf X}_{\bf m}$ or $\widehat{{\bf X}}_{\bf
m}^{(b)}$ that might be realized. The prototype example is a congestion-
related performance measure of a queueing system as time goes to infinity when
congestion increases without bound for some values of its interarrival-time
and service-time parameters. We refer to systems for which $\mu({\bf x})$ is
$\pm\infty$ for some values of ${\bf x}$ as potentially unstable.

Recall that $\mu({\bf x})={\rm E}[Y({\bf x})]$. A second problem arises when
for some values of ${\bf x}$ the random variable $Y({\bf x})$ is undefined.
The prototype example is a network for which we want to estimate some start-
to-finish performance measure, but the start and finish are not connected for
certain values of ${\bf x}$. We refer to such systems as potentially
undefined.

Below we use illustrative examples to describe what happens to metamodel-
assisted bootstrapping in each case, why we expect to be robust to unstable
systems, and what needs to be done for undefined systems. We assume that ${\bf
x}_{c}$ is an interior point of the space $\Psi$ for which $\mu({\bf x}_{c})$
is stable and $Y({\bf x}_{c})$ is defined so both problems disappear
asymptotically ($m\rightarrow\infty$), but they may occur when we apply the
metamodel-assisted bootstrapping approach to a finite sample of real-world
data.

####  5.3.1 Unstable Moments

Consider the simulation of an $M/M/1$ queue. Let $x_{1}$ and $x_{2}$ denote
the mean interarrival time and mean service time, respectively, and let ${\bf
x}=(x_{1},x_{2})^{\top}$. The unknown mean response $\mu({\bf x})$ is the
steady-state expected number of customers in the system. The true values
$x_{1}^{c}$ and $x_{2}^{c}$ are unknown and must be estimated from real-world
data; $x_{1}^{c}>x_{2}^{c}$ so the system is actually stable. We denote the
unstable and stable regions of ${\bf x}$ by $U=\\{(x_{1},x_{2}):0<x_{1}\leq
x_{2}\\}$ and $\bar{U}=\\{(x_{1},x_{2}):x_{1}>x_{2}>0\\}$, respectively, and
${\bf x}_{c}$ is an interior point of $\bar{U}$.

As described in the Appendix, we use an initial set of bootstrap resampled
moments to define an ellipsoid in which to embed our experiment design to fit
the metamodel, and then generate a second set at which we evaluate the
metamodel to form a two-sided, equal-probability bootstrap percentile
interval. The conditional probability that a bootstrap resampled moment
$\widehat{\mathbf{X}}_{\mathbf{m}}^{(b)}$ is located in the unstable region
given the real-world data is

$P_{U}\equiv\mbox{Pr}\left\\{\left.\widehat{\mathbf{X}}_{\mathbf{m}}^{(b)}\in
U\right|\mathbf{z}_{\mathbf{m}}^{(0)}\right\\}.$ (13)

For the $M/M/1$ queue we know $U$ so we know which bootstrap moments are in
the unstable region; therefore, we could map the mean response for unstable
moments (symbolically) to $\infty$ and only fit or use the metamodel to
predict the mean response at stable moments. If $P_{U}$ is large this could
lead to a one-sided lower confidence interval (infinite upper limit) that
would be interpreted as “the real system may be unstable;” this is an
appropriate conclusion if input uncertainty is substantial.

Unfortunately, in general stochastic systems it could be difficult or
impossible to determine which moments are in the unstable region either
analytically or empirically (Wieland ., 2003). Thus, in the experiment design
phase we might simulate the system at an ${\bf x}\in U$ to fit the metamodel,
and in the bootstrapping phase we might evaluate the resulting metamodel at an
${\bf x}\in U$ to estimate the CI. What is the effect of doing this when
$P_{U}>0$, possibly even large?

Suppose we start each replication of the $M/M/1$ queue with an empty and idle
system. Let $\mu({\bf x},t)$ denote the true expected number of customers in
the system at time $t\geq 0$. Then except for the case $x_{1}=x_{2}$, which we
ignore, it can be shown that $\mu({\bf x},t)$ satisfies the differential
equation

$\frac{d\mu({\bf x},t)}{dt}=\frac{1}{x_{1}}-\frac{1-p_{0}(t)}{x_{2}}$

where $p_{0}(t)$ is the probability that the system is empty at $t$ and
$p_{0}(0)=1$. If ${\bf x}\in\bar{U}$ then $p_{0}(t)\rightarrow 1-x_{2}/x_{1}$
as $t\rightarrow\infty$; however, if ${\bf x}\in U$ then $p_{0}(t)\rightarrow
0$ as $t\rightarrow\infty$. Thus, for large $t$,

$\frac{d\mu({\bf x},t)}{dt}\approx\left\\{\begin{array}[]{ll}0,&\mbox{if ${\bf
x}\in\bar{U}$}\\\ \frac{1}{x_{1}}-\frac{1}{x_{2}}>0,&\mbox{if ${\bf x}\in
U$}.\end{array}\right.$ (14)

For any finite run length $T$ and warm-up period $T_{0}<T$ the simulation
provides an unbiased estimator of

$\bar{\mu}({\bf x},T_{0},T)=\frac{1}{T-T_{0}}\int_{T_{0}}^{T}\mu({\bf
x},t)\,dt.$ (15)

Notice that this quantity is finite for any positive values of $x_{1}$ and
$x_{2}$, whether stable or not. However, if ${\bf x}\in\bar{U}$ then
$\bar{\mu}({\bf x},T_{0},T)$ converges to $\mu({\bf x})$ for large $T$; while
if ${\bf x}\in U$ then $\bar{\mu}({\bf x},T_{0},T)$ is increasing in $T$ for
$T$ large enough.

The key point is this: The expected value of any simulation-based estimator
will be finite, even if the true steady-state mean is not. Further, the
expected value of the simulation estimator at unstable ${\bf x}$ will tend to
be larger than at near-by stable ${\bf x}$. This means that the simulation
estimates corresponding to unstable ${\bf x}$ will tend to be the largest ones
observed, but still not infinite.

Consider the design points or bootstrap resampled moments that are in $U$.
When $P_{U}<\alpha/2$ and the run length is long enough, the unstable design
points used to fit the metamodel, or bootstrap moments at which it is
evaluated, tend not to adversely affect either the metamodel or the interval
estimate because they are in the right tail beyond the $\alpha/2$ quantile. On
the other hand, when $P_{U}\geq\alpha/2$ the large estimates corresponding to
unstable design points or bootstrap moments tend to lengthen the interval
estimate beyond what is required to cover $\mu({\bf x}_{c})$; this causes
overcoverage rather than undercoverage. Thus, metamodel-assisted bootstrapping
will often be robust to unstable moments in the sense of not being too short,
but possibly too long; we demonstrate this empirically in Section 6.

####  5.3.2 Undefined Moments

Consider the queueing network example in Figure 1. For simplicity, suppose
that the routing probabilities $p_{1},p_{2},p_{3}$ are the only input
parameters. Let ${\bf
x}=(x_{1},x_{2},x_{3})^{\top}=(p_{1},p_{2},p_{3})^{\top}$. The true parameters
${\bf x}_{c}$ are unknown and estimated by finite samples from Bernoulli
distributions ($1$ if the customer takes a particular route, $0$ otherwise).
Suppose that the mean response of interest, $\mu({\bf x}),$ is the steady-
state expected time for a customer to traverse the network, which exists and
is well-defined at ${\bf x}_{c}$. Unfortunately, $Y({\bf x})$ may not be
defined for every possible bootstrap resampled moment ${\bf x}$. For instance,
if $\widehat{{\bf X}}_{{\bf m}}^{(b)}=(0,0.665,0)^{\top}$ then Stations 1
(start) and 4 (end) are disconnected and no simulation output for time to
traverse the network will ever be generated. Thus, the system corresponding to
this bootstrap moment is undefined.

In practical problems for which we can obtain real-world input data, we should
know a priori that the system performance measure is well defined (e.g., we
would not include a route unless we actually observed a customer take it).
Further, it should not be difficult to detect moments for which the system
output is undefined, either because we understand the system logic (as in this
example) or because the simulation simply fails to run. Therefore, a
reasonable solution to the problem of undefined systems is to reject (and
sample again) bootstrap moments $\widehat{{\bf X}}_{{\bf m}}^{(b)}$ that imply
an undefined output. This makes our assessment of uncertainty conditional on
the system performance measure being defined, which makes sense.

##  6 Empirical Study

In this section we use the queueing network described in Section 3 to evaluate
the performance of our metamodel-assisted bootstrapping approach. The
performance measure is the steady-state expected number of customers in the
system. Both interarrival and service times follow gamma distributions and the
routing decisions follow Bernoulli distributions. Thus, it is a 13-dimensional
problem with $L=8$ input processes that include both continuous and discrete
distributions. The true parameters of the input distributions are
$\alpha_{A}=1,\beta_{A}=0.25$, $\alpha_{S_{i}}=1$, $\beta_{S_{i}}=0.2$ for
$i=1,2,3,4$ and $p_{1}=p_{2}=0.5,p_{3}=0.75$. These parameter values imply a
tractable Jackson network with steady-state number of customers in system
$\mu({\bf x}_{c})=12.67$. The maximum traffic intensity at any station is
$0.8$.

In the experiments we assume that all parameters for all input distributions
are unknown and are estimated from a finite sample of real-world data. Notice
that $\alpha_{A},\beta_{A}$, $\alpha_{S_{i}}$, $\beta_{S_{i}}$ for $i=1,2,3,4$
are estimated from continuous measurements, while the routing probabilities
$p_{1},p_{2},p_{3}$ are estimated from 0 or 1 observations that would
correspond to customer routing decisions. The model with estimated input
parameters is almost surely not a Jackson network and it could be unstable.
Our measure of uncertainty is a $95\%$ CI for $\mu({\bf x}_{c})$ as defined by
(4) because this is the objective desired in practice.

To evaluate the robustness of the metamodel-assisted bootstrapping approach,
we systematically examine the effect of the quantity of real-world data and
the number of design points and replications per design point used to fit the
metamodel; see Figure 2. We consider a wide range for the quantity of real-
world data $m=50,500,5000$, letting $m_{\ell}=m$ for $\ell=1,2,\ldots L$. The
levels for the number of design points are $k=20,40,80,130$. For a
13-dimensional problem $k=20$ is a very small design. The studies by Jones .
(1998) and Loeppky . (2009) recommend that the number of design points should
be $10$ times the dimension of the problem for kriging; we take this as the
maximum number of design points. The same number of replications are assigned
to all design points and we try $n=10,50,100$.

Figure 2: Experiment design.

Barton . (2014) demonstrated that $\mbox{CI}_{0}$ has good performance when
the impact of metamodel uncertainty is negligible. In this empirical study we
focus on situations where metamodel uncertainty may be significant. However,
rather than creating a problem that actually takes hours or days to run, we
instead construct a problem with high metamodel uncertainty by using short run
lengths for each replication: 20 time units after the warm up, which is
roughly equivalent to 80 finished customers. To avoid the influence from
initial bias, all simulations start loaded with the number of customers at
each station being their steady-state expected values (rounded) under ${\bf
x}_{c}$. Furthermore, a long warmup period of 200 time units is used. The net
effect is that the point estimators of the steady-state number in the network
have low bias, but may be quite variable.

To make the description of the empirical results easy to follow, we start with
overall conclusions:

  1. 1.

The new ACI $\mbox{CI}_{+}$ is robust to different levels of real-world data
$m$, number of design points $k$ and number of replications $n$.

  2. 2.

When metamodel uncertainty is significant, $\mbox{CI}_{0}$ tends to have
undercoverage that becomes more serious as $m$ increases. Since
$\mbox{CI}_{+}$ accounts for metamodel uncertainty, it does not exhibit this
degradation although it sometimes has slight overcoverage.

  3. 3.

Metamodel-assisted bootstrapping continues to deliver at least the nominal
coverage when the probability of an unstable system $P_{U}$ is large.

  4. 4.

The ratio $\widehat{\sigma}_{I}/\widehat{\sigma}_{T}$ is a useful measure of
the relative contribution of input to overall statistical uncertainty.

As discussed in Section 5.3.1, metamodel-assisted bootstrapping might behave
differently when $P_{U}<\alpha/2$ vs. $P_{U}\geq\alpha/2$. Since $P_{U}$ only
depends on $m$ and ${\bf x}_{c}$, we ran a side experiment to estimate it
using

$\widehat{P}_{U}=\frac{1}{B}\sum_{b=1}^{B}\mbox{I}\left(\widehat{\mathbf{X}}_{\mathbf{m}}^{(b)}\in
U\right),$ (16)

where $\mbox{I}(\cdot)$ is the indicator function. The means and standard
deviations (SD) of $\widehat{P}_{U}$ for $m=50,500,5000$ were estimated based
on 1000 macro-replications and are displayed in Table 1. In each macro-
replication we independently generated a sample of size $m$ of “real-world
data.” Then, conditional on these data, we drew $B=2000$ bootstrap resampled
moments. Finally, we calculated the estimate of $\widehat{P}_{U}$ using
Equation (16).

As $m$ increases the bootstrap resampled moments become more closely centered
around ${\bf x}_{c}$. Thus, both the mean and SD of $\widehat{P}_{U}$ decrease
with increasing $m$ as shown in Table 1. When $m=50$, ${P}_{U}$ appears to be
much larger than $\alpha/2$ so the bootstrap moments
$\widehat{\mathbf{X}}_{\mathbf{m}}^{(b)}$ that correspond to the upper
confidence bound are located in the unstable region $U$ with high probability.
When $m=500$, ${P}_{U}$ appears to be close to $\alpha/2=2.5\%$, while when
$m=5000$ there is little chance of getting unstable bootstrap moments.

Table 1: Percentage of unstable bootstrap resampled moments. | $m=50$ | $m=500$ | $m=5000$
---|---|---|---
mean of $\widehat{P}_{U}$ | 44.4% | 2.3% | 0
SD of $\widehat{P}_{U}$ | 31.7% | 7.9% | 0

In the following sections we describe the overall performance of
$\mbox{CI}_{0}$ and $\mbox{CI}_{+}$, including the situation where $P_{U}>0$,
and analyze the finite-sample performance of
$\widehat{\sigma}_{I}/\widehat{\sigma}_{T}$ as a measure of the relative
contribution of input to overall uncertainty.

###  6.1 Performance of CIs

Tables 2–3 show the results for $\mbox{CI}_{0}$ and $\mbox{CI}_{+}$ when
$m=50,500,5000$, including the probability of covering $\mu({\bf x}_{c})$, and
the mean and SD of the interval widths. All results are based on $1000$ macro-
replications.

When $m=50$, $P_{U}$ is much greater than $\alpha/2$ according to Table 1.
This explains the very large CI widths in Table 2. Nevertheless, both
$\mbox{CI}_{0}$ and $\mbox{CI}_{+}$ have reasonable coverage overall, an
observation we explore further in Section 6.3. Notice that $\mbox{CI}_{0}$
does exhibit undercoverage when we use a very small experiment design of
$k=20$ points, while the coverage of $\mbox{CI}_{+}$ is much closer to the
nominal value of $95\%$ in this case. If we fix the number of replications $n$
and increase the number of design points $k$, the coverage of $\mbox{CI}_{0}$
improves. For a fixed $k$ the effect of increasing $n$ is not as obvious.

Table 3 shows the results for $m=500,5000$. Compared with the results for
$m=50$, the mean and SD of the interval widths drop dramatically. The effects
of $k$ and $n$ are easier to discern especially when $m=5000$, which has no
unstable bootstrap moments. Specifically, for a fixed quantity of real-world
data $m$, if either the number of design points $k$ or replications per design
point $n$ is small then $\mbox{CI}_{0}$ tends to have undercoverage because it
fails to account for substantial metamodel uncertainty, unlike
$\mbox{CI}_{+}$. However, because $\mbox{CI}_{+}$ does incorporate metamodel
uncertainty it sometimes has slight overcoverage.

The most troubling observation about $\mbox{CI}_{0}$ is that, for fixed
$(n,k)$, as the amount of input data $m$ increases its undercoverage becomes
more serious. The diminished coverage occurs because as $m\rightarrow\infty$
the width of $\mbox{CI}_{0}$ shrinks to zero, which is not appropriate when
there is still metamodel uncertainty. Again, $\mbox{CI}_{+}$ does not exhibit
this degradation. As $n$ and $k$ increase, the coverages of $\mbox{CI}_{0}$
and $\mbox{CI}_{+}$ become closer to each other.

The behavior of $\mbox{CI}_{0}$ is what we would expect based on Barton .
(2014), which introduced $\mbox{CI}_{0}$. Their procedure continued to add
simulation effort (design points and replications) until its effect on the
confidence interval was negligible. Compared to $\mbox{CI}_{0}$, the new
interval, $\mbox{CI}_{+}$, is able to account for the effect of the remaining
simulation estimation error. Therefore, it can work under more general
situations where the simulated systems are complex and the simulation budget
is tight.

Table 2: Results for $\mbox{CI}_{0}$, $\mbox{CI}_{+}$ and $\widehat{\sigma}_{I}/\widehat{\sigma}_{T}$ when $m=50$.  $m=50$ | $k=20$ | $k=40$
---|---|---
| $n=10$ | $n=50$ | $n=100$ | $n=10$ | $n=50$ | $n=100$
Coverage of $\mbox{CI}_{0}$ | 91.9% | 92.3% | 91.5% | 93.8% | 94.4% | 93.4%
Coverage of $\mbox{CI}_{+}$ | 93.9% | 94.9% | 93.7% | 94.9% | 95.6% | 95.9%
$\mbox{CI}_{0}$ Width (mean) | 326.4 | 332.4 | 339.5 | 319.1 | 328.6 | 326.5
$\mbox{CI}_{+}$ Width (mean) | 344.1 | 348.8 | 357.1 | 332.3 | 342.3 | 341.2
$\mbox{CI}_{0}$ Width (SD) | 183.1 | 173.6 | 180.7 | 176.4 | 167.6 | 175
$\mbox{CI}_{+}$ Width (SD) | 188 | 175.7 | 183.8 | 178.2 | 169.2 | 176.1
$\widehat{\sigma}_{I}/\widehat{\sigma}_{T}$ | 0.963 | 0.965 | 0.964 | 0.973 | 0.973 | 0.971
$m=50$ | $k=80$ | $k=130$
| $n=10$ | $n=50$ | $n=100$ | $n=10$ | $n=50$ | $n=100$
Coverage of $\mbox{CI}_{0}$ | 94.6% | 96.3% | 95.4% | 94.2% | 95.1% | 95.4%
Coverage of $\mbox{CI}_{+}$ | 95.9% | 96.7% | 96.1% | 94.5% | 96% | 96.1%
$\mbox{CI}_{0}$ Width (mean) | 312.1 | 314.8 | 322.7 | 322 | 321.86 | 320
$\mbox{CI}_{+}$ Width (mean) | 322 | 325.7 | 334 | 330.2 | 331 | 329.4
$\mbox{CI}_{0}$ Width (SD) | 169.7 | 159.1 | 164.7 | 171.5 | 169.3 | 172.3
$\mbox{CI}_{+}$ Width (SD) | 171.2 | 159.4 | 165 | 172.7 | 169.5 | 172.7
$\widehat{\sigma}_{I}/\widehat{\sigma}_{T}$ | 0.982 | 0.98 | 0.978 | 0.985 | 0.985 | 0.983
Table 3: Results for $\mbox{CI}_{0}$, $\mbox{CI}_{+}$ and $\widehat{\sigma}_{I}/\widehat{\sigma}_{T}$ when $m=500$ and $m=5000$.  $m=500$ | $k=20$ | $k=40$
---|---|---
| $n=10$ | $n=50$ | $n=100$ | $n=10$ | $n=50$ | $n=100$
Coverage of $\mbox{CI}_{0}$ | 90.5% | 94.6% | 95.1% | 94.9% | 96.7% | 96.4%
Coverage of $\mbox{CI}_{+}$ | 95.7% | 97.7% | 97.8% | 96.6% | 98.3% | 97.8%
$\mbox{CI}_{0}$ Width (mean) | 24.8 | 28.1 | 29.4 | 27.1 | 28.5 | 28.7
$\mbox{CI}_{+}$ Width (mean) | 28.9 | 30.8 | 32.2 | 29.6 | 30.3 | 30.5
$\mbox{CI}_{0}$ Width (SD) | 19.9 | 19.4 | 20.6 | 19.1 | 19.2 | 19.9
$\mbox{CI}_{+}$ Width (SD) | 20.6 | 20.4 | 21.7 | 19.7 | 19.9 | 20.6
$\widehat{\sigma}_{I}/\widehat{\sigma}_{T}$ | 0.88 | 0.932 | 0.933 | 0.932 | 0.957 | 0.958
$m=500$ | $k=80$ | $k=130$
| $n=10$ | $n=50$ | $n=100$ | $n=10$ | $n=50$ | $n=100$
Coverage of $\mbox{CI}_{0}$ | 96.5% | 97.5% | 95.8% | 95.4% | 96.5% | 95.9%
Coverage of $\mbox{CI}_{+}$ | 98% | 98.3% | 97.3% | 97.5% | 97.1% | 96.9%
$\mbox{CI}_{0}$ Width (mean) | 26.3 | 28 | 28.7 | 26.4 | 27.9 | 27.6
$\mbox{CI}_{+}$ Width (mean) | 28 | 29 | 29.7 | 27.9 | 28.6 | 28.2
$\mbox{CI}_{0}$ Width (SD) | 17.4 | 18 | 19.3 | 18.8 | 19.6 | 19.3
$\mbox{CI}_{+}$ Width (SD) | 17.7 | 18.4 | 19.6 | 18.9 | 19.9 | 19.5
$\widehat{\sigma}_{I}/\widehat{\sigma}_{T}$ | 0.952 | 0.977 | 0.978 | 0.957 | 0.984 | 0.987
$m=5000$ | $k=20$ | $k=40$
| $n=10$ | $n=50$ | $n=100$ | $n=10$ | $n=50$ | $n=100$
Coverage of $\mbox{CI}_{0}$ | 70.7% | 89.2% | 93.1% | 81.5% | 94.3% | 94.8%
Coverage of $\mbox{CI}_{+}$ | 91.3% | 96.3% | 95.6% | 96.5% | 96.1% | 96.3%
$\mbox{CI}_{0}$ Width (mean) | 3.29 | 3.97 | 4.14 | 3.93 | 4.23 | 4.3
$\mbox{CI}_{+}$ Width (mean) | 5.85 | 4.8 | 4.56 | 6.08 | 4.64 | 4.52
$\mbox{CI}_{0}$ Width (SD) | 1.89 | 1.2 | 1 | 1.64 | 0.87 | 0.83
$\mbox{CI}_{+}$ Width (SD) | 2.12 | 1.13 | 1 | 1.52 | 0.89 | 0.85
$\widehat{\sigma}_{I}/\widehat{\sigma}_{T}$ | 0.588 | 0.85 | 0.924 | 0.664 | 0.924 | 0.959
$m=5000$ | $k=80$ | $k=130$
| $n=10$ | $n=50$ | $n=100$ | $n=10$ | $n=50$ | $n=100$
Coverage of $\mbox{CI}_{0}$ | 88.9% | 93.6% | 94.9% | 89.5% | 93.7% | 94.8%
Coverage of $\mbox{CI}_{+}$ | 98.1% | 95% | 96% | 98% | 95.6% | 95.5%
$\mbox{CI}_{0}$ Width (mean) | 4.54 | 4.29 | 4.29 | 4.52 | 4.35 | 4.32
$\mbox{CI}_{+}$ Width (mean) | 6.1 | 4.56 | 4.42 | 5.98 | 4.64 | 4.45
$\mbox{CI}_{0}$ Width (SD) | 1.37 | 0.85 | 0.77 | 1.28 | 0.9 | 0.79
$\mbox{CI}_{+}$ Width (SD) | 1.27 | 0.85 | 0.78 | 1.13 | 0.87 | 0.77
$\widehat{\sigma}_{I}/\widehat{\sigma}_{T}$ | 0.757 | 0.946 | 0.974 | 0.766 | 0.945 | 0.974

###  6.2 Performance of $\widehat{\sigma}_{I}/\widehat{\sigma}_{T}$

Tables 2–3 also demonstrate that $\widehat{\sigma}_{I}/\widehat{\sigma}_{T}$
provides a good measure of the relative contribution of input to overall
uncertainty, and behaves as it should:

  * •

For a fixed amount of real-world data $m$, increasing the number of design
points and replications $(n,k)$ drives
$\widehat{\sigma}_{I}/\widehat{\sigma}_{T}$ toward 1, indicating a decrease in
metamodel uncertainty.

  * •

For fixed simulation effort $(n,k)$, increasing the amount of real-world data
$m$ decreases $\widehat{\sigma}_{I}/\widehat{\sigma}_{T}$, indicating that
there is relatively less input uncertainty. Notice, however, that the
relationship is not simple because as $m$ increases the design space over
which we fit the metamodel becomes smaller, so that even with the same
simulation effort the absolute level of metamodel uncertainty will decrease
somewhat.

  * •

When $\widehat{\sigma}_{I}/\widehat{\sigma}_{T}$ is near $1$, the behaviors
(coverage and width) of $\mbox{CI}_{0}$ and $\mbox{CI}_{+}$ are similar and
both have coverage close to the nominal level; this is illustrated in Figure
3. Recall that $\mbox{CI}_{0}$ does not account for metamodel uncertainty, and
that $\widehat{\sigma}_{I}/\widehat{\sigma}_{T}\approx 1$ indicates that input
uncertainty is large relative to metamodel uncertainty, which is when
$\mbox{CI}_{0}$ will do best. Figure 3 also illustrates the general robustness
of $\mbox{CI}_{+}$.

Figure 3: The coverage errors for $\mbox{CI}_{0}$ and $\mbox{CI}_{+}$ vs.
$\widehat{\sigma}_{I}/\widehat{\sigma}_{T}$ when $m=5000$ across all values of
$n$ and $k$.

###  6.3 Robustness to Unstable Moments

Recall the observation from Table 2 that when there is a small quantity of
real-world data ($m=50$), resulting in a large probability of unstable
bootstrap moments, then both $\mbox{CI}_{0}$ and $\mbox{CI}_{+}$ had large
mean and SD of their widths, yet provided reasonable coverage. Examining the
results, we found that most of the intervals that fail to cover the mean do so
because the lower confidence bound is above $\mu({\bf x}_{c})$; this is the
case for both $\mbox{CI}_{0}$ and $\mbox{CI}_{+}$. Using all 1000 macro-
replications, the estimated probability that the lower confidence bound is
above the mean (greater than $\mu({\bf x}_{c})$, meaning too large) is 4.4%
for $\mbox{CI}_{0}$ and 3.8% for $\mbox{CI}_{+}$, while the estimated
probability that the upper confidence bound is below the mean (less than
$\mu({\bf x}_{c})$, meaning too small) is only $0.2\%$ for $\mbox{CI}_{0}$ and
$0.1\%$ for $\mbox{CI}_{+}$; for two-sided equal percentile intervals we would
expect these to be around $2.5\%$.

We conclude that even though the metamodel predicts a finite mean when it
should be infinite, this will still tend to lead to overcoverage rather than
undercoverage, and therefore is conservative.

##  7 Conclusions

In this paper, a metamodel-assisted bootstrapping approach is used for
statistical uncertainty analysis. Input uncertainty is approximated by the
bootstrap, an equation-based stochastic kriging metamodel is used to propagate
the input uncertainty to the output mean, and the metamodel uncertainty is
derived using properties of stochastic kriging. This approach delivers an
interval estimator that accounts for all statistical uncertainty, both
simulation and input. The asymptotic consistency of this interval is proved
under the assumption that the true response surface is a realization of a
Gaussian process and certain parameters are known.

An empirical study on a difficult problem demonstrates that our approach can
have good finite-sample performance even when there are several input
distributions (both discrete-valued and continuous-valued), a tight
computational budget, and bootstrap moments corresponding to unstable systems.
Thus, the new interval does not require a sequential experiment to make
metamodel uncertainty negligible, as in Barton . (2014).

If $\mbox{CI}_{+}$ is too wide, then it is important to know the relative
contributions from input and metamodel uncertainty as a guide to either
collecting more real-world data or doing more simulation or both. We give a
measure of the relative contribution of input to overall statistical
uncertainty by using a variance decomposition and analyze its asymptotic
properties.

## Acknowledgments

This paper is based upon work supported by the National Science Foundation
under Grant Nos. CMMI-0900354 and CMMI-1068473. The authors thank Ohad Perry,
the associate editor and two anonymous referees for their help in the
presentation of results and certain technical issues.

## References

  * Adler (2010) Adler_2010Adler, RJ.  2010\.  The Geometry of Random Fields The geometry of random fields.  Philadelphia, PASIAM.
  * Ankenman Nelson (2012) Ankenman_Nelson_2012Ankenman, BE. Nelson, BL.  2012\.  A Quick Assessment of Input Uncertainty A quick assessment of input uncertainty.  Proceedings of the 2012 Winter Simulation Conference Proceedings of the 2012 winter simulation conference ( 241-250).  Piscataway, New Jersey: Institute of Electrical and Electronics Engineers, Inc.
  * Ankenman . (2010) ankenman_nelson_staum_2010Ankenman, BE., Nelson, BL. Staum, J.  2010\.  Stochastic Kriging for Simulation Metamodeling Stochastic kriging for simulation metamodeling.  Operations Research58371-382.
  * Bangerth . (2006) Bangerth_2006Bangerth, W., Klie, H., Wheeler, MF., Stoffa, PL. Sen, MK.  2006\.  On Optimization Algorithms for the Reservoir Oil Well Placement Problem On optimization algorithms for the reservoir oil well placement problem.  Computational Geosciences10303-319.
  * Barton (2007) Barton_2007Barton, RR.  2007\.  Presenting A More Complete Characterization of Uncertainty: Can It Be Done? Presenting a more complete characterization of uncertainty: Can it be done?  Proceedings of the 2007 INFORMS Simulation Society Research Workshop. Proceedings of the 2007 informs simulation society research workshop.  Fontainebleau.
  * Barton . (2014) barton_nelson_xie_2011Barton, RR., Nelson, BL. Xie, W.  2014\.  Quantifying Input Uncertainty via Simulation Confidence Interval Quantifying input uncertainty via simulation confidence interval.  Informs Journal on Computing2674-87.
  * Barton Schruben (1993) Barton_Schruben_1993Barton, RR. Schruben, LW.  1993\.  Uniform And Bootstrap Resampling of Input Distributions Uniform and bootstrap resampling of input distributions.  Proceedings of the 1993 Winter Simulation Conference Proceedings of the 1993 winter simulation conference ( 503-508).  Piscataway, New Jersey: Institute of Electrical and Electronics Engineers, Inc.
  * Barton Schruben (2001) Barton_Schruben_2001Barton, RR. Schruben, LW.  2001\.  Resampling Methods for Input Modeling Resampling methods for input modeling.  Proceedings of the 2001 Winter Simulation Conference Proceedings of the 2001 winter simulation conference ( 372-378).  Piscataway, New Jersey: Institute of Electrical and Electronics Engineers, Inc.
  * Billingsley (1995) Billingsley_1995Billingsley, P.  1995\.  Probability and Measure Probability and measure.  New YorkWiley-Interscience.
  * Chen . (2012) chen_ankenman_nelson_2010Chen, X., Ankenman, BE. Nelson, BL.  2012\.  The Effects of Common Random Numbers on Stochastic Kriging Metamodels The effects of common random numbers on stochastic kriging metamodels.  ACM Transactions on Modeling and Computer Simulation227.
  * Cheng Holland (1997) Cheng_Holland_1997Cheng, RCH. Holland, W.  1997\.  Sensitivity of Computer Simulation Experiments to Errors in Input Data Sensitivity of computer simulation experiments to errors in input data.  Journal of Statistical Computation and Simulation57219-241.
  * Cheng Holland (2004) Cheng_Holland_2004Cheng, RCH. Holland, W.  2004\.  Calculation of Confidence Intervals for Simulation Output Calculation of confidence intervals for simulation output.  ACM Transactions on Modeling and Computer Simulation14344-362.
  * Chick (2001) Chick_2001Chick, SE.  2001\.  Input Distribution Selection for Simulation Experiments: Accounting for Input Uncertainty Input distribution selection for simulation experiments: Accounting for input uncertainty.  Operations Research49744-758.
  * Chick Ng (2002) Chick_Ng_2002Chick, SE. Ng, SH.  2002\.  Joint Criterion for Factor Identification and Parameter Estimation Joint criterion for factor identification and parameter estimation.  Proceedings of the 2002 Winter Simulation Conference Proceedings of the 2002 winter simulation conference ( 400-406).  Piscataway, New Jersey: Institute of Electrical and Electronics Engineers, Inc.
  * Den Hertog . (2006) DenHertog_2006Den Hertog, D., Kleijnen, JPC. Siem, AYD.  2006\.  The correct Kriging variance estimated by bootstrapping The correct Kriging variance estimated by bootstrapping.  Journal of the Operational Research Society57400-409.
  * Fowler Rose (2004) fowler2004grandFowler, JW. Rose, O.  2004\.  Grand challenges in modeling and simulation of complex manufacturing systems Grand challenges in modeling and simulation of complex manufacturing systems.  Simulation80469–476.
  * Ghosh . (2011) Ghosh_2011Ghosh, S., Matsuoka, Y., Asai, Y., Hsin, KY. Kitano, H.  2011\.  Software for Systems Biology: From Tools to Integrated Platforms Software for systems biology: From tools to integrated platforms.  Nature Reviews Genetics12821-832.
  * Horowitz (2001) Horowitz_2001Horowitz, JL.  2001\.  The Bootstrap The bootstrap ( 5).  Oxford, UKNorth Holland.
  * Jones . (1998) Jones_1998Jones, D., Schonlau, M. Welch, W.  1998\.  Efficient Global Optimization of Expensive Black-Box Functions Efficient global optimization of expensive black-box functions.  Journal of Global Optimization13455-492.
  * Kastner . (2002) Kastner_2002Kastner, J., Solomon, J. Fraser, S.  2002\.  Modeling a Hox Gene Network in Silico using a Stochastic Simulation Algorithm Modeling a hox gene network in silico using a stochastic simulation algorithm.  Developmental Biology246122-131.
  * Kleijnen (2008) Kleijnen_2008Kleijnen, JPC.  2008\.  Design and Analysis of Simulation Experiements Design and analysis of simulation experiements.  New YorkSpringer.
  * Lehmann Casella (1998) Lehmann_Casella_1998Lehmann, E. Casella, G.  1998\.  Theory of Point Estimation Theory of point estimation.  New YorkSpringer-Verlag.
  * Loeppky . (2009) Loeppky_2009Loeppky, JL., Sachs, J. Welch, WJ.  2009\.  Choosing the sample size of a computer experiment: A practical guide Choosing the sample size of a computer experiment: A practical guide.  Technometrics51366-376.
  * Ng Chick (2006) Ng_Chick_2006Ng, SH. Chick, SE.  2006\.  Reducing Parameter Uncertainty for Stochastic Systems Reducing parameter uncertainty for stochastic systems.  ACM Transactions on Modeling and Computer Simulation1626-51.
  * Picheny . (2010) Picheny_2010Picheny, V., Ginsbourger, D., Roustant, O., Haftka, RT. Kim, N.  2010\.  Adaptive Designs of Experiments for Accurate Approximation of a Target Region Adaptive designs of experiments for accurate approximation of a target region.  Journal of Mechanical Design132071008.
  * Santner . (2003) Santer_2003Santner, TJ., Williams, BJ. Notz, WI.  2003\.  The Design and Analysis of Computer Experiments The design and analysis of computer experiments.  New YorkSpringer.
  * Serfling (2002) Serfling_2002Serfling, RJ.  2002\.  Approximation Theorems of Mathematical Statistics Approximation theorems of mathematical statistics.  New YorkWiley.
  * Severini (2005) Severini_2005Severini, T.  2005\.  Elements of Distribution Theory Elements of distribution theory.  New YorkCambridge University Press.
  * Shao Tu (1995) Shao_1995Shao, J. Tu, D.  1995\.  The Jackknife and Bootstrap The jackknife and bootstrap.  New YorkSpringer.
  * Song Nelson (2013) SongNelson_2013Song, E. Nelson, BL.  2013\.  A Quicker Assessment of Input Uncertainty A quicker assessment of input uncertainty.  Proceedings of the 2013 Winter Simulation Conference Proceedings of the 2013 winter simulation conference ( 474-485).  Piscataway, New Jersey: Institute of Electrical and Electronics Engineers, Inc.
  * Sun Farooq (2002) SunFarooq_2002Sun, H. Farooq, M.  2002\.  Note on the generation of random points uniformly distributed in hyper-ellipsoids Note on the generation of random points uniformly distributed in hyper-ellipsoids.  Proceedings of the Fifth International Conference on Information Fusion Proceedings of the fifth international conference on information fusion ( 489-496).
  * Van Der Vaart (1998) Vaart_1998Van Der Vaart, AW.  1998\.  Asymptotic Statistics Asymptotic statistics.  Cambridge, UKCambridge University Press.
  * Wang . (2012) Wang_2012Wang, H., Ciaurri, DE., Durlofsky, LJ. Cominelli, A.  2012\.  Optimal Well Placement under Uncertainty using a Retrospective Optimization Framework Optimal well placement under uncertainty using a retrospective optimization framework.  SPE Journal17112-121.
  * Wieland . (2003) Wieland_Pasupathy_Schmeiser_2003Wieland, JR., Pasupathy, R. Schmeiser, BW.  2003\.  Queueing-Network Stability: Simulation-Based Checking Queueing-network stability: Simulation-based checking.  Proceedings of the 2003 Winter Simulation Conference Proceedings of the 2003 winter simulation conference ( 520-527).  Piscataway, New Jersey: Institute of Electrical and Electronics Engineers, Inc.
  * Xie . (2010) xie_nelson_staum_2010Xie, W., Nelson, BL. Staum, J.  2010\.  The Influence of Correlation Functions on Stochastic Kriging Metamodels The influence of correlation functions on stochastic kriging metamodels.  Proceedings of the 2010 Winter Simulation Conference Proceedings of the 2010 winter simulation conference ( 1067-1078).  Piscataway, New Jersey: Institute of Electrical and Electronics Engineers, Inc.
  * Zouaoui Wilson (2003) Zouaoui_Wilson_2003Zouaoui, F. Wilson, JR.  2003\.  Accounting for Parameter Uncertainty in Simulation Input Modeling Accounting for parameter uncertainty in simulation input modeling.  IIE Transactions35781-792.
  * Zouaoui Wilson (2004) Zouaoui_Wilson_2004Zouaoui, F. Wilson, JR.  2004\.  Accounting for Input-Model and Input-Parameter Uncertainties in Simulation Accounting for input-model and input-parameter uncertainties in simulation.  IIE Transactions361135-1151.

##  8 Appendix (Intended for an Online Companion)

In this appendix we prove Theorems 1–3 and provide a brief description of the
experiment design used to build stochastic kriging metamodels.

To be self-contained, we first state some definitions, lemmas and theorems
that are used in the proofs. Let $\stackrel{{\scriptstyle D}}{{\rightarrow}}$
denote convergence in distribution.

  * •

Borel-Cantelli Lemma (Billingsley, 1995): For events $A_{1},A_{2},\ldots$, if
$\sum_{n=1}^{\infty}\Pr(A_{n})$ converges, then

$\Pr\left(\limsup_{n}A_{n}\right)=0$

where

$\limsup_{n}A_{n}=\cap_{n=1}^{\infty}\cup_{k=n}^{\infty}A_{k}$

is the set of outcomes that occur infinitely many times.

  * •

Lemma 2.11 (Van Der Vaart, 1998): Suppose that
$\mathbf{X}_{n}\stackrel{{\scriptstyle D}}{{\rightarrow}}\mathbf{X}$ for a
random vector X with a continuous distribution function. Then the distribution
function of $\textbf{X}_{n}$ converges uniformly to that of X: $\parallel
F_{\textbf{X}_{n}}-F_{\textbf{X}}\parallel_{\infty}\rightarrow 0$, where
$\parallel h\parallel_{\infty}$ is the sup-norm of $h$ on $\Re$, $\parallel
h\parallel_{\infty}=\sup_{t}|h(t)|$.

  * •

Portmanteau Lemma (Van Der Vaart, 1998): For any random vectors
$\mathbf{X}_{n}$ and $\mathbf{X}$ the following statements are equivalent.

    1. 1.

$\mathbf{X}_{n}\stackrel{{\scriptstyle D}}{{\rightarrow}}\mathbf{X}.$

    2. 2.

$\mbox{E}[f(\mathbf{X}_{n})]\rightarrow\mbox{E}[f(\mathbf{X})]$ for all
bounded, continuous functions $f$.

  * •

Theorem 2.3 (Van Der Vaart, 1998): Let $g:\Re^{k}\rightarrow\Re^{m}$ be
continuous at every point in a set $\mathcal{C}$ such that
$\Pr\\{X\in\mathcal{C}\\}=1$. Then

    1. 1.

If $X_{n}\stackrel{{\scriptstyle D}}{{\longrightarrow}}X$ then
$g(X_{n})\stackrel{{\scriptstyle D}}{{\longrightarrow}}g(X)$.

    2. 2.

If $X_{n}\stackrel{{\scriptstyle P}}{{\longrightarrow}}X$ then
$g(X_{n})\stackrel{{\scriptstyle P}}{{\longrightarrow}}g(X)$.

    3. 3.

If $X_{n}\stackrel{{\scriptstyle a.s.}}{{\longrightarrow}}X$ then
$g(X_{n})\stackrel{{\scriptstyle a.s.}}{{\longrightarrow}}g(X)$.

In the proofs when we refer to the “continuous mapping theorem” we will mean
Theorem 2.3.

  * •

Glivenko-Cantelli Theorem (Van Der Vaart, 1998): If $X_{1},X_{2},\ldots,X_{n}$
are i.i.d. random variables with distribution function $F$ and $F_{n}$ is the
empirical cdf of $X_{1},X_{2},\ldots,X_{n}$, then $\parallel
F_{n}-F\parallel_{\infty}\stackrel{{\scriptstyle a.s.}}{{\longrightarrow}}0$
as $n\rightarrow\infty$.

  * •

Lemma 21.2 (Van Der Vaart, 1998): For cdf $F$, define the inverse cdf to be

$F^{-1}(p)=\inf\\{t{:}\ F(t)\geq p\\}.$

Then a sequence of cdfs $F_{n}(t)\rightarrow F(t)$ for every $t$ where $F$ is
continuous if and only if $F_{n}^{-1}(p)\rightarrow F^{-1}(p)$ for every $p$
where $F^{-1}$ is continuous.

  * •

Theorem 13.1 (Severini, 2005): Let $\mathbf{X}_{1},\mathbf{X}_{2},\ldots$
denote a sequence of $d$-dimensional random vectors such that, for some vector
$\mu$,

$\sqrt{n}(\mathbf{X}_{n}-\mbox{\boldmath$\mu$})\stackrel{{\scriptstyle
D}}{{\rightarrow}}\mbox{N}(\mathbf{0}_{d\times 1},\Sigma)\mbox{ as
}n\rightarrow\infty,$

where $\Sigma$ is a $d\times d$ positive definite matrix with
$|\Sigma|<\infty$. Let $g{:}\ \Re^{d}\rightarrow\Re^{k}$ denote a continuously
differentiable function and let $\nabla g(\mathbf{x})$ denote the $d\times k$
matrix of partial derivatives of $g$ with respect to $\mathbf{x}$. Then

$\sqrt{n}(g(\mathbf{X}_{n})-g(\mbox{\boldmath$\mu$}))\stackrel{{\scriptstyle
D}}{{\rightarrow}}\mbox{N}(\mathbf{0}_{k\times 1},\nabla
g(\mbox{\boldmath$\mu$})^{\top}\Sigma\nabla g(\mbox{\boldmath$\mu$}))\mbox{ as
}n\rightarrow\infty.$

  * •

Theorem 3.8 (Shao Tu, 1995): Let
$\mathbf{X}_{1},\mathbf{X}_{2},\ldots,\mathbf{X}_{m}$ denote $d$-dimensional
i.i.d. random vectors and
$\bar{\mathbf{X}}_{m}=m^{-1}\sum_{i=1}^{m}\mathbf{X}_{i}$. Let
$\bar{\mathbf{X}}_{m}^{*}=m^{-1}\sum_{i=1}\mathbf{X}^{*}_{i}$ where
$\\{\mathbf{X}_{1}^{*},\mathbf{X}_{2}^{*},\ldots,\mathbf{X}_{m}^{*}\\}$ are
randomly and independently drawn with replacement from
$\\{\mathbf{X}_{1},\mathbf{X}_{2},\ldots,\mathbf{X}_{m}\\}$. Let $g{:}\
\Re^{d}\rightarrow\Re^{k}$ denote a continuously differentiable function and
$\nabla g({\bf x})$ denote the $d\times k$ matrix of partial derivatives of
$g$ with respect to ${\bf x}$. Let $T_{m}=g(\bar{\mathbf{X}}_{m})$ and denote
the bootstrap variance estimator for $T_{m}$ by
$v_{m}^{*}=\mbox{Var}_{*}[g(\bar{\mathbf{X}}_{m}^{*})]$.

Suppose that $\mbox{E}[\mathbf{X}_{1}^{\top}\mathbf{X}_{1}]<\infty$ and
$\nabla g(\mbox{\boldmath$\mu$})\neq\mathbf{0}_{d\times k}$ where
$\mbox{\boldmath$\mu$}=\mbox{E}[\mathbf{X}_{1}]$. Suppose further that

$\max_{i_{1},\ldots,i_{m}}|T_{m}(\mathbf{X}_{i_{1}},\ldots,\mathbf{X}_{i_{m}})-T_{m}|/\tau_{m}\stackrel{{\scriptstyle
a.s.}}{{\rightarrow}}0,$ (17)

where the maximum is taken over all integers $i_{1},\ldots,i_{m}$ satisfying
$1\leq i_{1}\leq\cdots\leq i_{m}\leq m$, and $\\{\tau_{m}\\}$ is a sequence of
positive numbers satisfying $\lim\inf_{m}\tau_{m}>0$ and
$\tau_{m}=O(e^{m^{q}})$ with a $q\in(0,1/2).$ Then $v_{m}^{*}$ is strongly
consistent, i.e., $v_{m}^{*}/\sigma_{m}^{2}\stackrel{{\scriptstyle
a.s.}}{{\rightarrow}}1$, where $\sigma_{m}^{2}=m^{-1}\nabla
g(\mbox{\boldmath$\mu$})^{\top}\Sigma\nabla g(\mbox{\boldmath$\mu$})$ and
$\Sigma=\mbox{Var}(\mathbf{X}_{1}).$

  * •

Theorem 1.1 (Lehmann Casella (1998), Chapter 6): Let
$X_{1},X_{2},\ldots,X_{m}$ be i.i.d. with $\mbox{E}(X_{1})=\mu$,
$\mbox{Var}(X_{1})=\sigma^{2}$, and finite fourth moment, and suppose $h$ is a
function of a real variable whose first four derivatives
$h^{\prime}(x),h^{\prime\prime}(x),h^{(3)}(x)$ and $h^{(4)}(x)$ exist for all
$x\in I$, where $I$ is an interval with $\Pr(X_{1}\in I)=1$. Furthermore,
suppose that $|h^{(4)}(x)|\leq M$ for all $x\in I$, for some $M<\infty$. Then

$\mbox{E}[h(\bar{X})]=h(\mu)+\frac{\sigma^{2}}{2m}h^{\prime\prime}(\mu)+\mathcal{R}_{m}.$

If, in addition, the fourth derivative of $h^{2}$ is also bounded, then

$\mbox{Var}[h(\bar{X})]=\frac{\sigma^{2}}{m}[h^{\prime}(\mu)]^{2}+\mathcal{R}_{m}.$

In both cases the remainder $\mathcal{R}_{m}$ is $O(1/m^{2})$.

  * •

Multivariate Taylor Formula (Serfling (2002), page 44): Let the function $g$
defined on $\Re^{d}$ posses continuous partial derivatives of order $n$ at
each point of an open set $S\subset\Re^{d}$. Let ${\bf x}\in S$. For each
point $\mathbf{y}$, $\mathbf{y}\neq{\bf x}$, such that the line segment
$L({\bf x},\mathbf{y})$ joining ${\bf x}$ and $\mathbf{y}$ lies in $S$, there
exists a point $\mathbf{z}$ in the interior of $L(\mathbf{x},\mathbf{y})$ such
that

$\displaystyle g(\mathbf{y})=g({\bf
x})+\left.\sum_{k=1}^{n-1}\frac{1}{k!}\sum_{i_{1}=1}^{d}\cdots\sum_{i_{k}=1}^{d}\frac{\partial^{k}g(t_{1},\ldots,t_{d})}{\partial
t_{i_{1}}\cdots\partial t_{i_{k}}}\right|_{\mathbf{t}={\bf
x}}\cdot\prod_{j=1}^{k}(y_{i_{j}}-x_{i_{j}})$
$\displaystyle+\left.\frac{1}{n!}\sum_{i_{1}=1}^{d}\cdots\sum_{i_{n}=1}^{d}\frac{\partial^{n}g(t_{1},\ldots,t_{d})}{\partial
t_{i_{1}}\cdots\partial
t_{i_{n}}}\right|_{\mathbf{t}=\mathbf{z}}\cdot\prod_{j=1}^{n}(y_{i_{j}}-x_{i_{j}}).$

###  8.1 Asymptotic Consistency of $\mbox{CI}_{+}$

To prove Theorem 1, we first establish three supporting lemmas.

######  Lemma 8.1.

Suppose that Assumptions 1–2 hold. Then the bootstrap resampled moments
converge almost surely to the true moments
$\widehat{\mathbf{X}}_{m}\stackrel{{\scriptstyle
a.s.}}{{\rightarrow}}\mathbf{x}_{c}$ as $m\rightarrow\infty$.

Proof: Since all of the input processes are independent, we establish the
result for one input distribution $F^{c}$ without loss of generality. We prove
the result for $\mathbf{x}_{c}$ being the generic $h$th-order moment,
$\alpha_{h}\equiv\mbox{E}(Z^{h})<\infty$, for $Z\sim F^{c}$.

The $h$th-order bootstrap resampled moment is

$\widehat{X}_{m}=\frac{1}{m}\sum_{j=1}^{m}({\color[rgb]{0,0,0}Z^{(j;m)}})^{h}\mbox{
with }{\color[rgb]{0,0,0}Z^{(j;m)}}\stackrel{{\scriptstyle
i.i.d}}{{\sim}}\mathbf{Z}_{m}^{(0)}$ (18)

where “${\color[rgb]{0,0,0}Z^{(j;m)}}\sim\mathbf{Z}_{m}^{(0)}$” denotes the
$j$th independent sample with replacement from $\mathbf{Z}_{m}^{(0)}$. We use
the Chebychev Inequality and the Borel-Cantelli Lemma to prove the result.

By the Chebychev Inequality, for every $\epsilon>0$, we have

$\mbox{Pr}\left\\{|\widehat{X}_{m}-\alpha_{h}|>\epsilon\right\\}\leq\frac{\mbox{E}\left[(\widehat{X}_{m}-\alpha_{h})^{4}\right]}{\epsilon^{4}}.$
(19)

Notice that

$\mbox{E}\left[(\widehat{X}_{m}-\alpha_{h})^{4}\right]=\mbox{E}\left[\widehat{X}_{m}^{4}\right]-4\alpha_{h}\mbox{E}\left[\widehat{X}_{m}^{3}\right]+6\alpha_{h}^{2}\mbox{E}\left[\widehat{X}_{m}^{2}\right]-4\alpha_{h}^{3}\mbox{E}\left[\widehat{X}_{m}\right]+\alpha_{h}^{4}.$
(20)

We will analyze each term in Equation (20). First, we show that any $i$th
bootstrap resampled moment, denoted as $\widehat{\alpha}_{i}$, is unbiased,

$\displaystyle\mbox{E}\left[\widehat{\alpha}_{i}\right]\equiv\mbox{E}\left[\frac{1}{m}\sum_{j=1}^{m}({\color[rgb]{0,0,0}Z^{(j;m)}})^{i}\right]$
$\displaystyle=$
$\displaystyle\mbox{E}\left[\mbox{E}\left[({\color[rgb]{0,0,0}Z^{(j;m)}})^{i}|Z_{1}^{(0)},\ldots,Z_{m}^{(0)}\right]\right]$
$\displaystyle=$
$\displaystyle\mbox{E}\bigg{[}\frac{1}{m}\sum_{j=1}^{m}(Z_{j}^{(0)})^{i}\bigg{]}$
$\displaystyle=$ $\displaystyle\alpha_{i}.$

Thus, $\mbox{E}[\widehat{X}_{m}]=\alpha_{h}$. Notice that

$\displaystyle\mbox{E}\left[\widehat{X}_{m}^{2}\right]=\mbox{E}\left[\left(\frac{1}{m}\sum_{j=1}^{m}({\color[rgb]{0,0,0}Z^{(j;m)}})^{h}\right)^{2}\right]$
$\displaystyle=$
$\displaystyle\frac{1}{m^{2}}\mbox{E}\left[\left(\sum_{j=1}^{m}({\color[rgb]{0,0,0}Z^{(j;m)}})^{h}\right)^{2}\right]$
$\displaystyle=$
$\displaystyle\frac{1}{m^{2}}\mbox{E}\left[\sum_{j=1}^{m}({\color[rgb]{0,0,0}Z^{(j;m)}})^{2h}+\sum_{i\neq
j}({\color[rgb]{0,0,0}Z^{(i;m)}})^{h}({\color[rgb]{0,0,0}Z^{(j;m)}})^{h}\right]$
$\displaystyle=$
$\displaystyle\frac{1}{m^{2}}\left(m\alpha_{2h}+m(m-1)\mbox{E}\left[\mbox{E}[({\color[rgb]{0,0,0}Z^{(i;m)}})^{h}|Z_{1}^{(0)},\ldots,Z_{m}^{(0)}]\cdot\mbox{E}[({\color[rgb]{0,0,0}Z^{(j;m)}})^{h}|Z_{1}^{(0)},\ldots,Z_{m}^{(0)}]\right]\right)$
$\displaystyle=$
$\displaystyle\frac{1}{m^{2}}\left(m\alpha_{2h}+m(m-1)\mbox{E}\bigg{[}\bigg{(}\frac{1}{m}\sum_{i=1}^{m}(Z_{i}^{(0)})^{h}\bigg{)}^{2}\bigg{]}\right)$
$\displaystyle=$
$\displaystyle\frac{1}{m^{2}}\left(m\alpha_{2h}+\frac{m(m-1)}{m^{2}}\mbox{E}\bigg{[}\sum_{i=1}^{m}(Z_{i}^{(0)})^{2h}+\sum_{i\neq
j}(Z_{i}^{(0)})^{h}(Z_{j}^{(0)})^{h}\bigg{]}\right)$ $\displaystyle=$
$\displaystyle\frac{1}{m^{2}}\left(m\alpha_{2h}+\frac{m(m-1)}{m^{2}}(m\alpha_{2h}+m(m-1)\alpha_{h}^{2})\right)$
$\displaystyle=$
$\displaystyle\frac{1}{m^{2}}[(2m-1)\alpha_{2h}+(m-1)^{2}\alpha_{h}^{2}]$
$\displaystyle=$
$\displaystyle\frac{2}{m}\alpha_{2h}+\left(1-\frac{2}{m}\right)\alpha_{h}^{2}+O(m^{-2})$

where $O(m^{-2})$ means terms at most order $1/m^{2}$. Similar derivations
show that

$\displaystyle\mbox{E}\left[\widehat{X}_{m}^{3}\right]=\frac{1}{m^{4}}\Big{(}[m(4m-3)+(m-1)(m-2)]\alpha_{3h}$
$\displaystyle+\
[3m(m-1)^{2}+3(m-1)^{2}(m-2)]\alpha_{h}\alpha_{2h}+(m-1)^{2}(m-2)^{2}\alpha_{h}^{3}\Big{)}$
$\displaystyle=$
$\displaystyle\frac{6}{m}\alpha_{h}\alpha_{2h}+\left(1-\frac{6}{m}\right)\alpha_{h}^{3}+O(m^{-2})$

and

$\mbox{E}\left[\widehat{X}_{m}^{4}\right]=\frac{12}{m}\alpha_{h}^{2}\alpha_{2h}+\left(1-\frac{12}{m}\right)\alpha_{h}^{4}+O(m^{-2}).$

Thus,

$\displaystyle\mbox{E}\left[(\widehat{X}_{m}-\alpha_{h})^{4}\right]=\mbox{E}[\widehat{X}_{m}^{4}]-4\alpha_{h}\mbox{E}[\widehat{X}_{m}^{3}]+6\alpha_{h}^{2}\mbox{E}[\widehat{X}_{m}^{2}]-4\alpha_{h}^{3}\mbox{E}[\widehat{X}_{m}]+\alpha_{h}^{4}$
$\displaystyle=$
$\displaystyle\frac{12}{m}\alpha_{h}^{2}\alpha_{2h}+\left(1-\frac{12}{m}\right)\alpha_{h}^{4}-4\alpha_{h}\left[\frac{6}{m}\alpha_{h}\alpha_{2h}+\left(1-\frac{6}{m}\right)\alpha_{h}^{3}\right]$
$\displaystyle+\
6\alpha_{h}^{2}\left[\frac{2}{m}\alpha_{2h}+\left(1-\frac{2}{m}\right)\alpha_{h}^{2}\right]-3\alpha_{h}^{4}+O(m^{-2})$
$\displaystyle=$ $\displaystyle 0+O(m^{-2})$

because all of the $O(m^{-1})$ terms cancel. Therefore, combining Equations
(19), (20) and (8.1), we have

$\sum_{m=1}^{\infty}\mbox{Pr}\\{|\widehat{X}_{m}-\alpha_{h}|>\epsilon\\}\leq\sum_{m=1}^{\infty}\frac{c}{m^{2}\epsilon^{4}}<\infty$

where $c$ is some finite constant. Thus, if $\alpha_{4h}<\infty$, then
$\widehat{X}_{m}\stackrel{{\scriptstyle a.s.}}{{\rightarrow}}\alpha_{h}$ by
the first Borel-Cantelli Lemma in Section 4 of Billingsley (1995).

Since Assumption 2 guarantees $m_{\ell}\rightarrow\infty$ for each moment
associated with the $\ell$th input distribution, we can generalize the almost
sure convergence to a vector of moments by applying the converging together
lemma. Therefore, we have
$\widehat{\mathbf{X}}_{\mathbf{m}}\stackrel{{\scriptstyle
a.s.}}{{\rightarrow}}\mathbf{x}_{c}$.

Remark: The independent variables in our stochastic kriging metamodel consist
of central moments and standardized central moments. Since standardized
moments are continuous functions of raw moments, we can use the continuous
mapping theorem to obtain corresponding almost sure convergence of the
standardized moments.

Given a fixed and finite number of design points ${\bf x}_{1},{\bf
x}_{2},\ldots,{\bf x}_{k}$, let $\mathbf{M}=(M({\bf x}_{1}),M({\bf
x}_{2}),\ldots,M({\bf x}_{k}))^{\top}$. The simulation error at design point
${\bf x}_{i}$ is $\epsilon({\bf x}_{i})$, so let $\bar{\epsilon}({\bf
x}_{i})=\sum_{j=1}^{n_{i}}\epsilon({\bf x}_{i})/n_{i}$ for $i=1,2,\ldots,k$
denote the average. Therefore, the sample means of simulation outputs at all
design points can be represented as
$\bar{\mathbf{Y}}_{\mathcal{D}}=\mathbf{M}+\bar{\boldmath{\epsilon}}$, where
$\bar{\boldmath{\epsilon}}=(\bar{\epsilon}({\bf x}_{1}),\bar{\epsilon}({\bf
x}_{2}),\ldots,\bar{\epsilon}({\bf x}_{k}))^{\top}$. Finally, let
$M_{p}(\cdot)$ be a GP having the conditional distribution of $M(\cdot)$ given
$\bar{\mathbf{Y}}_{\mathcal{D}}$.

######  Lemma 8.2.

Suppose Assumptions 3–4 hold. Then $M_{p}(\cdot)$ has continuous sample paths
almost surely.

Proof: Let $(\Omega_{M},P_{M})$ be the underlying probability space for the GP
$M(\cdot)$, and $(\Omega_{\epsilon},P_{\epsilon})$ be the underlying
probability space for $\bar{\boldmath{\epsilon}}$. Notice that
$(\Omega_{\epsilon},P_{\epsilon})$ depends on the particular design points
${\bf x}_{1},{\bf x}_{k},\ldots,{\bf x}_{k}$ and corresponding numbers of
replications $n_{1},n_{2},\ldots,n_{k}$ which we consider fixed and given,
while $(\Omega_{M},P_{M})$ does not.

Let $\omega_{M}\in\Omega_{M}$ be an elementary outcome and
$M(\cdot,\omega_{M})$ the resulting random function. For notational
convenience, let $\mathbf{M}(\omega_{M})=\left(M({\bf
x}_{1},\omega_{M}),M({\bf x}_{2},\omega_{M}),\ldots,M({\bf
x}_{k},\omega_{M})\right)^{\top}$ the random function evaluated at ${\bf
x}_{1},{\bf x}_{2},\ldots,{\bf x}_{k}$. Similarly,
$\bar{\boldmath{\epsilon}}=\bar{\boldmath{\epsilon}}(\omega_{\epsilon})$ for
elementary outcome $\omega_{\epsilon}\in\Omega_{\epsilon}$. Notice that under
Assumption 3, $\bar{\boldmath{\epsilon}}(\omega_{\epsilon})$ has a
multivariate normal distribution.

Theorem 3.4.1 of Adler (2010) asserts that there is a $P_{M}$-measurable set
$\Omega_{M}^{c}\subset\Omega_{M}$ such that
$\Pr\\{\omega_{M}\in\Omega_{M}^{c}\\}=P_{M}(\Omega_{M}^{c})=1$, and for every
$\omega_{M}\in\Omega_{M}^{c}$ the function $M(\cdot,\omega_{M})$ is
continuous.

The random variable $\bar{{\bf Y}}_{\mathcal{D}}$ maps
$\Omega=\Omega_{M}\times\Omega_{\epsilon}\rightarrow\Re^{k}$ as $\bar{{\bf
Y}}_{\mathcal{D}}(\omega)=\mathbf{M}(\omega_{M})+\bar{\mathbf{\epsilon}}(\omega_{\epsilon})$
for $\omega=(\omega_{M},\omega_{\epsilon})\in\Omega$ with probability measure
$P=P_{M}\cdot P_{\epsilon}$ since they are independent. Our goal is to prove
that

$\Pr\\{\omega_{M}\in\bar{\Omega}_{M}^{c}|\bar{{\bf Y}}_{D}\\}=0$ (23)

almost surely.

We know that $0\leq\Pr\\{\omega_{M}\in\bar{\Omega}_{M}^{c}|\bar{{\bf
Y}}_{D}\\}\leq 1$ with probability $1$. But also

$0=\Pr\\{\omega_{M}\in\bar{\Omega}_{M}^{c}\\}={\rm
E}\left[\Pr\\{\omega_{M}\in\bar{\Omega}_{M}^{c}|\bar{{\bf Y}}_{D}\\}\right].$

Therefore, (23) must hold.

######  Lemma 8.3.

Suppose that Assumptions 1–4 hold. Then
$M_{p}(\widehat{\mathbf{X}}_{\mathbf{m}})\stackrel{{\scriptstyle
a.s.}}{{\rightarrow}}M_{p}(\mathbf{x}_{c})$ as $m\rightarrow\infty$.

Proof: Under Assumption 3, the GP $M(\cdot)$ has continuous sample paths
almost surely; applying Lemma 2, $M_{p}(\cdot)$ also has continuous sample
paths almost surely. Under Assumptions 1–2,
$\widehat{\mathbf{X}}_{\mathbf{m}}\stackrel{{\scriptstyle
a.s.}}{{\rightarrow}}\mathbf{x}_{c}$ as $m\rightarrow\infty$ by Lemma 1. And
$M_{p}(\cdot)$ and $\widehat{\mathbf{X}}_{\mathbf{m}}$ are independent. The
result follows by applying the continuous mapping theorem.

######  Theorem 8.1.

Suppose that Assumptions 1–4 hold. Then the interval $[M_{(\lceil
B\frac{\alpha}{2}\rceil)},M_{(\lceil B(1-\frac{\alpha}{2})\rceil)}]$ is
asymptotically consistent, meaning

$\lim_{m\rightarrow\infty}\lim_{B\rightarrow\infty}\Pr\\{M_{(\lceil
B\alpha/2\rceil)}\leq M_{p}(\mathbf{x}_{c})\leq M_{(\lceil
B(1-\alpha/2)\rceil)}\\}=1-\alpha.$ (24)

Proof: Define
$K_{\mathbf{m}}(t)\equiv\Pr\left\\{M_{p}(\widehat{\mathbf{X}}_{\mathbf{m}})\leq
t\right\\}.$ Notice that the distribution $K_{\mathbf{m}}(t)$ depends on both
the distributions of $M_{p}(\cdot)$ and $\widehat{\mathbf{X}}_{\mathbf{m}}$.
Specifically,

$\displaystyle K_{\mathbf{m}}(t)$ $\displaystyle=$
$\displaystyle\int\Pr\left\\{M_{p}({\bf x})\leq t|\widehat{{\bf X}}_{\bf
m}={\bf x}\right\\}d\widehat{F}_{{\bf X}_{\bf m}}({\bf x}|{\bf z}_{\bf
m}^{(0)})$ $\displaystyle=$ $\displaystyle\int\Phi\left(\frac{t-m_{p}({\bf
x})}{\sigma_{p}({\bf x})}\right)d\widehat{F}_{{\bf X}_{\bf m}}({\bf x}|{\bf
z}_{\bf m}^{(0)}).$

Thus, $K_{\mathbf{m}}(t)$ is a continuous distribution almost surely. Let
$\widehat{K}_{\bf m}$ be the empirical cdf of $M_{1},M_{2},\ldots,M_{B}$,
which are i.i.d. from $K_{\mathbf{m}}(t)$. Notice that $\widehat{K}_{\bf
m}^{-1}(\gamma)=M_{(\lceil B\gamma\rceil)}$ for $\gamma=\alpha/2$ and
$1-\alpha/2$.

By the Glivenko-Cantelli Theorem (Van Der Vaart, 1998), $||K_{\bf
m}-\widehat{K}_{\bf m}||_{\infty}\stackrel{{\scriptstyle
a.s.}}{{\longrightarrow}}0$ as $B\rightarrow\infty$. Therefore, by Lemma 21.2
of Van Der Vaart (1998),

$|M_{(\lceil B\gamma\rceil)}-K_{\bf m}^{-1}(\gamma)|\stackrel{{\scriptstyle
a.s.}}{{\longrightarrow}}0$

as $B\rightarrow\infty$ for $\gamma=\alpha/2,1-\alpha/2$. As a result,

$\lim_{B\rightarrow\infty}\Pr\\{M_{(\lceil B\alpha/2\rceil)}\leq
M_{p}(\mathbf{x}_{c})\leq M_{(\lceil B(1-\alpha/2)\rceil)}\\}=\Pr\\{K_{\bf
m}^{-1}(\alpha/2)\leq M_{p}(\mathbf{x}_{c})\leq K_{\bf
m}^{-1}(1-\alpha/2)\\}.$

Therefore, Equation (24) becomes

$\lim_{m\rightarrow\infty}\Pr\\{K_{\mathbf{m}}^{-1}(\alpha/2)\leq
M_{p}(\mathbf{x}_{c})\leq K_{\mathbf{m}}^{-1}(1-\alpha/2)\\}=1-\alpha.$ (25)

To show Equation (25), we only need to show that

$\lim_{m\rightarrow\infty}\Pr\\{K_{\mathbf{m}}^{-1}(\alpha/2)>M_{p}(\mathbf{x}_{c})\\}=\alpha/2$

because the proof of the upper bound is similar.

Since, conditional on $\bar{{\bf Y}}_{\mathcal{D}}$,
$M_{p}(\mathbf{x}_{c})\sim\mbox{N}(m_{p}(\mathbf{x}_{c}),\sigma_{p}^{2}(\mathbf{x}_{c}))$,
the cdf $H(t)\equiv\Pr\\{M_{p}(\mathbf{x}_{c})\leq t\\}$ is continuous. By
Lemma 3 and Lemma 2.11 in Van Der Vaart (1998),

$\displaystyle\sup_{t}|\Pr\\{M_{p}(\widehat{\mathbf{X}}_{\mathbf{m}})\leq
t\\}-\Pr\\{M_{p}(\mathbf{x}_{c})\leq t\\}|$ $\displaystyle=$
$\displaystyle\parallel K_{\mathbf{m}}-H\parallel_{\infty}\rightarrow 0\mbox{
as }m\rightarrow\infty.$

Therefore,

$\displaystyle\Pr\\{K_{\mathbf{m}}^{-1}(\alpha/2)>M_{p}(\mathbf{x}_{c})\\}$
$\displaystyle=$ $\displaystyle\Pr\\{\alpha/2\geq
K_{\mathbf{m}}(M_{p}(\mathbf{x}_{c}))\\}$ $\displaystyle=$
$\displaystyle\Pr\\{\alpha/2\geq
H(M_{p}(\mathbf{x}_{c}))\\}+{\color[rgb]{0,0,0}o(1)}$ $\displaystyle=$
$\displaystyle\Pr\\{M_{p}(\mathbf{x}_{c})\leq
H^{-1}(\alpha/2)\\}+{\color[rgb]{0,0,0}o(1)}$ $\displaystyle=$
$\displaystyle\alpha/2+{\color[rgb]{0,0,0}o(1)}.$

Equation (8.1) is obtained because

$|K_{\mathbf{m}}(M_{p}(\mathbf{x}_{c}))-H(M_{p}(\mathbf{x}_{c}))|\leq\parallel
K_{\mathbf{m}}-H\parallel_{\infty}\rightarrow 0\mbox{ as
$m\rightarrow\infty$}.$

Thus, we have

$\displaystyle\lim_{m\rightarrow\infty}\Pr\\{K_{\mathbf{m}}^{-1}(\alpha/2)\leq
M_{p}(\mathbf{x}_{c})\leq K_{\mathbf{m}}^{-1}(1-\alpha/2)\\}$ $\displaystyle=$
$\displaystyle\lim_{m\rightarrow\infty}\Pr\\{M_{p}(\mathbf{x}_{c})\leq
K_{\mathbf{m}}^{-1}(1-\alpha/2)\\}-\lim_{m\rightarrow\infty}\Pr\\{M_{p}(\mathbf{x}_{c})<K_{\mathbf{m}}^{-1}(\alpha/2)\\}$
$\displaystyle=$ $\displaystyle(1-\alpha/2)-\alpha/2$ $\displaystyle=$
$\displaystyle 1-\alpha.$

###  8.2 Asymptotic Analysis of Variance Component Estimators

######  Theorem 8.2.

Suppose that Assumptions 1–4 hold. Then the variance component estimators
$\widehat{\sigma}^{2}_{M},\widehat{\sigma}^{2}_{I},\widehat{\sigma}^{2}_{T}$
are consistent as $m,B\rightarrow\infty$.

Proof: When a GP $M(\cdot)$ has a continuous correlation function with all
parameters finite, the SK predictor

$m_{p}({\bf x})=\widehat{\beta}_{0}+\tau^{2}R({\bf
x})^{\top}[\Sigma+C]^{-1}(\bar{{\bf Y}}_{\mathcal{D}}-\widehat{\beta}_{0}\cdot
1_{k\times 1}),$ (27)

and corresponding variance

$\sigma_{p}^{2}({\bf x})=\tau^{2}-\tau^{4}R({\bf
x})^{\top}[\Sigma+C]^{-1}R({\bf x})+\mathbf{\eta}^{\top}[1_{k\times
1}^{\top}(\Sigma+C)^{-1}1_{k\times 1}]^{-1}\mathbf{\eta}$

where $R({\bf x})^{\top}=(r({\bf x}-{\bf x}_{1}),r({\bf x}-{\bf
x}_{2}),\ldots,r({\bf x}-{\bf x}_{k}))$ and $\eta=1-1_{k\times
1}^{\top}(\Sigma+C)^{-1}\tau^{2}R({\bf x})$, are continuous and bounded
functions of ${\bf x}$.

By the Strong Law of Large Numbers, the raw moment estimator
$\mathbf{X}_{\mathbf{m}}\stackrel{{\scriptstyle a.s.}}{{\rightarrow}}{\bf
x}_{c}$ as $m\rightarrow\infty$ under Assumptions 1–2. This almost sure
convergence can be extended to central moments and standardized central
moments by the continuous mapping theorem. By applying the Portmanteau Lemma
in Van Der Vaart (1998), we have

$\lim_{m\rightarrow\infty}\sigma_{M}^{2}=\lim_{m\rightarrow\infty}\int\sigma^{2}_{p}(\mathbf{x})\,dF^{c}_{\mathbf{X}_{\mathbf{m}}}(\mathbf{x})=\lim_{m\rightarrow\infty}\mbox{E}[\sigma^{2}_{p}(\mathbf{X}_{\mathbf{m}})]=\sigma_{p}^{2}({\bf
x}_{c}),$

and

$\displaystyle\lim_{m\rightarrow\infty}\sigma_{I}^{2}=\lim_{m\rightarrow\infty}\int(m_{p}(\mathbf{x})-\mu_{0})^{2}\,dF^{c}_{\mathbf{X}_{\mathbf{m}}}(\mathbf{x})$
$\displaystyle=$
$\displaystyle\lim_{m\rightarrow\infty}\mbox{E}\left[\big{(}m_{p}(\mathbf{X}_{\mathbf{m}})-\mbox{E}[m_{p}(\mathbf{X}_{\mathbf{m}})]\big{)}^{2}\right]$
$\displaystyle=$ $\displaystyle\big{(}m_{p}({\bf x}_{c})-m_{p}({\bf
x}_{c})\big{)}^{2}=0$

where
$\mu_{0}=\int\int\nu\,dF(\nu|\mathbf{x})\,dF^{c}_{\mathbf{X}_{\mathbf{m}}}(\mathbf{x})=\int
m_{p}(\mathbf{x})\,dF^{c}_{\mathbf{X}_{\mathbf{m}}}(\mathbf{x})$.

Recall that $F(\nu|{\bf x})$ is a normal distribution $\mbox{N}(m_{p}({\bf
x}),\sigma_{p}^{2}({\bf x}))$. Let $g({\bf
x})\equiv\int(\nu-\mu_{0})^{2}\,dF(\nu|{\bf x})$. Then

$\displaystyle\lim_{m\rightarrow\infty}\sigma_{T}^{2}=\lim_{m\rightarrow\infty}\int\int(\nu-\mu_{0})^{2}\,dF(\nu|\mathbf{x})\,dF^{c}_{\mathbf{X}_{\mathbf{m}}}(\mathbf{x})$
$\displaystyle=$ $\displaystyle\lim_{m\rightarrow\infty}\int g({\bf
x})\,dF^{c}_{\mathbf{X}_{\mathbf{m}}}(\mathbf{x}).$

However,

$\displaystyle g({\bf x})$ $\displaystyle=$ $\displaystyle\int(\nu-m_{p}({\bf
x})+m_{p}({\bf x})-\mu_{0})^{2}\,dF(\nu|\mathbf{x})$ $\displaystyle=$
$\displaystyle\int(\nu-m_{p}({\bf x}))^{2}\,dF(\nu|\mathbf{x})+(m_{p}({\bf
x})-\mu_{0})^{2}+(m_{p}({\bf x})-\mu_{0})\int(\nu-m_{p}({\bf
x}))\,dF(\nu|\mathbf{x})$ $\displaystyle=$ $\displaystyle\sigma_{p}^{2}({\bf
x})+(m_{p}({\bf x})-\mu_{0})^{2}+0.$

Since $m_{p}({\bf x})$ and $\sigma_{p}^{2}({\bf x})$ are continuous and
bounded functions, so is $g({\bf x})$. Therefore,

$\lim_{m\rightarrow\infty}\sigma_{T}^{2}=\lim_{m\rightarrow\infty}\mbox{E}[g(\mathbf{X}_{\mathbf{m}})]=g({\bf
x}_{c})=\sigma^{2}_{p}({\bf x}_{c})$

by applying the Portmanteau Lemma.

Next, we will show consistency of the variance estimators. By Lemma 1,
$\widehat{\mathbf{X}}_{\mathbf{m}}\stackrel{{\scriptstyle
a.s.}}{{\rightarrow}}{\bf x}_{c}.$ For the metamodel uncertainty estimator, we
have

$\displaystyle\lim_{m\rightarrow\infty}\lim_{B\rightarrow\infty}\widehat{\sigma}^{2}_{M}=\lim_{m\rightarrow\infty}\lim_{B\rightarrow\infty}\frac{1}{B}\sum_{b=1}^{B}\sigma^{2}_{p}(\widehat{\mathbf{X}}_{\mathbf{m}}^{(b)})$
$\displaystyle=$
$\displaystyle\lim_{m\rightarrow\infty}\mbox{E}\Big{[}\sigma_{p}^{2}(\widehat{\mathbf{X}}_{\mathbf{m}})|\mathbf{Z}_{\mathbf{m}}^{(0)}\Big{]}$
$\displaystyle=$ $\displaystyle\sigma_{p}^{2}({\bf x}_{c}).$

The last step follows by applying the Portmanteau Lemma.

For the input uncertainty estimator, we have

$\displaystyle\lim_{m\rightarrow\infty}\lim_{B\rightarrow\infty}\widehat{\sigma}^{2}_{I}=\lim_{m\rightarrow\infty}\lim_{B\rightarrow\infty}\frac{B}{B-1}\left[\frac{1}{B}\sum_{b=1}^{B}m_{p}^{2}(\widehat{\mathbf{X}}_{\mathbf{m}}^{(b)})-\bar{{\mu}}^{2}\right]$
$\displaystyle=$
$\displaystyle\lim_{m\rightarrow\infty}\bigg{(}\mbox{E}\Big{[}m_{p}^{2}(\widehat{\mathbf{X}}_{\mathbf{m}})|\mathbf{Z}_{\mathbf{m}}^{(0)}\Big{]}-\mbox{E}^{2}\Big{[}m_{p}(\widehat{\mathbf{X}}_{\mathbf{m}})|\mathbf{Z}_{\mathbf{m}}^{(0)}\Big{]}\bigg{)}$
$\displaystyle=$ $\displaystyle m_{p}^{2}({\bf x}_{c})-m_{p}^{2}({\bf
x}_{c})=0.$

The last step follows by applying Lemma 1 and the Portmanteau Lemma.

For the total variance estimator, we have

$\displaystyle\lim_{m\rightarrow\infty}\lim_{B\rightarrow\infty}\widehat{\sigma}^{2}_{T}=\lim_{m\rightarrow\infty}\lim_{B\rightarrow\infty}\frac{B}{B-1}\left(\frac{1}{B}\sum_{b=1}^{B}M_{b}^{2}-\bar{M}^{2}\right)$
$\displaystyle=$
$\displaystyle\lim_{m\rightarrow\infty}\mbox{E}\big{[}M_{p}^{2}(\widehat{\mathbf{X}}_{\mathbf{m}})|\mathbf{Z}_{\mathbf{m}}^{(0)}\big{]}-\lim_{m\rightarrow\infty}\mbox{E}^{2}\big{[}M_{p}(\widehat{\mathbf{X}}_{\mathbf{m}})|\mathbf{Z}_{\mathbf{m}}^{(0)}\big{]}$
$\displaystyle=$ $\displaystyle\mbox{E}[M_{p}^{2}({\bf
x}_{c})]-\mbox{E}^{2}[M_{p}({\bf x}_{c})]$ $\displaystyle=$ $\displaystyle
m_{p}^{2}({\bf x}_{c})+\sigma_{p}^{2}({\bf x}_{c})-m_{p}^{2}({\bf x}_{c})$
$\displaystyle=$ $\displaystyle\sigma^{2}_{p}({\bf x}_{c}).$

By Lemma 3, $M_{p}(\widehat{\mathbf{X}}_{\mathbf{m}})\stackrel{{\scriptstyle
a.s.}}{{\rightarrow}}M_{p}({\bf x}_{c})\mbox{ as }m\rightarrow\infty$. Then
Step (8.2) follows by applying Portmanteau Lemma.

######  Theorem 8.3.

Suppose that Assumptions 1–4 and the following additional assumptions hold:

  1. 5.

The first three derivatives of the correlation function of the GP $M({\bf x})$
exist and the third derivative is bounded; and

  2. 6.

$m_{\ell}/m\rightarrow 1$ for $\ell=1,2,\ldots,L$.

Then
$\lim_{m\rightarrow\infty}m\sigma_{I}^{2}=\lim_{m\rightarrow\infty}\lim_{B\rightarrow\infty}m\widehat{\sigma}^{2}_{I}=\sigma^{2}_{\mu}$
almost surely, where $\sigma^{2}_{\mu}$ is a positive constant.

Proof: Under Assumptions 1–2, and applying the multivariate central limit
theorem, we have as $m\rightarrow\infty$,

$\sqrt{m}(\mathbf{X}_{m}-\mathbf{x}_{c})\stackrel{{\scriptstyle
D}}{{\rightarrow}}\mbox{N}(\mathbf{0}_{d\times 1},\Lambda)$

where $\Lambda$ denotes the $d\times d$ positive definite asymptotic variance-
covariance matrix of ${\bf X}_{\mathbf{m}}$.

When a GP $M({\bf x})$ has a continuous correlation function with all
parameters finite, the SK predictor

$m_{p}({\bf x})=\widehat{\beta}_{0}+\tau^{2}R({\bf
x})^{\top}[\Sigma+C]^{-1}(\bar{{\bf Y}}_{\mathcal{D}}-\widehat{\beta}_{0}\cdot
1_{k\times 1}),$ (29)

given the simulation sample mean $\bar{{\bf Y}}_{\mathcal{D}}$, is continuous
and bounded. Under Assumption 5, the gradient $\nabla m_{p}({\bf x})$ exists
and is continuous. We will show that $\nabla m_{p}({\bf
x})\neq\mathbf{0}_{d\times 1}$ almost surely. By taking the derivative of
$m_{p}({\bf x})$ in Equation (29), we have

$\frac{\partial m_{p}({\bf x})}{\partial
x_{j}}=\underbrace{\frac{\partial{R}({\bf x})^{\top}}{\partial
x_{j}}\tau^{2}[\Sigma+C]^{-1}}_{\mathbf{A}}(\bar{{\bf
Y}}_{\mathcal{D}}-\widehat{\beta}_{0}\cdot 1_{k\times 1}).$ (30)

Since $\partial{R}({\bf x})^{\top}/\partial x_{j}=(\partial{r}({\bf x}-{\bf
x}_{1})/\partial x_{j},\partial{r}({\bf x}-{\bf x}_{2})/\partial
x_{j},\ldots,\partial{r}({\bf x}-{\bf x}_{k})/\partial
x_{j})\neq\mathbf{0}_{1\times k}$ and $\tau^{2}[\Sigma+C]^{-1}$ is positive
definite, $\mathbf{A}$ is a non-zero constant vector. Under Assumption 3,
$\mathbf{A}(\bar{{\bf Y}}_{\mathcal{D}}-\widehat{\beta}_{0}\cdot 1_{k\times
1})$ is a normal random variable that is equal to 0 with probability 0. Thus,
$\nabla m_{p}({\bf x})\neq\mathbf{0}_{d\times 1}$ almost surely. Applying
Theorem 13.1 in Severini (2005), we have

$\sqrt{m}(m_{p}(\mathbf{X}_{{m}})-m_{p}(\mathbf{x}_{c}))\stackrel{{\scriptstyle
D}}{{\rightarrow}}\mbox{N}(0,\sigma^{2}_{\mu})$

where $\sigma_{\mu}^{2}=\nabla m_{p}({\bf x}_{c})^{\top}\Lambda\nabla
m_{p}({\bf x}_{c})>0$. This establishes the constant.

Since $m_{p}(\cdot)$ is continuous and bounded, there always exists a finite
$M_{1}>0$ such that $|m_{p}({\bf x})|<M_{1}$ for all ${\bf x}\in\Re^{d}$.
Therefore, $\max_{{\bf x}\in\Re^{d}}|m_{p}({\bf x})-m_{p}({\bf
X}_{\mathbf{m}})|<2M_{1}$. Let $\tau_{m}=\mbox{e}^{m^{1/4}}$. Since
$2M_{1}/\tau_{m}\rightarrow 0$ as $m\rightarrow\infty$, Condition (17) of
Theorem 3.8 of Shao Tu (1995) holds. Thus, the bootstrap variance estimator is
strongly consistent:
$\lim_{m\rightarrow\infty}\lim_{B\rightarrow\infty}m\widehat{\sigma}_{I}^{2}=\sigma_{\mu}^{2}$
almost surely.

Next, we will show $\lim_{m\rightarrow\infty}m\sigma_{I}^{2}=\sigma_{\mu}^{2}$
by proving a multi-variate version of Theorem 1.1 in Lehmann Casella (1998),
Chapter 6. Let $L({\bf X}_{\mathbf{m}},{\bf x}_{c})$ denote the line segment
joining ${\bf X}_{\mathbf{m}}$ and ${\bf x}_{c}$. By the Multivariate Taylor
Formula (Serfling, 2002),

$m_{p}({\bf X}_{\mathbf{m}})=m_{p}({\bf x}_{c})+\nabla m_{p}({\bf
x}_{c})^{\top}({\bf X}_{\mathbf{m}}-{\bf x}_{c})+\frac{1}{2}({\bf
X}_{\mathbf{m}}-{\bf x}_{c})^{\top}\nabla^{2}m_{p}({\bf x}_{c})({\bf
X}_{\mathbf{m}}-{\bf x}_{c})+\mathcal{R}({\bf X}_{\mathbf{m}},{\bf x}_{c}).$

The remainder term

$\mathcal{R}({\bf X}_{\mathbf{m}},{\bf
x}_{c})=\frac{1}{3!}\sum_{i_{1}=1}^{d}\sum_{i_{2}=1}^{d}\sum_{i_{3}=1}^{d}\left.\frac{\partial^{3}m_{p}(x_{1},\ldots,x_{d})}{\partial
x_{i_{1}}\partial x_{i_{2}}\partial x_{i_{3}}}\right|_{{\bf x}={\bf
z}}\prod_{j=1}^{3}(X_{\mathbf{m},i_{j}}-x_{c,i_{j}})$

where ${\bf z}$ denotes a value in the interior of $L({\bf
X}_{\mathbf{m}},{\bf x}_{c})$, and $X_{\mathbf{m},i}$ and $x_{c,i}$ denote the
$i$th components of the vectors ${\bf X}_{\mathbf{m}}$ and ${\bf x}_{c}$. By
taking the expectation over both sides, we have

$\mbox{E}[m_{p}({\bf X}_{\mathbf{m}})]=m_{p}({\bf
x}_{c})+\frac{1}{2}\mbox{E}\left[({\bf X}_{\mathbf{m}}-{\bf
x}_{c})^{\top}\nabla^{2}m_{p}({\bf x}_{c})({\bf X}_{\mathbf{m}}-{\bf
x}_{c})\right]+\mbox{E}[\mathcal{R}(\mathbf{X}_{\mathbf{m}},{\bf x}_{c})]$
(31)

where $\nabla^{2}$ is the Hessian operator.

We will show that the second and third terms on the RHS of Equation (31) are
$O(m^{-1})$ and $O(m^{-2})$, respectively, under Assumption 5. Since all of
the input processes are independent, we establish the result for one input
distribution $F^{c}$ without loss of generality.

We prove the result for ${\bf x}_{c}$ being the generic $h$th-order moment,
$x_{c,h}=\mbox{E}(Z_{1}^{h})<\infty$, and
$X_{m,h}=m^{-1}\sum_{j=1}^{m}Z_{j}^{h}$ for
$Z_{j}\stackrel{{\scriptstyle\mbox{iid}}}{{\sim}}F^{c}$.

Let $C_{ij}\equiv\frac{1}{2}[\nabla^{2}m_{p}({\bf x}_{c})]_{i,j}$. We first
consider components of the second term on the RHS of Equation (31).

$\displaystyle\mbox{E}[C_{ij}(X_{m,i}-x_{c,i})(X_{m,j}-x_{c,j})]=C_{ij}\mbox{E}\left[\frac{1}{m}\sum_{k_{1}=1}^{m}(Z_{k_{1}}^{i}-x_{c,i})\cdot\frac{1}{m}\sum_{k_{2}=1}^{m}(Z_{k_{2}}^{j}-x_{c,j})\right]$
$\displaystyle=$
$\displaystyle\frac{C_{ij}}{m^{2}}\mbox{E}\left[\sum_{k=1}^{m}(Z_{k}^{i}-x_{c,i})(Z_{k}^{j}-x_{c,j})+\sum_{k_{1}\neq
k_{2}}(Z_{k_{1}}^{i}-x_{c,i})(Z_{k_{2}}^{j}-x_{c,j})\right]$ $\displaystyle=$
$\displaystyle\frac{C_{ij}}{m^{2}}\mbox{E}\left[\sum_{k=1}^{m}(Z_{k}^{i}-x_{c,i})(Z_{k}^{j}-x_{c,j})+0\right]=O\left(\frac{1}{m}\right).$

The last two steps follow because the $Z_{i}$ are i.i.d. and Assumption 5
holds. Thus, the second term on the RHS of Equation (31) is

$\frac{1}{2}\mbox{E}[({\bf X}_{\mathbf{m}}-{\bf
x}_{c})^{\top}\nabla^{2}m_{p}({\bf x}_{c})({\bf X}_{\mathbf{m}}-{\bf
x}_{c})]=\sum_{i=1}^{d}\sum_{j=1}^{d}\mbox{E}[C_{ij}(X_{m,i}-x_{c,i})(X_{m,j}-x_{c,j})]=O\left(\frac{1}{m}\right).$

Similarly, for the components of the third term of the RHS of Equation (31),
we have

$\displaystyle
D_{ijk}\mbox{E}[(X_{m,i}-x_{c,i})(X_{m,j}-x_{c,j})(X_{m,k}-x_{c,k})]$
$\displaystyle=$ $\displaystyle
D_{ijk}\mbox{E}\left[\frac{1}{m^{3}}\sum_{k_{1}=1}^{m}(Z_{k_{1}}^{i}-x_{c,i})\cdot\sum_{k_{2}=1}^{m}(Z_{k_{2}}^{j}-x_{c,j})\cdot\sum_{k_{3}=1}^{m}(Z_{k_{3}}^{k}-x_{c,k})\right]$
$\displaystyle=$
$\displaystyle\frac{D_{ijk}}{m^{3}}\mbox{E}\left[\sum_{k_{1}=1}^{m}(Z_{k_{1}}^{i}-x_{c,i})(Z_{k_{1}}^{j}-x_{c,j})(Z_{k_{1}}^{k}-x_{c,j})+0\right]=O\left(\frac{1}{m^{2}}\right).$

where

$D_{ijk}\equiv\left.\frac{1}{3!}\frac{\partial^{3}m_{p}(x_{1},\ldots,x_{d})}{\partial
x_{i}\partial x_{j}\partial x_{k}}\right|_{{\bf x}={\bf z}}.$

Again, the last two steps follow because the $Z_{i}$ are i.i.d. and Assumption
5 holds. Thus, the third term in Equation (31) is

$\mbox{E}[\mathcal{R}(\mathbf{X}_{\mathbf{m}},{\bf
x}_{c})]=\sum_{i=1}^{d}\sum_{j=1}^{d}\sum_{k=1}^{d}D_{ijk}\mbox{E}[(X_{m,i}-x_{c,i})(X_{m,j}-x_{c,j})(X_{m,k}-x_{c,k})]=O\left(\frac{1}{m^{2}}\right).$

Squaring both sides of Equation (31), we have

$[\mbox{E}(m_{p}({\bf X}_{\mathbf{m}}))]^{2}=m_{p}^{2}({\bf x}_{c})+m_{p}({\bf
x}_{c})\mbox{E}[({\bf X}_{\mathbf{m}}-{\bf x}_{c})^{\top}\nabla^{2}m_{p}({\bf
x}_{c})({\bf X}_{\mathbf{m}}-{\bf x}_{c})]+O\left(\frac{1}{m^{2}}\right).$
(32)

By repeating the same derivation that results in Equation (31) but using
$m_{p}^{2}(\cdot)$ instead of $m_{p}(\cdot)$, we obtain

$\displaystyle\mbox{E}[m_{p}^{2}({\bf X}_{\mathbf{m}})]=m_{p}^{2}({\bf
x}_{c})+\frac{1}{2}\mbox{E}[({\bf X}_{\mathbf{m}}-{\bf
x}_{c})^{\top}\nabla^{2}m_{p}^{2}({\bf x}_{c})({\bf X}_{\mathbf{m}}-{\bf
x}_{c})]+O\left(\frac{1}{m^{2}}\right)$ (33) $\displaystyle=$ $\displaystyle
m_{p}^{2}({\bf x}_{c})+\mbox{E}\Big{[}({\bf X}_{\mathbf{m}}-{\bf
x}_{c})^{\top}\nabla m_{p}({\bf x}_{c})\nabla m_{p}({\bf x}_{c})^{\top}({\bf
X}_{\mathbf{m}}-{\bf x}_{c})$ $\displaystyle+({\bf X}_{\mathbf{m}}-{\bf
x}_{c})^{\top}m_{p}({\bf x}_{c})\nabla^{2}m_{p}({\bf x}_{c})({\bf
X}_{\mathbf{m}}-{\bf x}_{c})\Big{]}+O\left(\frac{1}{m^{2}}\right).$

Then,

$\displaystyle\mbox{Var}[m_{p}({\bf X}_{\mathbf{m}})]=\mbox{E}[m_{p}^{2}({\bf
X}_{\mathbf{m}})]-\Big{(}\mbox{E}[m_{p}({\bf X}_{\mathbf{m}})]\Big{)}^{2}$
$\displaystyle=$ $\displaystyle\mbox{E}[({\bf X}_{\mathbf{m}}-{\bf
x}_{c})^{\top}\nabla m_{p}({\bf x}_{c})\nabla m_{p}({\bf x}_{c})^{\top}({\bf
X}_{\mathbf{m}}-{\bf x}_{c})]+O\left(\frac{1}{m^{2}}\right)$ $\displaystyle=$
$\displaystyle\mbox{E}[\nabla m_{p}({\bf x}_{c})^{\top}({\bf
X}_{\mathbf{m}}-{\bf x}_{c})({\bf X}_{\mathbf{m}}-{\bf x}_{c})^{\top}\nabla
m_{p}({\bf x}_{c})]+O\left(\frac{1}{m^{2}}\right)$ $\displaystyle=$
$\displaystyle\frac{1}{m}\nabla m_{p}({\bf x}_{c})^{\top}\Lambda\nabla
m_{p}({\bf x}_{c})+O\left(\frac{1}{m^{2}}\right).$

Step (8.2) follows because $\nabla m_{p}({\bf x}_{c})^{\top}({\bf
X}_{\mathbf{m}}-{\bf x}_{c})$ is a scalar. Thus, we have
$\lim_{m\rightarrow\infty}m\sigma_{I}^{2}=\sigma_{\mu}^{2}$.

Remark: The independent variables in our stochastic kriging metamodel consist
of central moments and standardized central moments, rather than raw moments.
However, Theorem 3 can easily be extended to central and standardized central
moments as follows.

Since standardized moments are continuous functions of raw moments, denoted
generically as $g(\cdot)$, we can consider the composite function $(m_{p}\circ
g)(\cdot)$ and follow steps analogous to those in the proof of Theorem 3. Up
to the third derivatives we have

$\displaystyle(m_{p}\circ g)^{\prime}(t)$ $\displaystyle=$ $\displaystyle
m_{p}^{\prime}(g(t))g^{\prime}(t)$ $\displaystyle(m_{p}\circ
g)^{\prime\prime}(t)$ $\displaystyle=$ $\displaystyle
m_{p}^{\prime\prime}(g(t))[g^{\prime}(t)]^{2}+m_{p}^{\prime}(g(t))g^{\prime\prime}(t)$
$\displaystyle(m_{p}\circ g)^{(3)}(t)$ $\displaystyle=$ $\displaystyle