Proceedings of the 2018 Winter Simulation Conference
M. Rabe, A. A. Juan, N. Mustafee, A. Skoogh, S. Jain, and B. Johansson, eds.

REVISITING DIRECT BOOTSTRAP RESAMPLING FOR INPUT MODEL UNCERTAINTY
Russell R. Barton

Henry Lam

Department of SCISThe Pennsylvania State University
210 Business Building
University Park, PA 16803, USA

Department of IEOR
Columbia University
500 W. 120th Street
New York, NY 10027, USA

Eunhye SongDepartment of IME
The Pennsylvania State University
310 Leonhard Building
University Park, PA 16802, USA

ABSTRACTMetamodel-based bootstrap methods for characterizing input model uncertainty have disadvantages forsettings where there are a large number of input distributions, or when using empirical distributions to drive
the simulation. Early direct bootstrapping of empirical distributions did not take into account the distinction
between intrinsic and extrinsic variations in the resampled quantities. When the intrinsic uncertainty islarge, the result is overcoverage of the bootstrap percentile intervals. We explore ways of accounting
for both sources in direct bootstrap characterization of input model uncertainty, and study the impact on
confidence interval (CI) coverage. Four new bootstrap-based CIs for the expected simulation output under
the unknown true distribution are proposed, basic shrinkage CI, percentile shrinkage CI, basic hierarchicalbootstrap CI, and percentile hierarchical bootstrap CI, and their empirical performances are demonstratedusing an example.
1

INTRODUCTION

When input models are fitted to data, the finiteness of the data introduces sampling error to the fitted
input distributions. Hence, any simulation analysis that assumes the fitted input distributions are the “true”
real-world distribution representing randomness in data fails to capture the uncertainty in the simulation
output caused by such sampling error. In particular, confidence interval coverage can be severely affectedby input uncertainty. Barton and Schruben (2001) show the coverage of nominal 90% confidence intervals
for the expected waiting time of capacitated queues to be no better than 20% when the sample size ofthe real-world data is 500 and the sampling error of the fitted input distributions is ignored. Input model
uncertainty analysis characterizes the impact of such sampling error on simulation output performancemeasure. For instance, one can compute the confidence interval (CI) for the performance measure that
accounts for the sampling error in the fitted input distributions. Clearly, the sampling distribution of thefitted distributions is unknown, however, we may approximate the sampling distribution by bootstrapping
hinging on the idea that bootstrapped samples resemble the statistical characteristics of real-world samples.
Since first applied by Barton and Schruben (1993), bootstrapping has been a popular tool for quantifyingBarton, Lam, and Song
input model uncertainty (Cheng and Holland 1997; Barton and Schruben 2001; Barton et al. 2014; Songand Nelson 2015).
Input model uncertainty analysis via bootstrapping involves three steps: resampling the input data,
computing the output (simulation run or runs) and using the accumulated outputs to quantify input model
uncertainty by computing relevant statistics. Methods fall into three categories based on how the first two
steps are completed: i) direct resampling of the empirical distributions, followed by simulation runs; ii)
direct parametric bootstrap resampling, assuming that the input distribution parametric families are known,
followed by simulation runs, and iii) parametric bootstrap resampling, evaluated via a previously fitted
simulation metamodel of the response surface rather than via simulation runs. Bayesian approaches can beapplied in settings ii) and iii), where the resampling is replaced by sampling from the posterior distribution
given the sample data.
Approaches i) and ii) may exhibit overcoverage due to improper accounting for the “intrinsic error”,
namely the statistical error arising from the finiteness of the simulation effort in estimating the performancemeasure (Barton 2007). Metamodels used in the approach iii) greatly reduce this phenomenon, and when
stochastic kriging is used with the metamodeling strategy, prediction error can be captured (Barton et al.
2014).
The metamodeling approach iii) has drawbacks as well. When the number of input model parametersis in the hundreds to thousands, fitting a metamodel for this number of variables becomes impractical for
two reasons. First, the computational complexity of estimating the parameters of the metamodel increases
as a function of the number of variables. Second, the number of simulation runs to fit such a model canbe an order of magnitude greater than the number of variables, even in a favorable case. This means thatthe third approach may require more simulation runs than the direct bootstrapping used in the first two
approaches. This problem is compounded if the metamodel must use as input an empirical distributions
rather than parametric ones. For each input distribution, rather than one or two parameters, the metamodel
must take as input the empirical cumulative distribution function (ecdf) values at tens or hundreds of points
corresponding to the original data.
In this paper, we focus on a nonparametric framework for quantifying input uncertainty. This is motivatedfrom the fact that the assumption of known parametric distribution families for input models introduces
its own source of error, whose effect on simulation output can sometimes be hard to quantify. Empirical
distributions for input models avoid this problem, at least in the asymptotic sense. This nonparametricframework necessitates the use of approach i). However, existing methods in i) encounter the problem ofovercoverage when the intrinsic error is relatively large. To account for this intrinsic error, naive remedies
would call for more demanding simulation runs, essentially to wash out the Monte Carlo noise, or else
suffer from statistical inefficiency, as we will describe in the sequel.
Our main contribution is to introduce two new approaches to produce direct bootstrap CIs that accountfor both intrinsic and extrinsic errors in terms of correct coverage, and also with satisfactory statistical
efficiency using affordable computational effort. To attain the latter properties, both of our introduced
approaches hinge on properly adjusting the bootstrap samples from few simulation runs to exhibit statisticalbehaviors as if there were more runs. The first approach uses a shrinkage on the bootstrap samples thatdeflate their variances, while the second approach uses a hierarchical bootstrap that de-biases the quantile
estimates used to construct the CI.
The remainder of this paper is as follows. Section 2 presents our mathematical formulation and thestandard approaches to conduct direct bootstrapping. Sections 3 and 4 explain the problems with theseapproaches and present remedies using shrinkage and the hierarchical bootstrap respectively. Section 5
provides a limited computational comparison. Section 6 concludes with an assessment of implications andfuture work.

Barton, Lam, and Song
2

MATHEMATICAL FRAMEWORK

We follow the notation in Song and Nelson (2015). Suppose the objective of simulation analysis is
to estimate a performance measure of the system, which can be represented as an expected value of asimulation output. The input distribution (cdf) that drives the simulation is denoted as F, which we allow
to be multivariate. The simulation output from the rth replication of the simulation is
Yr (F) = η(F) + εr (F), r = 1, 2, . . . , R,
where η(F) , E[Yr (F)|F] and ε is a random variable representing the between-replication output differences
due to the finiteness of the simulation runs with zero mean and finite variance σ 2 (F). Input model uncertainty
is introduced when we use a fitted distribution Fb in place of the true distribution F c assuming such ab When Fb isdistribution exists. The input distribution fitted to a size-n real-world sample is denoted by F.
used to run simulations,b = η(F)b + εr (F)b = η(F c ) + η(F)b − η(F c ) + εr (F).
b
Yr (F)(1)
Note that Fb is a realization of a random process but F c is deterministic. In that context, (1) shows twob − η(F c ) has been called “bias error” by Cheng
sources of error in the simulation output. The first error η(F)and Holland (2004), “extrinsic error” by Barton et al. (2014), and Song et al. (2014) define its variance as
b comes from the finiteness of the simulation effort, which is“input uncertainty.” The second error, εr (F),referred to as “variance error” by Cheng (1994), “intrinsic error” by Barton et al. (2014), and “simulationb ∼ N(0, σ 2 (F))
b but Barton
error” by Song et al. (2014). It is reasonable in many cases to assume that εr (F)and Schruben (2001) show that extrinsic error can be quite non-Gaussian.
2.1 Basic and Percentile Bootstraps
Our objective is to provide an asymptotically correct 1 − α confidence interval (CI) for η(F c ) by finding
qα/2 and q1−α/2 such that
lim Pr qα/2 ≤ η(F c ) ≤ q1−α/2 = 1 − α,n→∞while estimated [qα/2 , q1−α/2 ] may show overcoverage given finite n and simulation effort. We wish tocharacterize input model error for a limited sample, not with a large sample assumption that permits a delta
method approach (Cheng and Holland 2004; Morgan et al. 2017). Further, we are interested in approaches
appropriate for input models using empirical distributions, as opposed to the assumption of parametric
distributions appearing in most prior work.
Barton and Schruben (2001) use the following direct bootstrap approach to estimate qα/2 and q1−α/2
nonparametrically. Suppose that the size-n sample data result in empirical distribution Fb0 . The bootstrap
assumption is that, given Fb0 , the quantity Fbb constructed by resampling n values from the original datawith replacement satisfies
D
η(Fbb ) − η(Fb0 ) −→ η(Fb0 ) − η(F c )(2)
as n → ∞. However, in the stochastic simulation context, η(F) cannot be computed exactly and can onlyRbe estimated by ȲR (F) = ∑ Yr (Fb0 )/R given R replications. Therefore, instead of (2) we hope
r=1

DȲR (Fbb ) − η(Fb0 ) −→ ȲR (Fb0 ) − η(F c )(3)

as n → ∞. However, in practice η(Fb0 ) in the left-hand-side of (3) is also unknown. Therefore, we use
ȲR0 (Fb0 ) to approximate η(Fb0 ) in (3) for some large R0 . If we define the α/2 and 1 − α/2 sample quantiles

Barton, Lam, and Songof ȲR (Fb1 ) − ȲR0 (Fb0 ), ȲR (Fb2 ) − ȲR0 (Fb0 ), . . . , ȲR (FbB ) − ȲR0 (Fb0 ) as τbα/2 and τb1−α/2 , respectively, then a 1 − αconfidence interval for η(F c ) is constructed asȲR (Fb0 ) − τb1−α/2 ≤ η(F c ) ≤ ȲR (Fb0 ) − τbα/2 .

(4)Note ȲR (Fb0 ) 6= ȲR0 (Fb0 ), if R < R0 . This scheme is justified when R0 is chosen big enough (more precisely,that R and the number of bootstrap iterations are also large, but these can be of smaller order than R0 ). Theinterval (4) corresponds to the so-called basic bootstrap in the statistics literature (Davison and Hinkley1997).
Alternately, one may attempt to directly use the bootstrap percentile interval
qbα/2 ≤ η(F c ) ≤ qb1−α/2 .

(5)where qbα/2 and qb1−α/2 are the α/2 and 1 − α/2 sample quantiles of ȲR (Fbb )’s. The validity of this scheme,
however, relies on choosing R = R0 and the assumption that the distribution of ȲR (Fbb )− ȲR0 (Fb0 ) is symmetric.
However, when R is small, (5) may give an incorrect coverage.
To illustrate the last point, consider the basic bootstrap in (4). To make our discussion more concrete,
we write τbp (ξ ) as the p sample quantile of a generic random variable ξ , so that τb1−α/2 defined in (4) canbe written as τb1−α/2 (ȲR (Fbb ) − ȲR0 (Fb0 )). Here the quantile is taken with respect to the randomness of Fbb
given Fb0 . The lower bound in (4) is thereforeȲR (Fb0 ) − τb1−α/2 (ȲR (Fbb ) − ȲR0 (Fb0 )).

(6)

If the distribution of ȲR (Fbb ) − ȲR0 (Fb0 ) is symmetric, then −τb1−α/2 (ȲR (Fbb ) − ȲR0 (Fb0 )) = τbα/2 (ȲR (Fbb ) −ȲR0 (Fb0 )), and (6) becomesȲR (Fb0 ) + τbα/2 (ȲR (Fbb ) − ȲR0 (Fb0 )).
(7)
In situations where no simulation is needed to approximate the output, ȲR (Fb0 ) is the same as ȲR0 (Fb0 ) andhence (7) is equal to τbα/2 (ȲR (Fbb )); similarly for the upper bound. This leads to the bootstrap percentile
interval. However, if we have R 6= R0 , then we should use (7), which unfortunately does not lead to a cleanbootstrap percentile formula.
Note that both the percentile and basic bootstrap intervals can have problems with coverage if thebootstrap distribution is highly asymmetric. See Davison and Hinkley (1997) for alternatives for this case.
2.2 Asymptotic Normality Based Interval
The bootstrap intervals discussed above can be compared to methods based on asymptotic normality directly.
b where theCheng and Holland (1997) suggest to use the latter to construct CIs focusing on parametric F,
standard error can be estimated using the bootstrap. More specifically, they use
s
s
2
σ
σ2
≤ η(F c ) ≤ ȲR0 (Fb0 ) + z1−α/2 V 2 +(8)
ȲR0 (Fb0 ) − z1−α/2 V 2 +R0
R0
as a 1 − α confidence interval, where z1−α/2 is the (1 − α/2)-quantile of the standard normal distribution.
The quantity V 2 is the variance contributed from the input noise, namely Var(η(Fb0 )), and σ 2 is the variancefrom the simulation noise that, at least asymptotically, is given by σ 2 (F) defined previously.
The variance V 2 can be approximated by the bootstrap input variance, Var(η(Fbb )), where Var(·) is withrespect to the sampling distribution of Fbb given Fb0 . Introducing the notation Yb,r ≡ Yr (Fbb ) for convenience,
we define

2

2
1 B
B
R
B
b
b
b
SSW = ∑b=1 ∑r=1 Yb,r − ȲR (Fb ) and SSB = ∑b=1 ȲR (Fb ) − ∑i=1 ȲR (Fi ) .
(9)
BBarton, Lam, and Song
Then an estimate of V 2 is

SSB
SSW
−
.
B − 1 BR(R − 1)

(10)This formula is an unbiased estimator for Var(η(Fbb )), which consists of Var(ȲR (Fbb )) estimated by the firstterm in (10), and adjusted for the bias introduced by E[σ 2 (Fbb )/R] that is estimated by the second term in (10).
When B and R are small, (10) may result in a negative value in which case Vb 2 = 0 (Ankenman and Nelson2012). On the other hand, the simulation variance σ 2 typically can be estimated more straightforwardly,
for instance by taking the sample variance from all the simulation runs. These quantities have been studiedin Ankenman and Nelson (2012), who focus on the ratio Var(η(Fbb ))/E[σ 2 (Fbb )/R] under the assumptionof homoscedastic simulation error, i.e., σ 2 (Fbb ) = σ 2 for all b.
3

SHRINKAGE ADJUSTMENT

To motivate our investigation, let us first compare the performance of the basic bootstrap (4) and theasymptotic-normality-based interval (8). Both of them require R0 simulation runs to estimate a point
estimate, and B bootstrap iterations, each with R simulation runs, to estimate the quantile or standard
error. The basic bootstrap requires in addition a separate simulation of size R to estimate another point
estimate, namely the first term in (4). For simplicity, let us consider the case where R0 is very big (inthe sense that the Monte Carlo error of the corresponding point estimate is outwashed). Then, givena large number of bootstrap samples, the overall computational efforts for the basic bootstrap and theasymptotic-normality-based approach are roughly comparable.
Under the above presumptions, the asymptotic-normality-based interval is approximately
ȲR0 (Fb0 ) − z1−α/2V ≤ η(F c ) ≤ ȲR0 (Fb0 ) + z1−α/2V,
which has a half-width z1−α/2V . On the other hand, the half-width of the basic bootstrap interval in (4)pdepends on the standard error of ȲR (Fbb ) − ȲR0 (Fb0 ), which is approximately z1−α/2 V 2 + σ 2 /R (recall that
R0 is assumed big). This half-width is longer than the asymptotic-normality-based approach in general,thus making the basic bootstrap less efficient. Our main investigation is to lift the construction of aquantile-based interval, using the basic or the percentile bootstrap, to the same level of efficiency as
the asymptotic-normality-based approach. This is motivated from documented studies that quantile-based
approach has better finite-sample coverage properties.
We consider two approaches to address the above efficiency issue. The first is a shrinkage operation
described in Davison and Hinkley (1997) that aims to remove the “excess variation” in ȲR (Fbb )’s. To
explain, recall how we have obtained (4). We started with a point estimate ȲR (Fb0 ), and we approximated the
distribution of ȲR (Fb0 ) − η(F c ) p
with a bootstrapped version ȲR (Fbb ) − η(Fb0 ). This ends up giving us a halfwidth of approximately z1−α/2 V 2 + σ 2 /R. Suppose that we could use ȲR0 (Fb0 ) as our point estimate, andthat we could efficiently bootstrap the corresponding distribution of ȲR0 (Fbb ) − η(Fb0 ), then we could shrink
the half-width to z1−α/2V (again, assuming R0 is a big number). But bootstrapping ȲR0 (Fbb ) − η(Fb0 ) directlyis computationally intense because R0 is big (this in particular also means placing more computational
effort than the asymptotic-normality-based approach). Instead, we adjust the samples of ȲR (Fbb ) − η(Fb0 ) sothat they exhibit the same variance as ȲR0 (Fbb ) − η(Fb0 ). To do so, we define
YbR (Fbb ) = cȲ¯R + (1 − c) ȲR (Fbb ), 1 ≤ b ≤ B,

(11)

where Ȳ¯R is the grand mean given by B1 ∑Bi=1 ȲR (Fbi ), and c satisfies(1 − c)2 =B
Var(η(Fbb )).
B − 1 Var(η(Fbb )) + σ 2 /R

(12)

Barton, Lam, and SongTo verify our claim, we rewrite the variance of (11) (under Fb0 ) as follows. Defining Ȳ¯R,−b as the sample
mean of ȲR (Fbi ), ∀i 6= b,Var(YbR (Fbb )) = Var cȲ¯R + (1 − c) ȲR (Fbb )c((B − 1)Ȳ¯R,−b + ȲR (Fbb )) + (1 − c) ȲR (Fbb )
= Var


B
c
c
(B − 1)Ȳ¯R,−b ++ (1 − c) ȲR (Fbb )
= Var
B
B


c
2
2
2
c (B − 1)¯
Var(ȲR (Fbb ))
Var
Ȳ
+
+
(1
−
c)
=
R,−b
B2
B

 

2 
c2 (B − 1)2 1
σ2
σ2
c
b
b
=
Var(η(Fb )) ++ (1 − c)Var(η(Fb )) +
+
B2
B−1
R
B
R

 2

2
2
2
c
c
c
2c(1 − c)σ−+++ (1 − c)2=
Var(η(Fbb )) +B B2 B2
B
R



2
B−1 1
σ
=
(1 − c)2+
Var(η(Fbb )) +.
B
B
R
Thus, if we choose c satisfying



σ2
B−1 1
+
Var(η(Fbb )) += Var(η(Fbb )),(1 − c)2
B
B
R

(13)

then we match the variance of YbR (Fbb ) with ȲR0 (Fbb ) for large R0 . The relation (13) gives(1 − c)2 =B
Var(η(Fbb ))
1
−
.
2
b
B − 1 Var(η(Fb )) + σ /R B − 1(14)We have ignored the small second term in (14) when choosing the c in (12). The above derivation follows
the idea in Davison and Hinkley (1997), but there they assume a linear model and in that case (12) canbe shown to give a variance that is unbiased under the true distribution. Here, since we have a nonlinear
simulation model to begin with, we are contended that the variance of YbR (Fbb ) matches that of the bootstrap
variance Var(η(Fbb )).
Note that the choice of c depends on Var(η(Fbb )) and σ 2 , which are unknown in general. We plug inestimates for them using ideas similar to (10), giving
s


B
SSW
1 − cb = max 0,−.
(15)
B − 1 R (R − 1) SSB
See (9) for the expressions for SSW and SSB .
To sum up, in the shrinkage basic bootstrap, we first obtain a point estimate ȲR0 (Fb0 ) by running R0
simulation runs under the empirical distribution Fb0 . We then obtain each ȲR (Fbb ) by resampling from Fb0
and running R simulation runs, like in the basic bootstrap before. Then we form the new bootstrap output
sample YbR (Fbb ) using (11) with c defined in (15). The 1 − α CI for η(F c ) is constructed as
ȲR0 (Fb0 ) − τ̃1−α/2 ≤ η(F c ) ≤ ȲR0 (Fb0 ) − τ̃α/2 ,where τ̃1−α/2 and τ̃α/2 are now the 1 − α/2 and α/2 sample quantiles of YbR (Fbb ) − ȲR0 (Fb0 ).

(16)

Barton, Lam, and SongThe same idea could be applied to construct a percentile interval. Starting from the shrinkage basic
bootstrap above, we have, using the notations in (6),ȲR0 (Fb0 ) − τb1−α/2 (YbR (Fbb ) − ȲR0 (Fb0 )) ≤ η(F c ) ≤ ȲR0 (Fb0 ) − τbα/2 (YbR (Fbb ) − ȲR0 (Fb0 ))being a valid CI that has a half-width roughly z1−α/2V . Now, if ȲR (Fbb )− ȲR0 (Fb0 ) is approximately symmetric,
so is YbR (Fbb ) − ȲR0 (Fb0 ), and thus we haveȲR0 (Fb0 ) − τb1−α/2 (YbR (Fbb ) − ȲR0 (Fb0 )) ≈ ȲR0 (Fb0 ) + τbα/2 (YbR (Fbb ) − ȲR0 (Fb0 )) = τbα/2 (YbR (Fbb )).
The upper bound can be derived similarly. Hence, in the shrinkage percentile bootstrap, the CI isq̃α/2 ≤ η(F c ) ≤ q̃1−α/2 ,(17)where q̃α/2 and q̃1−α/2 are the α/2 and 1 − α/2 sample quantiles of YbR (Fbb ). Note that, in either the basicor the percentile bootstrap, computing cb and constructing the shrinkage interval requires no additional
simulation runs.
Below, we present the complete algorithm to construct the shrinkage basic and percentile CIs.
Algorithm: Shrinkage CIs0
1. Run R0 replications of the simulator using Fb0 as the input model to obtain ȲR0 (Fb0 ) = ∑Rr=1
Yr (Fb0 )/R0 .
2. For b = 1, 2, . . . , B(a) Generate Fbb by resampling Fb0 n times.
(b) Using Fbb , run R replications, Y1 (Fbb ),Y2 (Fbb ), . . . ,YR (Fbb ), and compute ȲR (Fbb ) = ∑Rr=1 Yr (Fbb )/R.
3. Using Yr (Fbb ), b = 1, 2 . . . , B, r = 1, 2, . . . , R, calculate cb from (15).
4. For b = 1, . . . , B, compute YbR (Fbb ) = cbȲ¯R + (1 − cb) ȲR (Fbb ), where Ȳ¯ = ∑Bb=1 YR (Fbb )/B.
5. CI construction:(a) (Basic CI) Find the empirical α/2 and 1 − α/2 quantiles of YbR (Fb1 ) − ȲR0 (Fb0 ), YbR (Fb2 ) −ȲR0 (Fb0 ), . . . , YbR (FbB ) − ȲR0 (Fb0 ), τ̃α/2 and τ̃1−α/2 , respectively, and construct the basic shrinkage
CI in (16).
(b) (Percentile CI) Find the empirical α/2 and 1 − α/2 quantiles of YbR (Fb1 ), YbR (Fb2 ), . . . , YbR (FbB ),q̃α/2 and q̃1−α/2 , respectively, and construct the percentile shrinkage CI in (17).

4HIERARCHICAL BOOTSTRAP

Next we discuss a hierarchical bootstrap approach that accounts for both the extrinsic and intrinsic errors inthe basic and the percentile bootstraps that, like the shrinkage, allows a comparable computational effort andstatistical efficiency as the normality-based approach. We start by discussing the basic bootstrap. Again,
let us assume that R0 is large. Consider the scheme
ȲR0 (Fb0 ) − τb1−α/2 (ȲR0 (Fbb ) − η(Fb0 )) ≤ η(F c ) ≤ ȲR0 (Fb0 ) − τbα/2 (ȲR0 (Fbb ) − η(Fb0 )).
As discussed before, though this is valid, it requires huge computational effort since approximatingτb1−α/2 (ȲR0 (Fbb ) − η(Fb0 )) and τbα/2 (ȲR0 (Fbb ) − η(Fb0 )) requires many iterations each with R0 simulation runs.
Our hierarchical bootstrap approach aims to replace these quantities with τb1−α/2 (ȲR (Fbb ) − η(Fb0 )) andτbα/2 (ȲR (Fbb ) − η(Fb0 )), and subsequently de-bias them back to τb1−α/2 (ȲR0 (Fbb ) − η(Fb0 )) and τbα/2 (ȲR0 (Fbb ) −η(Fb0 )). In other words, we want to findτb1−α/2 (ȲR (Fbb ) − η(Fb0 )) − τb1−α/2 (ȲR0 (Fbb ) − η(Fb0 )) = τb1−α/2 (ȲR (Fbb )) − τb1−α/2 (ȲR0 (Fbb ))(18)

Barton, Lam, and Songand
τbα/2 (ȲR (Fbb ) − η(Fb0 )) − τbα/2 (ȲR0 (Fbb ) − η(Fb0 )) = τbα/2 (ȲR (Fbb )) − τbα/2 (ȲR0 (Fbb )).

(19)

Since the differences in (18) and (19) stem from the difference in simulation size, we propose the following
scheme. Consider the existing bootstrapped outputs Yb∗ = ȲR (Fbb ), b = 1, . . . , B, where the simulation runs
∗ = Y ∗ + ε ∗ , r = 1, 2, . . . , R, where ε ∗ isfrom each resample are averaged. For each Yb∗ , compute Yb,r
b,r
b,r
b
2
2
2
∗
∗
b
b
b
sampled from εb,r ∼ N(0, S (Fb )) and S (Fb ) is an estimator of σ (Fb ). Then take the average of Yb,r acrossb∗∗r, for each b, to obtain Yb∗∗ . Use these Yb∗∗ to form the quantiles, say qb∗∗
α/2 . Correspondingly,
1−α/2 and q
∗∗
∗
∗∗
∗
∗
∗
let qb1−α/2 and qbα/2 be the quantiles of Yb . We then use qb1−α/2 − qb1−α/2 and qbα/2 − qb∗α/2 as estimates of(18) and (19). In this way, we approximate τb1−α/2 (ȲR0 (Fbb ) − η(Fb0 )) withb∗1−α/2 )τb1−α/2 (ȲR (Fbb ) − ȲR0 (Fb0 )) − (bq∗∗1−α/2 − qand τbα/2 (ȲR0 (Fbb ) − η(Fb0 )) withb∗α/2 ).
τbα/2 (ȲR (Fbb ) − ȲR0 (Fb0 )) − (bq∗∗α/2 − qA 1 − α basic hierarchical bootstrap CI isb∗1−α/2 ) ≤ η(F c )ȲR0 (Fb0 ) − τb1−α/2 (ȲR (Fbb ) − ȲR0 (Fb0 )) + (bq∗∗1−α/2 − qb∗α/2 ).
≤ ȲR0 (Fb0 ) − τbα/2 (ȲR (Fbb ) − ȲR0 (Fb0 )) + (bq∗∗α/2 − qEquivalently, this iscb
2ȲR0 (Fb0 ) − 2bq∗1−α/2 + qb∗∗q∗α/2 + qb∗∗
1−α/2 ≤ η(F ) ≤ 2ȲR0 (F0 ) − 2b
α/2 .

(20)

b∗1−α/2 and qb∗∗b∗α/2 to capture the inflation in the quantiles due
Roughly speaking, (20) uses qb∗∗
1−α/2 − qα/2 − qto the noise from R simulation runs, using a normal approximation for the simulation noise.
Similar idea applies to the percentile bootstrap. In this case, we need to approximate τb1−α/2 (ȲR0 (Fbb )) andb∗1−α/2 and qb∗∗b∗α/2 as estimates of τb1−α/2 (ȲR (Fbb )) − τb1−α/2 (ȲR0 (Fbb ))τbα/2 (ȲR0 (Fbb )). Using qb∗∗
1−α/2 − qα/2 − qand τbα/2 (ȲR (Fbb )) − τbα/2 (ȲR0 (Fbb )), as for (18) and (19), and qb∗1−α/2 and qb∗α/2 for τb1−α/2 (ȲR (Fbb )) andτbα/2 (ȲR (Fbb )), we approximate τb1−α/2 (ȲR0 (Fbb )) with qb∗
− (b
q∗∗
− qb∗
) = 2bq∗− qb∗∗
, and
1−α/2

1−α/2

1−α/2

1−α/2

1−α/2

similarly τbα/2 (ȲR0 (Fbb )) with 2bq∗α/2 − qb∗∗
α/2 . A 1 − α hierarchical bootstrap percentile CI isc2bq∗α/2 − qb∗∗q∗1−α/2 − qb∗∗
α/2 ≤ η(F ) ≤ 2b1−α/2 .

(21)

Below, we present the complete algorithm to construct the hierarchical basic and percentile CIs.
Algorithm: Hierarchical Bootstrap CIs
0
Yr (Fb0 )/R0 .
1. Run R0 replications of the simulator using Fb0 as the input model to obtain ȲR0 (Fb0 ) = ∑Rr=1
2. For b = 1, 2, . . . , B(a) Generate Fbb by resampling Fb0 n times.
(b) Using Fbb , run R replications, Y1 (Fbb ),Y2 (Fbb ), . . . ,YR (Fbb ), and compute ȲR (Fbb ) = ∑Rr=1 Yr (Fbb )/R
and S2 (Fbb ) = ∑Rr=1 (Yr (Fbb ) − ȲR (Fbb ))2 /(R − 1).
∗ ∼ N(0, S2 (Fbb )) and compute Y ∗ = ȲR (Fbb ) + ε ∗ .
(c) For r = 1, 2, . . . , R, generate εbr
br
br
R
∗∗
∗
(d) Let Yb = ∑r=1 Ybr /R.
3. Find the empirical α/2 and 1 − α/2 quantiles of ȲR (Fbb ), qb∗α/2 and qb∗1−α/2 , respectively.

Barton, Lam, and Songb∗∗
4. Find the empirical α/2 and 1 − α/2 quantiles of Yb∗∗ , qb∗∗
1−α/2 , respectively.
α/2 and q
5. CI construction:(a) (Basic CI) Construct the basic hierarchical boostrap CI in (20).
(b) (Percentile CI) Construct the percentile hierarchical boostrap CI in (21).
Note that, like the shrinkage procedure, the hierarchical approach for both the basic and the percentile
bootstraps does not require additional simulation runs.
5

EMPIRICAL COMPARISON OF METHODSIn this section, we compare the empirical coverage probabilities of the CIs we discussed: basic bootstrap
CI (4), percentile bootstrap CI (5), basic shrinkage CI (16), percentile shrinkage CI (17), basic hierarchicalbootstrap CI (20), and percentile hierarchical bootstrap CI (21). We also present the results of asymptotic
normality-based CI (8) and nominal t-interval that only accounts for intrinsic error for comparison.
We simulate an M/M/1/10 system whose inter-arrival and service time distributions have rates 0.8 and1, respectively. We assume these distributions are unknown, but n observations from each distribution
are available. The performance measure of interest is the steady-state expected time in system. For eachsimulation replication, we average 200 waiting times while deleting the initial 400 observations for warm-up.
We tested two different real-world sample sizes, n = 100 and n = 1,000. For the number of replications
at each bootstrap sample and the number of bootstrap samples, we tested R = 2 and R = 40, and B = 100
and B = 1,000, respectively. Table 1 show these eight combinations of (n, R, B). All cases are repeated
1,000 times and the target coverage was set to 1 − α = 0.9.
Figure 1 shows empirical coverage probability and the average width for each CI for each test case.
Notice that the nominal t-interval shows poor coverage when R is large as it fails to account for input
model uncertainty. The percentile bootstrap CI has a better coverage than the basic bootstrap CI in allcases. The former shows overcoverage when R is small (Case 1, 2, 5, and 6), while the latter showsundercoverage when n is small (Case 1, 2, 3, and 4). Notice that the basic and percentile bootstrap CIs arewider than other CIs (except for the t-intervals) as it does not compensate for the inflation due to finite
R. The undercoverage of the basic bootstrap CIs, yet with long interval widths, hints that the positionsof the CIs highly vary. When R and B are small, asymptotic normal CIs undercover on average as input
uncertainty is poorly estimated by (10). This becomes more severe when n is large as input uncertainty
becomes small relative to intrinsic error variance, which may make (10) negative. As mentioned earlier,
we set Vb 2 = 0 in this case, i.e., the resulting asymptotic normal CI accounts only for the intrinsic noise.
Similar observations can be made for the basic and percentile shrinkage CIs. Small B and R may resultin cb = 1 from (12), which makes all YbR (Fbb ) = Ȳ¯R . Then, both basic and percentile shrinkage CIs become
degenerate. For both n = 100 and n = 1,000, the percentile shrinkage CIs provide good coverage when B
and R are large, however, the basic shrinkage CIs undercover when n = 100. Hierarchical boostrap CIs are
the most robust under small R and B. In particular, when R is small, the hierarchical bootstrap inflates theCIs enough to provide the correct coverage, but not as much as the basic and percentile bootstrap CIs do.
6

CONCLUSION

Basic and percentile bootstrap CIs tend to overcover or exhibit high positioning variabilities because of
intrinsic noise that inflates the CIs given any finite R. Asymptotic normality-based CIs, which is typically
Table 1: Combinations of (n, R, B) used for testing the CIs in Figure 1.
Case
n
R
B

1
100
2
100

2
100
2
1,000

3
100
40
100

4
100
40
1,000

5
1,000
2
100

6
1,000
2
1,000

7
1,000
40
100

8
1,000
40
1,000

1.0
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0.0
-0.1
-0.2
-0.3
-0.4
-0.5
-0.6
-0.7
-0.8
-0.9
-1.0

-8.0
-6.0
-4.0
-2.0
0.0

1

Basic bootstrap
Percentile shrinkage

2

3

4

5

6

7

8
2.0
4.0
6.0

CI width

Coverage probability

Barton, Lam, and Song

8.0

Percentile bootstrap
Basic hierarchical

Asymptotic normal
Percentile hierarchical

Basic shrinkage
t-interval

Figure 1: Empirical coverage probability and the average width for the eight CIs in comparison computedfrom 1,000 repeats. The target coverage is 1 − α = 0.9.
b becomes asymptotically normallyused in the literature to correct for the overcoverage, relies on that η(F)distributed when real-world sample size n increases. Although it appears to work well for the example inb We propose four new
Section 5, bootstrap-based CIs tend to be more robust to non-normality of η(F).
cbootstrap-based CIs for η(F ). The percentile shrinkage CI perform well when n, R, and B are large while
the basic hierarchical bootstrap CI show good coverage across all combinations of (n, R, B).
For shrinkage CIs, we assumed homoscedasticity of intrinsic error for any Fbb . Relaxing this assumption
could increase their coverage probability. From the experiment, the basic hierarchical bootstrap CI improves
the basic bootstrap CI even for small n, while the percentile hierarchical bootstrap does not. Further analysis
is required to understand their asymptotic properties.
REFERENCES
Ankenman, B. E., and B. L. Nelson. 2012. “A Quick Assessment of Input Uncertainty”. In Proceedings ofthe 2012 Winter Simulation Conference, edited by C. Laroque et al., 1–10. Piscataway, New Jersey:
IEEE.
Barton, R. R., and L. W. Schruben. 1993. “Uniform and Bootstrap Resampling of Empirical Distributions”.
In Proceedings of the 1993 Winter Simulation Conference, edited by G. W. Evans et al., 503–508.
Piscataway, New Jersey: IEEE.
Barton, R. R., and L. W. Schruben. 2001. “Resampling Methods for Input Modeling”. In Proceedings of the2001 Winter Simulation Conference, edited by B. A. Peters et al., 372–378. Piscataway, New Jersey:
IEEE.
Barton, R. R. 2007. “Presenting a More Complete Characterization of Uncertainty: Can it be done?”. In
Proceedings of the 2007 INFORMS Simulation Society Research Workshop, edited by S. C. et al.
Barton, R. R., B. L. Nelson, and W. Xie. 2014. “Quantifying input uncertainty via simulation confidence
intervals”. INFORMS Journal on Computing 26(1):74–87.
Cheng, R. C. H. 1994. “Selecting Input Models”. In Proceedings of the 1994 Winter Simulation Conference,edited by J. D. Tew et al., 184–191. Piscataway, New Jersey: IEEE.
Cheng, R. C. H., and W. Holland. 1997. “Sensitivity of Computer Simulation Experiments to Errors in
Input Data”. Journal of Statistical Computation and Simulation 57(1-4):219–241.

Barton, Lam, and SongCheng, R. C. H., and W. Holland. 2004. “Calculation of Confidence Intervals for Simulation Output”. ACM
Transactions on Modeling and Computer Simulation 14:344–362.
Davison, A. C., and D. V. Hinkley. 1997. Bootstrap Methods and their Application. Cambridge Series in
Statistical and Probabilistic Mathematics. Cambridge University Press.
Morgan, L. E., B. L. Nelson, A. C. Titman, and D. J. Worthington. 2017. “Detecting Bias due to Input
Modelling in Computer Simulation”. In Proceedings of the 2017 Winter Simulation Conference, edited
by W. K. V. Chan et al., 1974–1985. Piscataway, New Jersey: IEEE.
Song, E., B. L. Nelson, and C. D. Pegden. 2014. “Advanced Tutorial: Input Uncertainty Quantification”. In
Proceedings of the 2017 Winter Simulation Conference, edited by A. Tolk et al., 162–176. Piscataway,
New Jersey: IEEE.
Song, E., and B. L. Nelson. 2015. “Quickly Assessing Contributions to Input Uncertainty”. IIE Transactions 47(9):893–909.
ACKNOWLEDGEMENT
Computations for this research were performed on the Penn State University’s Institute for CyberScience
Advanced CyberInfrastructure (ICS-ACI). We gratefully acknowledge support from the National Science
Foundation under grants CMMI-1542020 and CAREER CMMI-1653339/1834710.
AUTHOR BIOGRAPHIESRUSSELL R. BARTON is Distinguished Professor of Supply Chain and Information systems andProfessor of Industrial Engineering at the Pennsylvania State University. He received a B.S. degree inelectrical engineering from Princeton University and M.S. and Ph.D. degrees in operations research from
Cornell University. He serves as INFORMS Vice President for Sections and Societies. He is a fellowof IIE and a Certified Analytics Professional. His research interests include applications of statisticaland simulation methods to system design and to product design, manufacturing and delivery. His email
address is rbarton@psu.edu.
HENRY LAM is an Associate Professor in the Department of Industrial Engineering and Operations
Research at Columbia University. He received his Ph.D. degree in statistics at Harvard University in 2011,
and was on the faculty of Boston University and the University of Michigan before joining Columbia in
2017. His research focuses on Monte Carlo simulation, risk and uncertainty quantification, and stochastic
optimization. His email address is khl2114@columbia.edu.
EUNHYE SONG is Harold and Inge Marcus Early Career Assistant Professor in the Department of Industrial and Manufacturing Engineering at the Penn State University. Her research interests include simulation
design of experiments, simulation uncertainty and risk quantification, optimization via simulation under
input model risk and large-scale discrete optimization via simulation. Her email address is eus358@psu.edu.