Proceedings of the 2022 Winter Simulation Conference
B. Feng, G. Pedrielli, Y. Peng, S. Shashaani, E. Song, C.G. Corlu, L.H. Lee, E.P. Chew, T. Roeder, andP. Lendermann, eds.
DISTRIBUTIONAL INPUT UNCERTAINTY
Motong Chen
Zhenyuan Liu
Henry Lam
Department of Industrial Engineering andOperations Research
Columbia University
500 West 120th Street
New York, NY 10027, USA

ABSTRACTThe vast majority of the simulation input uncertainty literature focuses on estimating target output quantities
that are real-valued. However, outputs of simulation models are random and real-valued targets essentiallyserve only as summary statistics. In this paper, we study the input uncertainty problem from a distributionalview, namely we construct simultaneous confidence bands for the entire output distribution function. Our
approach utilizes a novel test statistic that consists of the supremum of the sum of a Brownian bridge
and a mean-zero Gaussian process whose covariance function is characterized by the influence functionof the true output distribution function, which generalizes the Kolmogorov-Smirnov statistic to accountfor input uncertainty. We demonstrate how subsampling helps estimate the covariance function of theGaussian process, thereby leading to an implementable estimation of the quantile of the test statistic anda statistically valid confidence band. Finally, we present some supporting numerical experiments.
1

INTRODUCTION

Stochastic simulation is widely used for operational decision-making and scientific modeling (Law et al.
(2007), Banks et al. (2014)). It consists of the generation of random input variates, fed into the system
logic, to produce outputs that are subsequently analyzed for downstream tasks. In many real-world settings,
the distributions that govern the input variates are unknown and need to be calibrated from external data. In
this case, proper output analysis needs to account for both the Monte Carlo noises in running the simulation
and also the input data noises that propagate to the outputs. This issue, which is beyond the handling of
only the Monte Carlo noises in traditional output analysis, is known commonly as the input uncertainty
problem in the simulation literature. For an overview on this topic, see, e.g., Henderson (2003), Barton(2012), Song et al. (2014), Lam (2016), Corlu et al. (2020), Barton et al. (2022), and Nelson (2013)
Chapter 7.
More concretely, the goal of handling input uncertainty can often be cast as the construction ofconfidence intervals that are statistically valid under both sources of noises from Monte Carlo and input
data. A closely related task is variance estimation, or the estimation of the different components of thevariance that attribute to each source of noise. There are several major methods in the literature. The
first is bootstrap resampling, which includes the variance bootstrap (Cheng and Holloand (1997), Songand Nelson (2015)), quantile-based bootstrap (Barton and Schruben (1993), Barton and Schruben (2001)),
and enhancements such as the metamodel-assisted bootstrap (Barton et al. (2014)) and subsampling (Lamand Qian (2021)). The second line of methods uses the delta method directly which involves estimating
gradient information of the model with respect to the inputs. This includes the two-point method (Cheng
and Holland (1998), Cheng and Holland (2004)), regression approach (Lin et al. (2015), Song and Nelson

Chen, Liu, and Lam(2019)), and subsampled infinitesimal jackknife (Lam and Qian (2019), Lam et al. (2022)). Other recent
approaches include the use of robust optimization (Ghosh and Lam (2019)) which relates to the empirical
likelihood method (Lam and Qian (2016), Lam and Qian (2017)). Finally, there is also a stream of work on
Bayesian approaches that rely on posterior sampling (Chick (2001), Zouaoui and Wilson (2003), Zouaoui
and Wilson (2004), Xie et al. (2014)).
Despite the growing literature in input uncertainty, most, if not all, of the works have been focusing
on the estimation of target quantities that are real-valued. However, in stochastic simulation the outputs
are random (hence the Monte Carlo noise) and the real-valued target quantities often only serve as asummary statistic. For instance, suppose the simulation model represents a queuing system and we areinterested in the waiting time or queue length at a certain time horizon. Then, in the input uncertainty
literature, existing approaches would consider a summary target such as the expected waiting time. Suchtype of summary statistic is arguably convenient but, depending on our downstream decision-making tasks,
could hide valuable information (in this example this could be the quantiles or tails of the waiting time
distribution). Having an approach to understand the interplay of Monte Carlo and input uncertainty at
the distributional level would therefore address this issue. In a sense, doing so would capture both theepistemic uncertainty, i.e., the uncertainty coming from the data noise, and aleatory uncertainty, i.e., theintrinsic uncertainty from the randomness in the stochastic model.
Our main contribution of this paper is to propose an approach to construct a confidence band for theentire output distribution that accounts for both Monte Carlo and input data noises. Our approach relies
on a novel statistic that comprises of the supremum of a sum of a Brownian bridge and a mean-zero
Gaussian process with a covariance structure governed by the so-called influence function of the target
output distribution. In this statistic, the Brownian bridge signifies the Monte Carlo noises and resembles
the limit of the Kolmogorov-Smirnov (KS) statistic. The mean-zero Gaussian process, on the other hand,signifies the input data noises, where the influence function that controls the covariance structure can be
regarded as the derivative of the target output distribution with respect to the input parameters and carrysensitivity information like in the standard delta method. Our approach can be viewed as a generalization of
the KS statistic, which can be used to construct confidence band under a single source of randomness from
data, to the simulation input uncertainty setting where both epistemic and aleatory uncertainties need to beaccounted for. In fact, we will illustrate that when there is input uncertainty, using the standard KS machinery
to construct confidence bands for the output distribution would lead to a systematic under-coverage. The
additional Gaussian process serves to inflate the variability to lift up the coverage by the right amount.
With the test statistic, we propose a subsampling approach to efficiently estimate the covariancefunction of the Gaussian process. Note that the influence function described above, while conceptually
understandable, is difficult to use in practice because it requires computing essentially infinite-dimensional
derivatives. Thus the bootstrap is used in lieu of the direct estimation of the influence function. However,standard variance bootstrap is inefficient because of the entanglement of the Monte Carlo and input datanoises, which necessitates washing away the Monte Carlo noises via a large amount of simulation runs. Thisphenomenon has been pointed out in Lam and Qian (2021) even when only real-valued target quantities
are considered. To address this computational challenge, we generalize the idea in Lam and Qian (2021)
to devise a subsampling method where we draw a resample of smaller size than the original input data
size, at the outer layer of a nested simulation scheme used to bootstrap estimate the covariance function.
The rest of this paper is as follows. Section 2 introduces the input uncertainty problem for the stochastic
simulation output distribution. Section 3 briefly discusses the theory of the input uncertainty quantification
for the distribution function. Section 4 then explains and presents the algorithms for estimating the covariancefunction of the Gaussian process and constructing the confidence bands. Finally Section 5 displays some
supporting numerical results.

Chen, Liu, and Lam
2

DISTRIBUTIONAL INPUT UNCERTAINTY PROBLEM

Consider a stochastic simulation model where each run outputs a target quantity, say Y , that follows a
distribution Q. This simulation model can depend on m input distributions P := (P1 , . . . , Pm ), so we may write
Q = Q(P) which highlights that the output distribution depends on the input distributions. For example,
think of a queuing system with two input distributions where P1 denotes the inter-arrival time distribution
and P2 denotes the service time distribution. The quantity Y can be the average waiting time of the first
10 customers in a queue, which follows some distribution Q(P).
Suppose we face input uncertainty for the stochastic simulation model, i.e., the true input distributions
are unknown. Instead, we have a finite number of data points for each distribution, say {Xi, j : i = 1, . . . , m, j =1, . . . , ni } where Xi, j , j = 1, . . . , ni are i.i.d. random variables from Pi and they are independent for different
i’s. In order to obtain an approximate distribution of Y , we need to drive the simulation model using an
estimate of P. A natural choice is the empirical distribution P̂ := (P̂1 , . . . , P̂m ) constructed from the dataP̂i (·) =1 ni∑ δXi, j (·),
ni j=1

where δx denotes the degenerate Dirac measure at x. Then, pretending that the empirical distributions areaccurate enough, we use them to drive Re simulation runs to obtain Re output samples such as the average
waiting time.
The output estimator we just mentioned contains statistical noises coming from both the external data
used to calibrate the input distributions, and the Monte Carlo noise in running the simulation model. To
conduct statistically valid inference for Y , we need to quantify both sources of noises. Existing literature
mainly focuses on inferring a real-value performance measure, e.g., the expected value E[Y ]. In this paper,
we are interested in the inference of the entire output distribution which, depending on the downstream
decision-making tasks, could give more complete and useful information about the aleatory uncertainty ofthe output.
More precisely, our goal is to construct simultaneous confidence bands for the output distributionfunction. For convenience, let Q(t, P) = PP (Y ≤ t) denote the true output distribution function where thesubscript P makes clear that Y is the output generated using the true unknown input distributions P. We
would like to find two functions L(t) and U(t) such that
P(L(t) ≤ Q(t, P) ≤ U(t) ∀t ∈ R) → 1 − α,(1)where 1 − α is a pre-specified confidence level, e.g., 97.5%, where the limit is taken as ni and Re get large.
To get a sense of how one can construct these L(t) and U(t), note that when we do not have noises
in calibrating the input distribution (i.e., suppose we know the input distributions completely), the only
uncertainty is from the Monte Carlo runs, and in this case we can use, for instance, the KS statistic. Suppose
the true distribution function Q(t, P) is continuous in t. Let Q̂(t, P) be the output empirical distribution
constructed from Re simulation runs driven by the input distribution P, i.e.,Q̂(t, P) =1 Re
∑ I(Yr ≤ t),
Re r=1

where Yr , r = 1, . . . , Re denote i.i.d. Monte Carlo outputs from the simulation model. Then the confidenceband (1) is accomplished with
q1−αL(t) = Q̂(t, P) − √ ,Re

q1−αU(t) = Q̂(t, P) + √ ,Rewhere q1−α is the (1 − α)-quantile of the asymptotic of the KS statistic, namely sups∈[0,1] |BB(s)| whereBB(·) is a standard Brownian bridge. The above L and U are piece-wise constant and can be computed

Chen, Liu, and Lam
readily. Thus, when there is no input uncertainty, the problem straightforwardly reduces to the classicalgoodness-of-fit problem.
However, when the input distributions need to be calibrated from external data, the above L and U
would under-cover the true output distribution Q(t, P), because they have ignored this additional source of
input noise. To remedy this, we need to properly inflate the variability of the considered statistic. This isdiscussed next.
3

A STATISTIC TO QUANTIFY DISTRIBUTIONAL INPUT UNCERTAINTY

We describe our procedure to construct confidence band for the entire output distribution. Along the way,
we explain the statistical underpinnings and briefly discuss the mathematical theory.
Let n = (1/m) ∑m
i=1 ni be the average data size. We assume the data sizes are balanced, i.e.,
lim ni /n = βi > 0, ∀i = 1, . . . , m.

(2)

n→∞

Then we have the following approximation results when n and Re are large:d G(t)
BB(Q(t, P))√,Q̂(t, P̂) − Q(t, P) = [Q(t, P̂) − Q(t, P)] + [Q̂(t, P̂) − Q(t, P̂)] ≈ √ +
n
Re

(3)

and consequentlydG(t) BB(Q(t, P))√
sup |Q̂(t, P̂) − Q(t, P)| ≈ sup √ +
,
(4)
n
Re
t∈R
t∈R
√
√
where G(t)/ n is an approximation of Q(t, P̂) − Q(t, P) and BB(Q(t, P))/ Re is an approximation ofQ̂(t, P̂) − Q(t, P̂). Here, BB(Q(t, P)) is a Brownian bridge with changed time scale Q(t, P),t ∈ R, and G(t)is a Gaussian process independent of BB(Q(t, P)) with domain on R, mean zero, and covariance structure
given bym
1
Cov(G(t), G(s)) = ∑ Cov(IFi (t, Xi ; P), IFi (s, Xi ; P)),i=1 βiwhere Xi follows the distribution Pi and IFi (t, x; P) is the so-called influence function of Q(t, P) for the i-th
input distribution defined as
Q(t, P1 , . . . , Pi−1 , (1 − ε)Pi + εδx , Pi+1 , . . . , Pm ) − Q(t, P).
ε↓0εIFi (t, x; P) = lim

The influence function IFi (t, x; P) signifies the functional gradient or the Gateaux derivative of Q(t, P)
with respect to the i-th input distribution Pi (see Wasserman (2006) Section 2.3 and Serfling (2009)
Chapter 6). It has mean zero under Pi , i.e., EPi [IFi (t, Xi ; P)] = 0. Now we briefly explain why theapproximation in (3) holds. For the Brownian bridge part, when n is large, we have P̂ ≈ P. Therefore,√dQ̂(t, P̂) − Q(t, P̂) ≈ Q̂(t, P) − Q(t, P) ≈ BB(Q(t, P))/ Re where the second approximation is due to the
Donsker Theorem (Van der Vaart (2000) Theorem 19.3). Then we explain how to obtain the Gaussian
process G(t), especially how we obtain the form of its covariance function. Since the influence function
can serve as the first order Taylor expansion of the statistical functional, we can obtainm Z

Q(t, P̂) − Q(t, P) ≈ ∑

i=1

1 ni
∑ IFi (t, Xi, j ; P),
i=1 ni j=1
m

IFi (t, x; P)d(P̂i (x) − Pi (x)) = ∑where the equality follows from the mean zero property of IFi (t, x; P) under Pi . Now using the independence
among Xi, j ’s, balanced data assumption in (2) and multivariate central limit theorem, we can easily show
that
!
m
√
1 ni
f .d.
n ∑ ∑ IFi (t, Xi, j ; P) → G(t),t ∈ R,
(5)
i=1 ni j=1

Chen, Liu, and Lam
f .d.

where → means any finite dimensional distribution of the left hand side (LHS) process converges in
distribution to the corresponding finite dimensional distribution of the right hand side process. The
argument is rigorized via proper checking of tightness, control of remainder terms in the linearization, andasymptotic independence between the components elicited by the Monte Carlo and input data noises.
Based on the approximation in (4), if we can estimate the (1 − α)-quantile ofG(t) BB(Q(t, P))√
sup √ +nRet∈R
say q̃1−α , then a statistically valid pair of confidence bands is given byL(t) = Q̂(t, P̂) − q̃1−α ,U(t) = Q̂(t, P̂) + q̃1−α ,(6)which are readily computable since they are piece-wise constant.
4

CONFIDENCE BAND CONSTRUCTION

In view of the approximation (4), in order to construct a (1 − α)-level confidence band for Q(t, P), the key
is to estimate the (1 − α)-quantile of the random variable
G(t) BB(Q(t, P))√
sup √ +nRe
t∈R
Exact simulation is hard to implement since it contains the supremum over the whole real line. To thisend, we discretize t on a grid over the real line, say t1 , . . . ,tk and replace the supremum over t ∈ R withthe maximum over t ∈ {t1 , . . . ,tk }, i.e.,
G(ti ) BB(Q(ti , P)) d
G(t) BB(Q(t, P))√√√ +≈ sup √ +.
i=1,...,k
n
n
Re
Re
t∈R
max

(7)

If we simulate the quantile of the LHS random variable in (7) and construct the confidence band like (6),
we actually end up giving a statistically valid confidence band only at these discrete values. But when thediscretization scale is refined enough, this discretized confidence band becomes an accurate approximationand its natural extrapolation should cover the true distribution function in R with probability close to
the nominal level. Therefore, we next aim at estimating the quantile of the LHS random variable in (7).
√BB(Q(ti , P)) can be easily approximated by BB(Q̂(ti , P̂)). But approximating G(ti ) or equivalently G(ti )/ n
is more complicated. Note that G(t) is a mean-zero Gaussian process. In order to simulate G(t), it suffices
to know its covariance function. However, the covariance function is hard to obtain since it involves theinfluence function which is not easy to compute in general. In this paper, we resort to the bootstrap toestimate the covariance function. In particular we will use subsampling to relieve the computation burden.
4.1 Subsampling to Estimate Covariance Function
Recall that we have ni data points Xi, j , j = 1, . . . , ni for the i-th input distribution Pi . We specifiy a subsample
ratio θ ∈ (0, 1], and compute si = ⌊θ ni ⌋, ∀i = 1, . . . , m as the subsample sizes. Based on the approximation(3), we havenCov(Q(t, P̂), Q(s, P̂)) ≈ Cov(G(t), G(s)).
(8)
At this point, we could consider using the bootstrap to approximate the LHS. However, note that thebootstrap requires a nested simulation where in the outer layer, we resample the input distributions and, in
the inner layer, we use the resample distribtuions to drive simulation runs to obtain the output distributions.
The issue with a naive use of the bootstrap in this context is that the total required simulation effort can

Chen, Liu, and Lambe huge (Lam and Qian (2021)), essentially an amount that can wash away the Monte Carlo noises. Moreconcretely, as shown in Lam and Qian (2021), this amount is of larger order than the data size itself.
To remedy this computational demand, we utilize a generalization of Lam and Qian (2021) to devise asubsampling procedure. This procedure aims to exploit the form of the covariance in (8) that scales withthe data size n and turns out to reduce the computational demand from a larger order than n to an orderthat is independent of n.
∗ , . . . , X ∗ be random resamples uniformly drawn with replacement
To explain our procedure, let Xi,1i,si
∗ = (1/s ) si δ ∗ be the corresponding subsample empirical distribution.
from Xi,1 , . . . , Xi,ni and let P̂i,s∑i
j=1 Xi, j
i∗∗ , . . . , P̂∗ ). Using the rationale of the bootstrap, we can approximate the LHS of (8) as
Denote P̂θ = (P̂1,sm,sm
1
follows∗∗
θ nCov∗ (Q(t, P̂θ ), Q(s, P̂θ )) ≈ nCov(Q(t, P̂), Q(s, P̂)) ≈ Cov(G(t), G(s)),

and consequently we obtain∗∗θ Cov∗ (Q(t, P̂θ ), Q(s, P̂θ )) ≈ Cov


G(t) G(s)
√ , √
,
n
n

(9)

where Cov∗ denotes the covariance conditional on the input data. Note that we rescale the covariance inthe LHS of (9) by θ because we are estimating a quantity of order 1/n. When we use the smaller averagedata size s = (1/m) ∑m
i=1 si in our subsampling, the estimated quantity we get is of order 1/s = 1/θ n. The
rescaling factor θ thus helps us retain correct estimation for the original quantity of interest.
To implement the approximation in (9), for b = 1, . . . , B, we uniformly sample with replacement∗b , . . . , X ∗b . Then we construct the subsample empirical distribution
from Xi,1 , . . . , Xi,ni si times to obtain Xi,1
i,si
∗b

∗b

∗b
∗b
∗b = (1/s ) si δP̂i,si ∑ j=1 X ∗b and let P̂θ = (P̂1,s1 , . . . , P̂m,sm ). We use P̂θ to drive R independent simulation runs
ii, j

to obtain R outputs Yr∗b , r = 1, . . . , R and then construct the empirical output distribution as follows∗bQ̂(t, P̂θ ) =1 R∑ I(Yr∗b ≤ t).
R r=1

Then we use the ANOVA-based estimator to estimate the covariance in (9):
θ

1 B
∗b
∗b
∑ (Q̂(t, P̂θ ) − Q̄(t))(Q̂(s, P̂θ ) − Q̄(s))
B − 1 b=1

!
B R
1
∗b
∗b
∗b
∗b
−
∑ ∑ (I(Yr ≤ t) − Q̂(t, P̂θ ))(I(Yr ≤ s) − Q̂(s, P̂θ )) ,BR(R − 1) b=1
r=1
∗b

∗b

where Q̄(·) = (1/B) ∑Bb=1 Q̂(·, P̂θ ) is the sample mean of Q̂(·, P̂θ ). We can show this estimator is unbiased
for estimating the LHS of (9) conditional√ on the input data. Algorithm 1 summarizes the procedure toestimate the covariance matrix of G(t)/ n for t = t1 , . . . ,tk . Importantly, this estimator allows to use lesssimulation effort than the standard bootstrap by choosing s to be much smaller than n, while still retaining
meaningful statistical consistency.
4.2 Confidence Band Construction for the Output Distribution Function
Once we obtain the estimated covariance matrix V returned by Algorithma k-dimensional√ 1, we can simulate√
Gaussian vector (Z1 , . . . , Zk ) ∼ N(0,V ) as the estimator for (G(t1 )/ n, . . . , G(tk )/ n). We can also simudlate (W1 , . . . ,Wk ) = (BB(Q̂(t1 , P̂)), . . . , BB(Q̂(tk , P̂))) as the estimator for (BB(Q(t1 , P)), . . . , BB(Q(tk , P))).
Since the processes G(t) and BB(Q(t, P)) are independent, the simulation procedures of (Z1 , . . . , Zk ) andChen, Liu, and Lam
Algorithm 1 Subsampling to Estimate Covariance Matrix
Inputs: B, R, observed inputs {Xi, j : i = 1, . . . , m, j = 1, . . . , ni }, subsample ratio θ , a grid {tl : l =1, . . . , k}Compute si = ⌊θ ni ⌋, i = 1, . . . , mfor b = 1, . . . , B do
∗b , . . . , X ∗b } uniformly with replacement from theFor each i = 1, . . . , m draw a subsample Xi∗b = {Xi,1i,si∗binput {Xi,1 , . . . , Xi,ni }. Form resample empirical distribution P̂i,si∗b∗b , . . . , P̂∗b )Let P̂θ = (P̂1,sm,sm
1
for r = 1 . . . R do∗b
Simulate Ŷr∗b = f (P̂θ ) from the stochastic simulation model f
end for∗b
Form empirical distribution for Q̂(·, P̂θ ) at points {tl : l = 1, . . . , k}:∗bQ̂(tl , P̂θ ) = R1 ∑Rr=1 I(Yr∗b ≤ tl ), l = 1, . . . , k
end for∗bQ̄(tl ) = B1 ∑Bb=1 Q̂(tl , P̂θ ), where l = 1, . . . , kCompute the estimated covariance matrix V ∈ Rk×k at points {tl : l = 1, . . . , k} from B empirical distributions constructed from the subsampled runs:

Vi, j = θ

1 B
∗b
∗b
∑ (Q̂(ti , P̂θ ) − Q̄(ti ))(Q̂(t j , P̂θ ) − Q̄(t j ))B − 1 b=1

!
B R
1
∗b
∗b
−
∑ ∑ (I(Yr∗b ≤ ti ) − Q̂(ti , P̂θ ))(I(Yr∗b ≤ t j ) − Q̂(t j , P̂θ )BR(R − 1) b=1
r=1
return V(W1 , . . . ,Wk ) should also be independent, conditional on the input data. Now in view of the approximation(7), we can deduce that
Wj dG(t) BB(Q(t, P))√
max Z j + √≈ sup √ +.
(10)
j=1,...,k
n
Re
Re
t∈R
Thus, the quantile of the LHS distribution in (10) can be taken as the approximation of the quantile of the
RHS distribution in (10).
To estimate the quantile of the LHS distribution in (10), we independently repeat the above Rq times
and output the (1 − α)-quantile ofW ji
max Z ij + √, i = 1, . . . , Rq
j=1,...,kRe
as q̃1−α . Based on (4), a statistically valid pair of confidence bands on the discretized time grid is given byL(tl ) = Q̂(tl , P̂) − q̃1−α ,U(tl ) = Q̂(tl , P̂) + q̃1−α ,∀l = 1, . . . , k,eI(Yr ≤ tl ) and Yr , r = 1, . . . , Re are i.i.d. outputs of the simulation model
where Q̂(tl , P̂) = (1/Re ) ∑Rr=1
driven by the empirical input distribution P̂. We summarize our confidence band construction procedurein Algorithm 2.

Chen, Liu, and Lam
Algorithm 2 Construction of Confidence Band for the Output Distribution Function
Inputs: Re , Rq , α, observed inputs {Xi, j : i = 1, . . . , m, j = 1, . . . , ni }, a grid {tl : l = 1, . . . , k}, estimatedcovariance matrix V
for r = 1, . . . , Re do
Yr = f (P̂) from the stochastic simulation model f
end forForm empirical distribution for Q̂(·, P̂) at points {tl : l = 1, . . . , k}:eQ̂(tl , P̂) = R1e ∑Rr=1
I(Yr ≤ tl ), l = 1, . . . , k
for i = 1, . . . , Rq do
Simulate Z i = (Z1i , . . . , Zki ) from the k-dimensional Gaussian random vector N(0,V )Simulate W i = (W1i , . . . ,Wki ) as a time transformed discretized brownian bridge BB(Q̂(tl , P̂)), l = 1, . . . , kWi

ei = max Z ij + √Rj
j=1,...,ke

end forFind q̃1−α = (1 − α)-quantile of e1 , . . . , eRq
return L(·) = Q̂(·, P̂) − q̃1−α and U(·) = Q̂(·, P̂) + q̃1−α
5

NUMERICAL EXPERIMENTS

We conduct experiments on an M/M/1 queue and report our experimental findings. This system has 2 input
distributions P = (P1 , P2 ) where P1 denotes the inter-arrival time distribution and P2 denotes the service
time distribution, for which we have n1 and n2 i.i.d. data points available respectively. We consider two
ground truths, one with arrival rate 0.5 and service rate 1, and another with arrival rate 0.9 and service
rate 1. To satisfy the balanced data condition (2), we choose n2 = 2n1 in both settings and we only reportmini ni for convenience. Our goal is to estimate the entire distribution of the average waiting time of thefirst 10 customers entering the system, more precisely the 97.5%-level confidence band for the distributionof W̄10 , the average waiting time of the first 10 customers. i.e., Q(t, P) = EP [1{W̄10 ≤ t}].
In the experiments, we validate the performance of our procedure and investigate the impact of proceduralparameters: subsample ratio θ , subsampling split between outer loop B and inner loop R when estimating
the covariance matrix by Algorithm 1, and provide practical guidelines for choosing these parameters.
Moreover, we investigate the budget split between Algorithms 1 and 2. Note that Algorithm 1 consists of
Rv = BR simulation model runs for the estimation of the covariance matrix while Algorithm 2 consists ofRe simulation model runs for the estimation of Q̂(·, P̂). Suppose we are given a total budget of N simulation
runs (the total number of times we can run the stochastic simulation model) and we allocate Rv to estimate
the covariance matrix V and Re = N − Rv to compute the CDF estimator Q̂(·, P̂). We investigate the optimal
split between Rv and Re . Here we do not take into account the Rq runs in Algorithm 2 because it does
not require running the simulation model, but rather just the generation of normal vectors and thus theruntime is relatively negligible. The performances of different configurations are measured by the coverage
probability as well as the width of the confidence band.
5.1 Performance Validation and Guidelines for Algorithmic Configuration
We validate and examine the performance of a wide range of selections of θ , B, R. For each example and for input data size ranging from 300 to 2000, we test different combinations of θ , B, R. The
subsample size θ mini ni ∈ {5, 30, 60, 90} and the budge allocation (B, R) ∈ {(10, 100), (12, 80), (16, 60),
(25, 40), (40, 25), (200, 5)} when Rv = 1000 and B = ⌊Rv /R⌋. We also consider the case when we havea smaller budget Rv = 500 and test the corresponding budget allocation (B, R) ∈ {(5, 100), (6, 80), (8, 60),
(12, 40), (20, 25), (100, 5)}.

Chen, Liu, and Lam
To calculate the coverage probability of the estimated confidence band, we first run the simulation modelwith the true input distributions (exponential distributions in the M/M/1 case) 100000 times to approximate
the true output distribution. When the number of simulation runs is large enough, the empirical output
distribution we constructed with the true input distributions is a good enough approximation of the true output
distribution. We then perform 200 independent runs of the confidence band construction procedure andreport the percentage of runs in which the constructed confidence band covers the true output distribution.
We first establish a guideline for choosing the optimal allocation between B and R given a subsample sizeθ mini ni . As shown in Figures 1-3, for different θ mini ni , the optimal R is generally between (2/3)θ mini niand θ mini ni since larger R results in undercoverage while smaller R produces a much wider confidenceband which is less informative and leads to overcoverage. As observed from Figures 1-3, when R is below(2/3)θ mini ni , most of the experiments have 100% coverage probability and thus these confidence bands
are too conservative. To verify that this selection of R is independent of the simulation budget Rv , we also
take a look at Rv = 500 in Figure 4 and observe a similar optimal choice of R. We therefore conclude thatR ∈ [(2/3)θ mini ni , θ mini ni ] is recommended as a practical choice in budget allocation.
0.7

95

Confidence Band Width

Coverage Probability

100

90
85

0.6

0.9/1 min ni = 300
0.5/1 min ni = 300
0.9/1 min ni = 500
0.5/1 min ni = 500
0.9/1 min ni = 1000
0.5/1 min ni = 1000
0.9/1 min ni = 2000
0.5/1 min ni = 2000

0.5
0.4
0.3
0.2

80
5

25

40

R

60

80

0.1

100

5

25

40

R

60

80

100

100

1.0

98

0.8

Confidence Band Width

Coverage Probability

Figure 1: Coverage probability and confidence band width for different R with Rv = 1000 and θ mini ni = 30.

96940.9/1 min ni = 300
0.5/1 min ni = 300
0.9/1 min ni = 500
0.5/1 min ni = 500
0.9/1 min ni = 1000
0.5/1 min ni = 1000
0.9/1 min ni = 2000
0.5/1 min ni = 2000

0.6
0.4
0.2

92
5

25

40

R

60

80

100

5

25

40

R

60

80

100

Figure 2: Coverage probability and confidence band width for different R with Rv = 1000 and θ mini ni = 60.
We then investigate the optimal subsample size. We compare the performance as θ mini ni varies with
optimally tuned R = θ mini ni . In Figure 5, we see that almost all cases have a good coverage probability
close to the nominal level 97.5%. On the other hand, a small subsample size such as 5 gives us a wider
confidence band compared to moderately larger subsample size. Choosing a moderate subsample size
such as 30 gives us a fairly accurate coverage with narrower confidence band. Besides, increasing the
subsample size from 30 does not result in a narrower confidence band. We prefer a small subsample size
because this helps reduce simulation budget consumption. We therefore recommend choosing subsample
size θ mini ni ≈ 30.

Chen, Liu, and Lam
100

1.2
Confidence Band Width

Coverage Probability

99
98
97
96
95
94

1.0

0.9/1 min ni = 300
0.5/1 min ni = 300
0.9/1 min ni = 500
0.5/1 min ni = 500
0.9/1 min ni = 1000
0.5/1 min ni = 1000
0.9/1 min ni = 2000
0.5/1 min ni = 2000

0.8
0.6
0.4
0.2

93
5

25

40

R

60

80

100

5

25

40

R

60

80

100

Figure 3: Coverage probability and confidence band width for different R with Rv = 1000 and θ mini ni = 90.
1.0

98

Confidence Band Width

Coverage Probability

100

96
94
92

0.9/1 min ni = 300
0.5/1 min ni = 300
0.9/1 min ni = 500
0.5/1 min ni = 500
0.9/1 min ni = 1000
0.5/1 min ni = 1000
0.9/1 min ni = 2000
0.5/1 min ni = 2000

0.8
0.6
0.4
0.2

90
5

25

40

R

60

80

100

5

25

40

R

60

80

100

Figure 4: Coverage probability and confidence band width for different R with Rv = 500 and θ mini ni = 60.
5.2 Budget Split between Covariance Matrix Estimation and Distribution Function Construction
We also investigate the optimal split between Rv and Re given a total budget N. Assume that we have total
budget N = 1000 to be allocated to Rv and Re . Based on the practical guidelines for choosing R and θ mini ni
in Section 5.1, we fix the optimal subsample size θ mini ni = 30 and optimal R = θ mini ni = 30 in all casesand only vary Rv . We try six different splits Rv = 90, 180, 300, 510, 690, 900 (accordingly, Re = 1000 − Rv ).
From Figure 6, we see that Rv = 180 or Rv = 300 is generally a good choice of budget assigned to covariancematrix estimation in all cases. Such an allocation gives us the tightest confidence band while having theproper coverage probability. Therefore, we suggest that given a fixed total budget N, we assign 20%-30%
to the estimation of covariance matrix and allocate the rest to estimate the output distribution function.
This practical guideline is also verified when the total budget is increased to N = 1500, from which weobserve that the optimal Rv is around 300-500 (we do not present a figure for N = 1500 due to page limit) .
ACKNOWLEDGMENTS
We gratefully acknowledge support from the National Science Foundation under grants CAREER CMMI1834710 and IIS-1849280.
REFERENCES
Banks, J., J. Carson, B. Nelson, and D. Nicol. 2014. Discrete Event System Simulation. Fifth ed. Harlow: Pearson.
Barton, R., H. Lam, and E. Song. 2022. “Input Uncertainty in Stochastic Simulation”. To appear in the Palgrave Handbookof Operations Research, Chapter 17.
Barton, R. R. 2012. “Tutorial: Input Uncertainty in Output Analysis”. In Proceedings of the 2012 Winter Simulation Conference,edited by C. Laroque, J. Himmelspach, R. Pasupathy, O. Rose, and A. Uhrmacher, 1–12. Piscataway, New Jersey: Instituteof Electrical and Electronics Engineers, Inc.
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.

100

0.300

99

0.275

Confidence Band Width

Coverage Probability

Chen, Liu, and Lam

98
97
96
95

0.9/1 min ni = 300
0.5/1 min ni = 300
0.9/1 min ni = 500
0.5/1 min ni = 500
0.9/1 min ni = 1000
0.5/1 min ni = 1000
0.9/1 min ni = 2000
0.5/1 min ni = 2000

0.250
0.225
0.200
0.175
0.150

94
5

30

*min(ni)

60

90

5

30

*min(ni)

60

90

Figure 5: Coverage probability and confidence band width for different θ mini ni with optimally tuned R.
100

0.325
Confidence Band Width

Coverage Probability

0.300
98
96
94

0.9/1 min ni = 300
0.5/1 min ni = 300
0.9/1 min ni = 500
0.5/1 min ni = 500
0.9/1 min ni = 1000
0.5/1 min ni = 1000
0.9/1 min ni = 2000
0.5/1 min ni = 2000

0.275
0.250
0.225
0.200
0.175
0.150

92

0.125
90

180

300

510
Rv

690

900

90

180

300

510
Rv

690

900

Figure 6: Coverage probability and confidence band width for different splits with N = 1000 and optimallytuned R and θ mini ni .
Barton, R. R., and L. W. Schruben. 1993. “Uniform and Bootstrap Resampling of Empirical Distributions”. In Proceedings ofthe 1993 Winter Simulation Conference, edited by G. W. Evans, M. Mollaghasemi, E. Russell, and W. Biles, 503–508.
Piscataway, New Jersey: Institute of Electrical and Electronics Engineers, Inc.
Barton, R. R., and L. W. Schruben. 2001. “Resampling Methods for Input Modeling”. In Proceedings of the 2001 Winter
Simulation Conference, edited by B. A. Peters, J. S. Smith, D. J. Medeiros, and M. W. Rohrer, 372–378. Piscataway, New
Jersey: Institute of Electrical and Electronics Engineers, Inc.
Cheng, R. C., and W. Holland. 1998. “Two-Point Methods for Assessing Variability in Simulation Output”. Journal of StatisticalComputation Simulation 60(3):183–205.
Cheng, R. C., and W. Holland. 2004. “Calculation of Confidence Intervals for Simulation Output”. ACM Transactions on
Modeling and Computer Simulation (TOMACS) 14(4):344–362.
Cheng, R. C., and W. Holloand. 1997. “Sensitivity of Computer Simulation Experiments to Errors in Input Data”. Journal of
Statistical Computation and Simulation 57(1-4):219–241.
Chick, S. E. 2001. “Input Distribution Selection for Simulation Experiments: Accounting For Input Uncertainty”. Operations
Research 49(5):744–758.
Corlu, C. G., A. Akcay, and W. Xie. 2020. “Stochastic Simulation Under Input Uncertainty: A Review”. Operations Research
Perspectives 7:100162.
Ghosh, S., and H. Lam. 2019. “Robust Analysis in Stochastic Simulation: Computation and Performance Guarantees”. Operations
Research 67(1):232–249.
Henderson, S. G. 2003. “Input Modeling: Input Model Uncertainty: Why Do We Care and What Should We Do About It?”. In
Proceedings of the 2003 Winter Simulation Conference, edited by S. Chick, P. J. Sánchez, D. Ferrin, and D. J. Morrice,
90–100. Piscataway, New Jersey: Institute of Electrical and Electronics Engineers, Inc.
Lam, H. 2016. “Advanced Tutorial: Input Uncertainty and Robust Analysis in Stochastic Simulation”. In Proceedings of the2016 Winter Simulation Conference, edited by T. M. K. Roeder, P. I. Frazier, R. Szechtman, E. Zhou, T. Huschka, andS. E. Chick, 178–192. Piscataway, New Jersey: Institute of Electrical and Electronics Engineers, Inc.
Lam, H., Z. Liu, and H. Qian. 2022. “Efficient Input Uncertainty Quantification Without Nested Simulation: A Subsampled
Infinitesimal Jackknife Approach”. Preprint.

Chen, Liu, and Lam
Lam, H., and H. Qian. 2016. “The Empirical Likelihood Approach to Simulation Input Uncertainty”. In Proceedings of the2016 Winter Simulation Conference, edited by T. M. K. Roeder, P. I. Frazier, R. Szechtman, E. Zhou, T. Huschka, andS. E. Chick, 791–802. Piscataway, New Jersey: Institute of Electrical and Electronics Engineers, Inc.
Lam, H., and H. Qian. 2017. “Optimization-Based Quantification of Simulation Input Uncertainty via Empirical Likelihood”.
arXiv preprint arXiv:1707.05917.
Lam, H., and H. Qian. 2019. “Random Perturbation and Bagging to Quantify Input Uncertainty”. In Proceedings of the 2019
Winter Simulation Conference, edited by N. Mustafee, K.-H. G. Bae, S. Lazarova-Molnar, M. Rabe, and C. Szabo, 320–331.
Piscataway, New Jersey: Institute of Electrical and Electronics Engineers, Inc.
Lam, H., and H. Qian. 2021. “Subsampling to Enhance Efficiency in Input Uncertainty Quantification”. Operations Research.
Law, A. M., W. D. Kelton, and W. D. Kelton. 2007. Simulation Modeling and Analysis, Volume 3. New York: Mcgraw-hill.
Lin, Y., E. Song, and B. L. Nelson. 2015. “Single-Experiment Input Uncertainty”. Journal of Simulation 9(3):249–259.
Nelson, B. 2013. Foundations and Methods of Stochastic Simulation: A First Course. New York: Springer Science & BusinessMedia.
Serfling, R. J. 2009. Approximation Theorems of Mathematical Statistics, Volume 162. New York: John Wiley & Sons.
Song, E., and B. L. Nelson. 2015. “Quickly Assessing Contributions to Input Uncertainty”. IIE Transactions 47(9):893–909.
Song, E., and B. L. Nelson. 2019. “Input-Output Uncertainty Comparisons for Discrete Optimization via Simulation”. Operations
Research 67(2):562–576.
Song, E., B. L. Nelson, and C. D. Pegden. 2014. “Advanced Tutorial: Input Uncertainty Quantification”. In Proceedings ofthe 2014 Winter Simulation Conference, edited by A. Tolk, S. Diallo, I. Ryzhov, L. Yilmaz, S. Buckley, and J. Miller,
162–176. Piscataway, New Jersey: Institute of Electrical and Electronics Engineers, Inc.
Van der Vaart, A. W. 2000. Asymptotic Statistics, Volume 3. Cambridge: Cambridge university press.
Wasserman, L. 2006. All of Nonparametric Statistics. New York: Springer Science & Business Media.
Xie, W., B. L. Nelson, and R. R. Barton. 2014. “A Bayesian Framework for Quantifying Uncertainty in Stochastic Simulation”.
Operations Research 62(6):1439–1452.
Zouaoui, F., and J. R. Wilson. 2003. “Accounting For Parameter Uncertainty in Simulation Input Modeling”. Iie Transactions 35(9):781–792.
Zouaoui, F., and J. R. Wilson. 2004. “Accounting For Input-Model and Input-Parameter Uncertainties in Simulation”. IIE
Transactions 36(11):1135–1151.

AUTHOR BIOGRAPHIESMOTONG CHEN is a Master of Financial Engineering student in the Department of Industrial Engineering and Operations
Research at Columbia University. She received her bachelor’s degree in Computer Science and Philosophy from Columbia
University, and was a rate options trader at Bank of America. Her research interests include Monte Carlo simulation anduncertainty quantification. Her email address is mc3917@columbia.edu.
ZHENYUAN LIU is a Ph.D. student in the Department of Industrial Engineering and Operations Research at Columbia
University. He received his bachelor’s degree from the Department of Financial Mathematics, School of Mathematical Sciences,
Peking University. His research focuses on data-driven stochastic optimization, uncertainty quantification, rare-event estimation
and applied probability. His email address is zl2817@columbia.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 from 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,
uncertainty quantification, risk analysis, and stochastic and robust optimization. His email address is henry.lam@columbia.edu.