Proceedings of the 2019 Winter Simulation Conference
N. Mustafee, K.-H.G. Bae, S. Lazarova-Molnar, M. Rabe, C. Szabo, P. Haas, and Y.-J. Son, eds.

RANDOM PERTURBATION AND BAGGING TO QUANTIFY INPUT UNCERTAINTYHenry Lam
Huajie Qian
Department of Industrial Engineering and Operations Research
Columbia University
500 W. 120th Street
New York, NY 10027

ABSTRACT
We consider the problem of estimating the output variance in simulation analysis that is contributed fromthe statistical errors in fitting the input models, the latter often known as the input uncertainty. Thisvariance contribution can be written in terms of the sensitivity estimate of the output and the variance of
the input distributions or parameters, via the delta method. We study the direct use of this representationin obtaining efficient estimators for the input-contributed variance, by using finite-difference and random
perturbation to approximate the gradient, focusing especially in the nonparametric case. In particular, weanalyze a particular type of random perturbation motivated from resampling that connects to an infinitesimal
jackknife estimator used in bagging. We illustrate the optimal simulation allocation and the simulation
effort complexity of this scheme, and show some supporting numerical results.
1

INTRODUCTION

Simulation analysis consists of the repeated generation of random variates from pre-specified input models
which, fed through the system logic, produces outputs that are used for prediction and decision-making.
When the input models are misspecified or noisily estimated from data, these input errors can propagate tothe outputs and lead to an incorrect conclusion. This calls for a need to quantify the impact of these errors,a problem often known as input uncertainty or extrinsic uncertainty, among other names in the literature.
For general overview, see, e.g., Barton et al. (2002), Henderson (2003), Chick (2006), Barton (2012), Songet al. (2014), Lam (2016), and Nelson (2013) Chapter 7.
A common way in quantifying input uncertainty is to estimate the variance of the output that iscontributed from the input noise. These variances carry information on the amount of variability comingfrom each individual input model or relative to the Monte Carlo noise in running the simulation replications.
Moreover, they can be used as constituent ingredients to construct confidence intervals that account for
both the noises from the inputs and the Monte Carlo runs. Established approaches to estimate these
input-contributed variances (or “input variance” for short) roughly fall into two categories. First is to usebootstrap resampling (e.g., Cheng and Holland 1997; Song and Nelson 2015). This consists of resampling
the input data to construct resampled input models, which are then used to drive simulation replications.
The input variance can be approximated using the bootstrapped variance and a removal of the bias induced
from the Monte Carlo noise via an analysis-of-variance. In the Bayesian contexts, this also leads to related
approaches utilizing posterior sampling (Zouaoui and Wilson 2003; Zouaoui and Wilson 2004). The secondline of approach uses delta-method approximation, which involves estimating the sensitivity or gradientof the outputs with respect to the input distributions or parameters, and combining with the estimation
variance of these input quantities. These approaches include the two-point method (Cheng and Holland
1998) and regression-based approaches (Wieland and Schmeiser (2006); Lin et al. 2015; Song and Nelson

Lam and Qian2019). The first line of methods has been studied under both parametric and nonparametric frameworks,
whereas the second line of methods has thus far focused on the parametric case.
A recurrent challenge in estimating input variance is the high computational cost, and relatedly thecomplexity in the simulation plan, in obtaining accurate estimates. On a high level, this is due to the
convolution of the input and Monte Carlo noises for which, to obtain accurate estimation, one would need todisintegrate these two effects or wash away the Monte Carlo noise by abundant simulation runs. In bootstrap
resampling, this challenge manifests in a nested sampling where one first resamples the input data (capturing
the input noise) and then runs a number of replications per resample (capturing the Monte Carlo noise). Bothlayers need adequate sample sizes and lead to a multiplicatively large overall required simulation size. In
delta-method approximation, the corresponding challenge comes from the gradient estimation. In general,
without structural information, this involves a finite-difference estimator that requires a perturbation size
selection to control both bias and variance, and maintaining a small overall error requires again a large
enough Monte Carlo size.
In this paper, we consider the delta-method approximation to estimate input variance. Our focus is theinvestigation of the required simulation replication size needed for the involved gradient estimation, andin turn consistent and accurate estimation of the variance. We consider especially the nonparametric case,which to our best knowledge has not been studied and also exhibits additional complications compared to theparametric counterpart. In the basic form, our approach to estimate gradient is to use random perturbation,
namely by first generating a random vector in the dimension of the parameters, and suitably projecting toeach dimension to obtain a finite-difference estimate for the derivative with respect to each parameter. On
a high level, this idea follows from the two-point method in Cheng and Holland (1998) and especially the
regression approaches in Wieland and Schmeiser (2006) and Song and Nelson (2019). On the other hand, italso bears some differences with these works. First, by choosing the distribution of the perturbation vector
properly, and using a “score function” multiplier like in simultaneous perturbation or smoothing used inzeroth-order stochastic optimization (e.g., Spall 1992; Nesterov and Spokoiny 2017; Flaxman et al. 2005),
the projection of each dimension in our scheme automatically gives rise to small-bias gradient estimate
thanks to the cancellation of higher-order errors. Second, in the nonparametric case, the input parameter isthe entire input distribution function (potentially an infinite-dimensional object), and correspondingly thegradient is measured via the so-called influence function (Hampel 1974) in nonparametric statistics theory,
which arises as a functional derivative. The length of the random perturbation vector or the “effective”
dimension of the problem then grows with the input data size, and leads to a different analysis and ultimate
simulation replication size requirement than the parametric case. Lastly, in this paper, we focus on aparticular type of random perturbation using a multinomial distribution that is motivated from resampling.
We demonstrate some connections of this scheme with the infinitesimal jackknife estimator for varianceestimation in bagging (Wager et al. 2014). We also show some results on the optimal simulation allocation,and demonstrate how a small overall simulation effort of this scheme can be attained.
2

SETTING AND OVERVIEW OF THE APPROACH

Consider a target performance measure ψ(F1 , . . . , Fm ) that depends on independent input probability distributions F1 , . . . , Fm (As an example, think of them as the interarrival and service times in a queueing
system). Suppose that, given these distributions, the simulator can generate unbiased copies for ψ.
We consider the situation where F1 , . . . , Fm are not known, but instead i.i.d. data for each Fi are available.
(i)(i)
Call them X1 , . . . , Xni where ni is the sample size for model i. Then one can form the empirical distributionniF̂i (·) = (1/ni ) ∑ j=1 δX (i) (·) (where δx (·) is the delta measure at x), and a point estimate of ψ(F1 , . . . , Fm )jfrom R simulation runs is given byψ̂(F̂1 , . . . , F̂m ) =1 R∑ ψ̂r (F̂1 , . . . , F̂m )R r=1

Lam and Qian
where ψ̂r (F̂1 , . . . , F̂m ) is an individual simulation run. The variance of ψ̂(F̂1 , . . . , F̂m ) ism

VarFi (IFi (X (i) )) τ 2
+
ni
R
i=1

Var(ψ̂(F̂1 , . . . , F̂m )) ≈ ∑where each term in the summation is the variance contributed from input model i (i.e., the input variance)
and the last term is the variance contributed from the simulation noise (see Proposition 8 in Glynn andLam 2018, and a parametric version in Theorem 1 in Cheng and Holland 1997). Here VarFi (·) refers to
the variance under Fi and X (i) is a generic variable generated from Fi . IFi (·) is the influence function for
input model i, which is defined as the directional derivativeψ(F1 , · · · , Fi−1 , (1 − ε)Fi + εδx , Fi+1 , . . . , Fm ) − ψ(F1 , · · · , Fm ).
εε→0IFi (x) = IFi (x; F1 , . . . , Fm ) = lim+

In other words, IFi (x) represents the infinitesimal effect in mixing Fi with the delta mass at x, δx . Thisfunction has mean zero under Fi , i.e., EFi [IFi (X (i) ; F1 , . . . , Fm )] = 0. The constant τ 2 = Var(ψ j (F1 , . . . , Fm ))is the variance of one simulation replication, under the true input distributions (or estimated distributions,which give rise to the same variance asymptotically under suitable smoothness conditions).
Our goal is to estimate the input variances VarFi (IFi (X (i) )). First of all, an asymptotically accurate
estimate based on the input data is the empirical influence function variance given by
1 ni(i)∑ IFi (X j ; F̂1 , . . . , F̂m )2 .
ni j=1(1)pThis estimate generally has an error of order O p (1/ mini=1,...,m ni ) in estimating VarFi (IFi (X (i) )). The
function IFi (x; F̂1 , . . . , F̂m ) needs simulation approximation, i.e., we use the estimate
1 ni c (i)∑ IF i (X j ; F̂1 , . . . , F̂m )2
ni j=1(2)where we use a random perturbation approach to estimate IFi . In other words, letting F = (F1 , . . . , Fm ) andF̂ = (F̂1 , . . . , F̂m ), we considerR

c i (X (i) ; F̂) = 1 ∑ ψ̂r (F̂1 , · · · , F̂i−1 , (1 − h)F̂i + hG, F̂i+1 , . . . , F̂m ) S j (G) for all i = 1, . . . , m, j = 1, . . . , ni
IF
j
R r=1
h
(i)

(i)

where G is a discrete distribution supported on the data {X1 , . . . , Xni } that is randomized, and S j (G) for each j
is some “score function” of the perturbation such that E[ψ(F̂1 , · · · , F̂i−1 , (1−h)F̂i +hG, F̂i+1 , . . . , F̂m )S j (G)/h|F̂] ≈(i)
IFi (X j ; F̂) for all i = 1, . . . , m and j = 1, . . . , ni . Note that, given the input data, every realization of G canbe identified with a vector of probability weights (G1 , . . . , Gni ) on the ni observations through the relation(i)
PG ({X j }) = G j for each j = 1, . . . , ni where PG denotes the probability under the distribution G, therefore
for convenience G will be referred to as a random discrete distribution or a random vector on the probabilityisimplex {(G1 , . . . , Gni ) ∈ Rni : ∑nj=1
G j = 1, G j ≥ 0 for all j} interchangeably. To get a meaningful estimate(i)for VarFi (IFi (X )), we would like to control the approximation error of (2) for (1) to be o p (1).
Since one can always estimate each input variance VarFi (IFi (X (i) )) separately by fixing all others at
the empirical distributions, for simplicity we focus on the case of a single input model F (m = 1) with ni.i.d. observations X1 , . . . , Xn and our goal is to estimate the input variance VarF (IF(X)). A convenient
formula for Si (G) can be derived for exchangeable random perturbations, i.e., G such that (G1 , . . . , Gn ) andLam and Qian
(Gπ(1) , . . . , Gπ(n) ) have the same distribution for all permutations π : {1, 2, . . . , n} → {1, 2, . . . , n}. Hence
we have E[G1 ] = 1n and ∑ni=1 (Gi − 1n ) = 0, which implies ∑ni=1 E[(G1 − 1n )(Gi − 1n )] = 0 = Var(G1 ) + (n −11)Cov(G1 , G2 ), i.e., Cov(G1 , G2 ) = − n−1
Var(G1 ). Therefore, we let
Si (G) =n−1
1(Gi − ).
nVar(G1 )n

For sufficiently smooth ψ(·), the perturbed performance measure admits a first order Taylor expansion
b + hEX∼G [IF(X; F)|b F]
b = ψ(F)b + h ∑ni=1 IF(Xi ; F)Gb i (see, e.g., Van der Vaart
ψ((1 − h)Fb + hG) ≈ ψ(F)b2000 Chapter 20) around the empirical distribution F. ThereforeEnb ψ(F) ψ((1 − h)Fb + hG)
b i0 − 1 )Si (G)|F̂
Si (G)|F̂ ≈ E
Si (G) + ∑ IF(Xi0 ; F)(G
h
h
n
i0 =1

=

n−1
n−1
b
b
IF(Xi0 ; F)Cov(GIF(Xi ; F)
i0 , Gi ) +∑nVar(G1 ) i0 6=in=−1
b = IF(Xi ; F)bb + n − 1 IF(Xi ; F)IF(Xi0 ; F)∑n i0 6=inb = 0.
recovering the influence function, where the last equality follows because ∑ni=1 IF(Xi ; F)
3

PERTURBATION MOTIVATED FROM RESAMPLING

In this section we analyze a special yet flexible class of random perturbations that are motivated fromresampling of the input data. Consider the multinomial distribution (s, n1 1) with 1 = (1, 1, . . . , 1) ∈ Rn thatcounts s independent and uniform draws from n objects. Letting si to be the count of the i-th object
drawn (∑ni=1 si = s), we set the i-th component of perturbation to be Gi = si /s, and h to be 1 so that1ψ((1 − h)Fb + hG) = ψ(G). It is straightforward to see that Var(G1 ) = sn
(1 − 1n ), hence the appropriate
1
multiplier Si (G) = sn(Gi − n ) for i = 1, . . . , n. This leads to the following estimator for VarF (IF(X))(assuming infinite simulation replications for now)
σs2 =

1 2
1 n 
Cov∗ (ψ(G), sn(Gi − ))
∑
n i=1
n

(3)

where Cov∗ denotes the covariance with respect to the perturbation G given the input data X1 , . . . , Xn (the dataserve as the support of the distribution G), and the subscript s indicates the dependence on s. Adaptationsof the estimator (3) to the case of multiple input models shall be discussed briefly at the end of this section.
3.1 Connection to Bagging and Infinitesimal Jackknife
Note that such a random perturbation G coincides with the empirical distribution induced by the resample{X1∗ , . . . , Xs∗ } of size s drawn uniformly from the data with replacement. In particular, when s = n it
coincides with the standard bootstrap resample. This connects the multinomial perturbation to bagging or
bootstrap aggregating, a popular ensemble technique (i.e., technique that combines multiple base learners)
in machine learning, and interprets the estimator (3) as a so-called Infinitesimal Jackknife (IJ) varianceestimator (Efron 2014; Wager et al. 2014).
Bagging has been used to reduce prediction error of base learners especially tree-based methods, suchas random forests (Breiman 2001). It possesses variance reduction properties (Bühlmann and Yu 2002) andis capable of smoothing statistical estimators that involve model selection (Efron 2014). Given a training
set for a prediction problem, bagging repeatedly draws bootstrap resamples from the data, trains a baseLam and Qian
learner from each resample, and at the end averages all individual base learners as the bagged learner.
For example, a random forest consists of bagged decision trees. The IJ variance estimator estimates thevariance of the bagged learner, and can be readily computed from the bagging procedure as a by-productwith negligible computation overhead.
b ofThe estimator (3) is connected to the IJ estimator as follows. Consider the “plug-in” estimator ψ(F)∗the performance measure as a base learner. With Fb representing a bootstrap resample, i.e., the perturbation
G when s = n, a bagged estimator of the performance measure isb = E∗ [ψ(Fb∗ )].
ψbag (F)b one can argue that theSupposing that ψ(·) is sufficiently smooth, then, by linearizing ψ(Fb∗ ) around F,b relative to ψ(F)b is of order 1/n, i.e., ψbag (F)b = ψ(F)b + O p (1/n). This observation also
bias of ψbag (F)reveals that the variance of the bagged estimator is close to that of the plug-in estimator in the sensebb + O p (1/n))
Var(ψbag (F))= Var(ψ(F)
b + 2Cov(ψ(F),b O p (1/n)) +Var(O p (1/n))= Var(ψ(F))

b + O 1/n3/2≈ Var(ψ(F))
b O p (1/n)) = O(1/n3/2 ) by Cauchy Schwartz
where O(1/n3/2 ) is the correct order because Cov(ψ(F),2b is in the form of an
inequality and Var(O p (1/n)) = O(1/n ). Note that the bagged estimator ψbag (F)bexpectation over the resamples, therefore is smooth in F. This IJ estimator for its variance is derived fromb The following shows the form of the IJ variance estimator (Efron
applying the delta method on ψbag (F).
2014):b isLemma 1 The IJ estimator of Var(ψbag (F))
nσIJ2 = ∑ Cov∗ (ψ(Fb∗ ), Ni∗ )2
i=1

where Ni∗ counts the number of replications of data point Xi in a bootstrap resample. Moreover, letting δx
be the delta measure at x, we havenCov∗ (ψ(Fb∗ ), Ni∗ ) = lim

1

ε→0+ εb(ψbag ((1 − ε)Fb + εδXi ) − ψbag (F).

To link our random-perturbation variance estimator (3) to the IJ estimator, consider the case s = n. Notethat G = Fb∗ and that (nGi −1) can be identified with Ni∗ , therefore nCov∗ (ψ(Fb∗ ), Ni∗ ) = Cov∗ (ψ(G), n2 (Gi −
1
2
2
b
b
n )) for each i and our estimator σn = nσIJ which reconfirms the relation Var(ψbag (F)) ≈ Var(ψ(F)) ≈VarF (IF(X))/n.
3.2 Estimation Procedure
The estimator (3) assumes a hypothetical infinite number of bootstrap samples. In reality, one needs to use
a finite number of perturbations and simulation size per perturbation. Algorithm 1 shows an estimator thatincorporates these simulation runs for σs2 . By associating the multinomial perturbation with resampling,
one can see that each Fbb = G, each Nib /s corresponds to the i-th component Gi , and nNib − s is the multipliersn(Gi − 1n ) in (3).
Next, we explain briefly why σ̂s2 should possess a central limit theorem, and in particular implies
a consistent estimator for σs2 . Our explanation is based on the observation that σ̂s2 is approximately aLam and Qian
U-statistic (Lee 1990; Serfling 2009; Van der Vaart 2000), a type of symmetric statistic that is constructedfrom aggregating all combinations (without replacement) of the data. To this end, we ignore the error of
the gross average ψ̄¯ in estimating the true mean E∗ [ψ(G)], so that the output of Algorithm 1 can be written
in the form
 b

0
1 n
1
1
b0
(
ψ̄
−
E
[ψ(G)])(
ψ̄
−E[ψ(G)])(nNib − s)(nNib − s) =h(ξb , ξb0 )σ̂s2 ≈∗∗∑∑
B(B − 1) b6=b0
n i=1
B(B − 1) b6∑=b0
where each ξb represents the random variates generated in the resampling and the simulation for the b-th
0
perturbation, and h is an appropriate symmetric kernel that maps the pair ξb , ξb0 to (ψ̄ b − E∗ [ψ(G)])(ψ̄ b −0E∗ [ψ(G)]) 1n ∑ni=1 (nNib − s)(nNib − s). This is precisely a U-statistic estimator for σs2 by noting that02E∗ [h(ξb , ξ√b0 )] = σs for all b 6= b . Since U-statistics admit central limit theorems (CLTs), a CLT of2
2
the form B(σ̂s − σs ) → N (0, ρ 2 ) is expected to hold for σ̂s2 when B is chosen large enough.
Algorithm 1 Multinomial Perturbation
Parameters: B ≥ 2, R, s ≤ nfor b = 1 to B do
Draw a resample {X1b , . . . , Xsb } of size s from {X1 , . . . , Xn }, inducing an empirical distribution Fbb
Compute Nib = number of Xi in {X1b , . . . , Xsb } for all ifor r = 1 to R do
Draw a ψ̂r (Fbb )end forCompute ψ̄ b = R1 ∑Rr=1 ψ̂r (Fbb )end forCompute ψ̄¯ = B1 ∑Bb=1 ψ̄ b
Output
 b

0
1
1 n
b0
¯
¯
(
ψ̄
−
ψ̄
)(
ψ̄
σ̂s2 =−
ψ̄)(nNib − s)(nNib − s)∑∑B(B − 1) b6=b0
n i=1
We note that the subsampling scheme in Algorithm 1 resembles the subbagging proposal in Mentch
and Hooker (2016). In the latter, subsampling is introduced to ensure central limit convergence, viewing
a random forest as an infinite-order U-statistic. This is contrary to our use of subsampling here as a tool
to alleviate the total simulation requirement as the next section shows.
3.3 Optimal Simulation Budget Allocation and Overall Effort
In this section we present theories on the budget allocation and the overall required simulation effort inusing Algorithm 1. We start by showing validity of the variance estimator (3) assuming absence of Monte
Carlo errors, and then examine the simulation effort required to achieve consistency in the presence ofMonte Carlo errors, as well as optimal allocation of a given simulation budget.
We shall frequently use the notations ω(·) and Θ(·). a = ω(b) means a/b → ∞ as n, s → ∞, anda = Θ(b) means that there exist constants c, c > 0 such that cb ≤ a ≤ cb for all n, s.
We assume the following smoothness conditions for ψ:
Assumption 1 (Smoothness at true input model)m Zb = ψ(F) + ∑ψ(F)IF(x; F)d(Fbi − Fi )(x) + ε

i=1

where the remainder satisfies E[ε 2 ] = o(n−1 ), and the influence functions IF satisfies VarF (IF(X; F)) > 0
and EF [IF 4 (X; F)] < ∞.

Lam and Qian
Assumption 2 (Smoothness at empirical input model) Assume the empirical influence function convergesb − IF(X1 ; F))4 ] → 0. Let Fbs∗ be the resample of size s. Assume
to the truth in the sense that E[(IF(X1 ; F)
the remainder in the Taylor expansion
m Z

b +∑
ψ(Fbs∗ ) = ψ(F)b Fbs∗ − Fbi )(x) + ε ∗IF(x; F)d(i=1


satisfies E∗ [(ε ∗ − E∗ [ε ∗ ])4 ] = o p s−3 .
We have the following consistency result of the estimator σs2 :Theorem 1 Under Assumptions 1 and 2, it holds
σs2 → σ 2 in probability as n → ∞, s → ∞where σ 2 := VarF (IF(X)) and s is allowed to grow independently of n.
The important message from Theorem 1 is that the resample size s is allowed to grow at an arbitrary
rate relative to the data size n.
Now we examine the simulation effort needed for Algorithm 1. We first need to figure out themean square error (MSE) of the variance estimate σ̂s2 output from Algorithm 1. To proceed, denote byτ 2 = Var(ψ̂(F)) the variance of a single simulation replication under the true input model. Correspondingly
b and τ̂ ∗2 = Var∗ (ψ̂(Fbs∗ )) be the variances when driven by the empirical input modellet τ̂ 2 = Var(ψ̂(F))
and the resampled model.
Assumption 3 (Convergence of simulation variance) τ̂ 2 → τ 2 in probability as n → ∞, and E∗ [(τ̂ ∗2 − τ̂ 2 )2 ] =o p (1) as both n, s → ∞.
Assumptions 1-3 may seem complicated at first glance. However, they can be readily checked to holdfor performance measures with finite time horizons under mild conditions.
Theorem 2 (Finite horizon measure) If ψ(F) = EF [h(X)], where X = (X(1), . . . , X(T )) consists of T
i.i.d. variables distributed under F, h is some performance function, and1.
2.

T
0 < VarF (∑t=1
EF [h(X)|X(t) = X) < ∞Let I = (I(1), . . . , I(T )) be a sequence of indices such that each 1 ≤ I(t) ≤ T , and denote by
XI = (X(I(1)), . . . , X(I(T ))). It holds maxI EF [|h(XI )|4 ] < ∞then Assumptions 1-3 hold.
The MSE of σ̂s2 can be characterized as:
Lemma 2 Under Assumptions 1-3, if the outer size B and resample size s are chosen to be ω(1), then thevariance estimate σ̂s2 output by Algorithm 1 has the following decompositionq
σ̂s2 = σs2 + E + o p ( E[E 2 ])where the leading error E satisfies E[E ] = 0 andE[E 2 ] =2s2 n σ 2 τ 2 2 4sσ 2  2σ 2 τ 2 
+
+
+
.
B2 s
R
B
s
R

(4)

Note that, given a total simulation budget N = BR, the error (4) takes the form
c2
c1
c3 c4 c5+
+ 2+ +
2B
BN N
B N
where ci , i = 1, . . . , 5 do not depend on B or R. Therefore, in order to minimize the error by adjusting B, R,
it is optimal to set B as large as possible, i.e., B = N. This is summarized below:

Lam and Qian
Theorem 3 (Optimal allocation) Under Assumptions 1-3, as the input data size n → ∞, the total simulationbudget N → ∞ and resample size s → ∞, the optimal B, R that minimize the leading term in the MSE ofσ̂s2 in estimating σs2 are
B∗ = N, R∗ = 1pwhich gives an error σ̂s2 − σs2 = E + o p ( E[E 2 ]) whereE[E 2 ] =
2s2 n σ 2
4sσ 2  2σ 2
2 2
2
+
τ
+
+
τ
.
N2 s
N
s

(5)

Theorem 3 is appealing in terms of the transparency of allocation rule. The allocation essentiallyremoves the nestedness of the simulation procedure by requiring only one simulation run per resample.
This property comes intuitively from the observation that the ψ in the covariance term in the estimator (3)is unchanged by replacing with merely one simulation run ψ̂r .
By forcing the MSE (4) to o(1), we identify the range of B, R, s that produces a statistically valid
variance estimate, giving rise to the overall required simulation effort:
Theorem 4 (Overall simulation effort requirement) Under Assumptions 1-3, if the parameters B, R and sof Algorithm 1 are chosen such that
B = ω n1/2 , BR = ω sn1/2 , s = ω(1)(6)
then σ̂s2 is consistent, i.e., σ̂s2 → σ 2 in probability. In particular, the minimum required total simulationbudget is N := BR = ω n1/2 when the resample size is set to s = ω(1).
The second condition in (6) forces the total simulation effort to be N = ω(sn1/2 ). On the other hand,s can grow at an arbitrary rate according to Theorem 1, thus leading to the minimum required effort
N = ω n1/2 . This shows the advantage, in terms of the required simulation effort, of random perturbationb by (ψ̂((1 − h)Fb + hδXi ) − ψ̂(F))/h.
b
over individualized estimation, i.e., estimating each IF(Xi ; F)Withbindividualized estimation, obtaining an estimate for each IF(Xi ; F) costs at least one simulation run, hence
the required√ simulation effort has to scale at least linearly in n, whereas with random perturbation an effortof order n would suffice.
Finally, we present the optimal resample size s given data size n and budget N to minimize the grosserror of σ̂s2 in estimating σ 2 . Note that σ̂s2 − σ 2 = (σ̂s2 − σs2 ) + (σs2 − σ 2 ). The first error is the Monte
Carlo error with a MSE leading term (5) if single-run allocation is used, and is increasing in s. The second
error however is decreasing in s because a larger s corresponds to a smaller perturbation hence better
approximation to the influence function. The optimal s is then found by balancing the two errors. We
have the following result concerning the optimal resample size (note that the optimal allocation is always
non-nested in this case thanks to Theorem 3):
Theorem 5 (Optimal resample size) Let ψ be the finite horizon performance measure from Theorem 2. For
a given simulation budget N = ω(n1/2 ) and data sizes n, if optimal resample size s and B, R of Algorithm
1 are(s∗ = Θ(N 1/2 n−1/4 )if N ≤ n3/2
Θ(n1/2 ) ≤ s∗ ≤ Θ((Nn−1 ) ∧ n) if N > n3/2
B∗ = N, R∗ = 1
then we have a gross error σ̂s2 − σ 2 = E + o p (N −1/2 n1/4 + n−1/2 ), where the leading term has an MSE√
n 1
2
E[E ] = Θ+ .
NnLam and Qian
We comment that in the case of multiple input models m > 1, one can estimate the individual variance
VarFi (IFi (X (i) )) by perturbing the i-th model only and using Algorithm 1, as mentioned before, arriving at
the following estimator that is similar to (3)VarFi\(IFi (X (i) )) =1 ni 
1 2
Cov∗ (ψ(F̂1 , . . . , F̂i−1 , G, F̂i+1 , . . . , F̂m ), sni (Gi − )) for all i = 1, . . . , m.
∑ni j=1ni

An alternative is to perturb all the models simultaneously, i.e., using ψ(G(1) , . . . , G(m) ) where each G(i) isa random perturbation for the i-th model, and estimate all the m influence functions simultaneously. This
leads to the estimator
VarFi\(IFi (X (i) )) =1 ni 
1 2(i)Cov∗ (ψ(G(1) , . . . , G(m) ), si ni (G j − )) for all i = 1, . . . , m∑
ni j=1ni

where si is the resample size used for the i-th model. This simultaneous perturbation is achieved by drawing
resamples from each of the m input models with possibly different resample sizes.
4

NUMERICAL RESULTS

We consider an M/M/1 queuing system with arrival rate 0.5 and service rate 1. Suppose the system is
empty at time zero. The performance measure of interest is the probability that the waiting time of the 20-th
arrival exceeds 2 units of time, whose true value is approximately 0.182. Specifically, the system has two
input distributions, i.e. the inter-arrival time distribution F1 = Exp(0.5) and the service time distribution
F2 = Exp(1), for which we have n1 and n2 i.i.d. data available respectively. If At is the inter-arrival time
between the t-th and (t + 1)-th arrivals, and St is the service time for the t-th arrival, then the system output
ψ(F1 , F2 ) = EF1 ,F2 [1{W20 > 2}]
where 1{W20 > 2} = 1 if W20 > 2 and 0 otherwise, and the waiting time W20 is calculated by the Lindley
recursion W1 = 0,Wt+1 = max{Wt + St − At , 0} for t = 1, 2, . . . , 19.
We test Algorithm 1 under different input data sizes. For each input data size, 1, 000 95%-levelCIs are constructed for the true performance measure ψ(F), each from an independently generated input
data set. We allocate 1, 000 simulation runs to compute σ̂12 , σ̂22 using Algorithm 1, and another 200
simulation runs driven by the empirical input distributions to compute the point estimator ψ̄(Fb1 , Fb2 ). Letτ̂ 2 be the sample
qvariance of the 200 simulation replications. Then the 95%-level CI is computed asψ̄(Fb1 , Fb2 ) ± 1.96 σ̂ 2 /n1 + σ̂ 2 /n2 + τ̂ 2 /200. Tables 1, 2 and 3 report the coverage probability, mean and12

standard deviation of the CI lengths, the estimated ratio between input and simulation standard deviations,and the number of times that the algorithm outputs a negative variance estimate (the estimate is notguaranteed to be always positive, and reset to 0 if it turns out negative).
Table 1: Algorithm 1 with B = 1, 000, R = 1, s1 = n1 , s2 = n2 .
input data sizes
n1 = 60, n2 = 30
n1 = 200, n2 = 100
n1 = 600, n2 = 300
n1 = 2, 000, n2 = 1, 000

coverageestimate
85.1%
83.7%
89.4%
95.1%

mean CI
length
0.427
0.252
0.212
0.242

std. CI
length
0.179
0.114
0.118
0.168

input se.
sim. se.

# neg var

4.25
2.42
2.04
2.56

6
176
419
479

The results show that, compared to standard bootstrap for which s = n (Table 1), using a relatively
small resample size s (Tables 2 and 3) can help reduce the variability of the variance estimate and generateLam and Qian
Table 2: Algorithm 1 with B = 1, 000, R = 1, s1 = s2 = 30.
input data sizes
n1 = 60, n2 = 30
n1 = 200, n2 = 100
n1 = 600, n2 = 300
n1 = 2, 000, n2 = 1, 000

coverageestimate
84.4%
92.8%
94.9%
95.2%

mean CI
length
0.420
0.264
0.179
0.134

std. CI
length
0.169
0.064
0.031
0.016

input se.
sim. se.

# neg var

4.20
2.30
1.34
0.72

2
2
9
82

√
−1/4 1/2
Table 3: Algorithm 1 with B = 1, 000, R = 1 and si = min( 1, 000ni, ni ), the optimal resample size
suggested by Theorem 5.
input data sizes
n1 = 60, n2 = 30
n1 = 200, n2 = 100
n1 = 600, n2 = 300
n1 = 2, 000, n2 = 1, 000

coverageestimate
87.2%
93.1%
93.5%
94.7%

mean CI
length
0.392
0.255
0.174
0.133

std. CI
length
0.103
0.043
0.016
0.007

input se.
sim. se.

# neg var

3.59
2.11
1.18
0.61

0
0
0
0

statistically more accurate CIs. More negative variance estimates suggests more variability. Table 1 has
much more negative variance estimates than the other two tables, therefore using s = n leads to much larger
variability in the estimate. Another sign of estimation variability is the much larger standard deviationsof the CI lengths in Table 1, especially for large data sizes. This is also reflected through the coverage
probability. CIs in Table 1 tend to undercover the truth, except for the case n1 = 2, 000, n2 = 1, 000 which
however have much wider CIs. This misleading phenomenon is probably because that manually resetting
negative variance estimates to zero makes the estimator upward biased which in turn gives rise to wider
confidence intervals. Actually the mean CI length for the case n1 = 2, 000, n2 = 1, 000 in Table 1 is 0.242,
almost twice of those in Tables 2 and 3.
Based on the trend of the coverage probability as the data sizes grow, our variance estimate for input
uncertainty appears statistically consistent. In both Tables 2 and 3, as the data size increases from 30 to
2, 000, the coverage probability approaches the nominal value 95%. Note that this happens under a fixed
the simulation budget 1, 000 for estimating the input variance, which demonstrates that our algorithm isgenerating statistically valid variance estimate even when the data size is of the same order as the simulationbudget, verifying Theorem 4.
Lastly, compared to the optimal resample size (Table 3) suggested by Theorem 5, the simple strategy
of using s1 = s2 = 30 (Table 2) seems to work equally well in the considered cases. Since the optimal
size from Theorem 5 contains hidden constants, in practice it could be useful to have simple and workablestrategies like this one.
5

CONCLUSION

We have studied a direct use of the delta method in estimating nonparametric input variance in simulation
estimation under input uncertainty. Our approach is based on a random perturbation to obtain a finitedifference estimator that approximates the gradient information needed in the delta method, exhibitedby the influence function in the nonparametric setting. By multiplying with a proper score function as
in simultaneous perturbation and other smoothing methods, the projection of our random perturbation
automatically cancels out higher-order approximation terms and leads to a small bias in estimating thegradient. We have specialized our perturbation distribution choice to a multinomial distribution over the
input data set, which connects our approach to the infinitesimal jackknife estimator used in approximating

Lam and Qian
the variance of bagging. We have presented the consistency, sample size requirement and optimal sampling
allocation strategy for our estimator, and provided numerical examples to support our theoretical findings.
Future work includes the full development and expansion of our analyses to more general classes of random
perturbation estimators and more elaborate numerical comparisons with other approaches in both parametricand nonparametric settings.
ACKNOWLEDGMENTS
We gratefully acknowledge support from the National Science Foundation under grants CAREER CMMI1653339/1834710 and IIS-1849280.
REFERENCES
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., S. E. Chick, R. C. Cheng, S. G. Henderson, A. M. Law, B. W. Schmeiser, L. M. Leemis, L. W. Schruben,
and J. R. Wilson. 2002. “Panel Discussion on Current Issues in Input Modeling”. In Proceedings of the 2002 Winter
Simulation Conference, edited by E. Yücesan, C.-H. Chen, J. L. Snowdon, and J. M. Charnes, 353–369. Piscataway, New
Jersey: Institute of Electrical and Electronics Engineers, Inc.
Breiman, L. 2001. “Random Forests”. Machine Learning 45(1):5–32.
Bühlmann, P., and B. Yu. 2002. “Analyzing Bagging”. The Annals of Statistics 30(4):927–961.
Cheng, R. C., 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.
Cheng, R. C., and W. Holland. 1998. “Two-Point Methods for Assessing Variability in Simulation Output”. Journal of StatisticalComputation Simulation 60(3):183–205.
Chick, S. E. 2006. “Bayesian Ideas and Discrete Event Simulation: Why, What and How”. In Proceedings of the 2006 Winter
Simulation Conference, edited by L. F. Perrone, F. P. Wieland, J. Liu, B. G. Lawson, D. M. Nicol, and R. M. Fujimoto,96–106. Piscataway, New Jersey: Institute of Electrical and Electronics Engineers, Inc.
Efron, B. 2014. “Estimation and Accuracy After Model Selection”. Journal of the American Statistical Association 109(507):991–
1007.
Flaxman, A. D., A. T. Kalai, A. T. Kalai, and H. B. McMahan. 2005. “Online Convex Optimization in the Bandit Setting: Gradient
Descent Without a Gradient”. In Proceedings of the Sixteenth Annual ACM-SIAM Symposium on Discrete Algorithms,
SODA ’05, 385–394. Philadelphia, PA, USA: Society for Industrial and Applied Mathematics.
Glynn, P. W., and H. Lam. 2018. “Constructing Simulation Output Intervals under Input Uncertainty via Data Sectioning”. In
Proceedings of the 2018 Winter Simulation Conference, edited by M. Rabe, A. A. Juan, N. Mustafee, A. Skoogh, S. Jain,
and B. Johansson, 1551–1562. Piscataway, New Jersey: Institute of Electrical and Electronics Engineers, Inc.
Hampel, F. R. 1974. “The Influence Curve and Its Role in Robust Estimation”. Journal of the American Statistical Association 69(346):383–393.
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 Jersy: Institute of Electrical and Electronics Engineers, Inc.
Lee, A. J. 1990. U-Statistics: Theory and Practice. New York: Taylor & Francis Group.
Lin, Y., E. Song, and B. Nelson. 2015. “Single-Experiment Input Uncertainty”. Journal of Simulation 9(3):249–259.
Mentch, L., and G. Hooker. 2016. “Quantifying Uncertainty in Random Forests via Confidence Intervals and Hypothesis Tests”.
The Journal of Machine Learning Research 17(1):841–881.
Nelson, B. 2013. Foundations and Methods of Stochastic Simulation: A First Course. New York: Springer Science & BusinessMedia.
Nesterov, Y., and V. Spokoiny. 2017. “Random Gradient-Free Minimization of Convex Functions”. Foundations of Computational
Mathematics 17(2):527–566.
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.

Lam and Qian
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. Y. Diallo, I. O. Ryzhov, L. Yilmaz, S. Buckley, and J. A.
Miller, 162–176. Piscataway, New Jersey: Institute of Electrical and Electronics Engineers, Inc.
Spall, J. C. 1992. “Multivariate Stochastic Approximation Using a Simultaneous Perturbation Gradient Approximation”. IEEE
Transactions on Automatic Control 37(3):332–341.
Van der Vaart, A. W. 2000. Asymptotic Statistics, Volume 3. Cambridge, UK: Cambridge University Press.
Wager, S., T. Hastie, and B. Efron. 2014. “Confidence Intervals for Random Forests: the Jackknife and the Infinitesimal
Jackknife.”. Journal of Machine Learning Research 15(1):1625–1651.
Wieland, J. R., and B. W. Schmeiser. 2006. “Stochastic Gradient Estimation Using a Single Design Point”. In Proceedings ofthe 2006 Winter Simulation Conference, edited by L. F. Perrone, F. P. Wieland, J. Liu, B. G. Lawson, D. M. Nicol, andR. M. Fujimoto, 390–397. Piscataway, New Jersey: Institute of Electrical and Electronics Engineers, Inc.
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 BIOGRAPHIESHENRY LAM is an Associate Professor in the Department of Industrial Engineering and Operations Research at Columbia
University. His research focuses on Monte Carlo simulation, uncertainty quantification, risk analysis, and stochastic and robust
optimization. His email address is khl2114@columbia.edu.
HUAJIE QIAN is a Ph.D. student in the Department of Industrial Engineering and Operations Research at Columbia University.
His research interest lies in simulation uncertainty quantification, data-driven simulation analysis, and stochastic optimization.
His email address is h.qian@columbia.edu.