Proceedings of the 2016 Winter Simulation Conference
T. M. K. Roeder, P. I. Frazier, R. Szechtman, E. Zhou, T. Huschka, and S. E. Chick, eds.

THE EMPIRICAL LIKELIHOOD APPROACH TO SIMULATION INPUT UNCERTAINTYHenry Lam

Huajie Qian

Department of Industrial and Operations Engineering
University of Michigan
1205 Beal Ave.
Ann Arbor, MI 48109, USA

Department of Mathematics
University of Michigan
2074 East Hall
530 Church Street
Ann Arbor, MI 48109, USA

ABSTRACT
We study the empirical likelihood method in constructing statistically accurate confidence bounds forstochastic simulation under nonparametric input uncertainty. The approach is based on positing a pair of
distributionally robust optimization, with a suitably averaged divergence constraint over the uncertain input
distributions, and calibrated with a χ 2 -quantile to provide asymptotic coverage guarantees. We present thetheory giving rise to the constraint and the calibration. We also analyze the performance of our stochastic
optimization algorithm. We numerically compare our approach with existing standard methods such as thebootstrap.
1

INTRODUCTION

Stochastic simulation often relies on input distributions that are not fully known but only observed via a finite
sample of input data. The quantification of the errors in performance analyses due to the misspecificationof input models is known as input uncertainty. This paper focuses particularly on situations where there areno specific assumptions on the parametric form of the input models, or nonparametric input uncertainty.
One major goal in this context is to construct confidence interval (CI) for the target performance measurethat accounts for input model errors on top of simulation errors.
This paper proposes an optimization-based method to construct such a CI. The upper and lower boundsof the CI are obtained as the optimal values of a pair of optimization programs posited over the space of
input probability distributions. These programs can be interpreted as finding the upper and lower worst-caseperformance measures, subject to a constraint on the unknown input distributions within a neighborhood of
the empirical distributions. We show that, by choosing the right size of this neighborhood and the appropriate
nonparametric distance to measure its size, these worst-case optimizations lead to asymptotically exact
confidence bounds for covering the true performance measure. We develop these statistical guarantees
by using the empirical likelihood (EL) method, which can be viewed as a nonparametric analog of the
maximum likelihood theory.
Our approach resembles, but also should be contrasted with, the literature of distributionally robust
optimization (DRO) (e.g., Delage and Ye 2010, Ben-Tal et al. 2013, Goh and Sim 2010). This literature aims
to find worst-case stochastic performance measures (or optimizes decisions over the worst-case) subject
to so-called uncertainty sets that represent the information or uncertainty about an underlying probability
distribution. One common constraint is a nonparametric neighborhood ball around a baseline model,
typically corresponding to the modeler’s best guess of the truth, measured by some nonparametric distancesuch as φ -divergence (Pardo 2005). When this ball contains the true distribution, the optimal values of theworst-case optimizations will cover the true measure. Thus, it is argued that one should calibrate the ballsize by estimating the divergence between the data and the baseline distribution, and consequently translate
the confidence of this estimation into a confidence bound for the performance measure. In practice, thisLam and Qian
implies calibration using a goodness-of-fit χ 2 -quantile, with the degree of freedom corresponding to theeffective number of categories in a discrete distribution, or divergence estimation that involves estimating
a functional of density. These methods could run into the statistical challenges in divergence estimation,
and also potentially a loss of coverage accuracy in translating the uncertainty set into a CI.
One main contribution of this paper is thus a new way of interpreting these DROs via the EL method.
Through this connection we locate the precise nonparametric distance one should use in constructing valid
CIs under nonparametric input uncertainty. Our distance deviates from those in the literature as it consistsof a weighted average of a collection of divergences, each applied on an independent input model needed in
the simulation. We also show that the size of the nonparametric neighborhood can be taken as a χ 2 -quantile,
with degree of freedom exactly one, independent of the number of input models and the continuity of theserandom variates. This is in sharp contrast with the proposed calibration methods in the DRO literature.
In general, the optimization programs we impose is simulation-based. As another contribution, weinvestigate the properties of mirror descent stochastic approximation (MDSA) (Nemirovski et al. 2009)used to locally solve these programs. This simulation-based optimization framework provides an alternate
approach to quantifying nonparametric input uncertainty compared to the current major technique of bootstrap
resampling. The latter is a sampling approach that repeatedly generates new empirical distributions to drive
the simulation runs. On a high level, our approach trades the computational load in resampling with the
optimization routine needed in the iteration of the SA algorithm. We investigate this tradeoff, and compare
the performances of the obtained CIs and their vulnerability to the simulation noise.
The rest of this paper is organized as follows. Section 2 lays out our optimization programs and explainsthem using the EL method. Section 3 presents and analyzes our MDSA algorithm. Section 4 show somenumerical results and comparison with the bootstrap.
2

EMPIRICAL-LIKELIHOOD-BASED CONFIDENCE INTERVAL

We consider a performance measure in the form
Z(P1 , . . . , Pm ) = EP1 ,...,Pm [h(X1 , . . . , Xm )] ,(1)

where each Xi = (Xi (1), . . . , Xi (Ti )) is a sequence of Ti i.i.d. random variables under an independent input
model (or distribution) Pi , and Ti is a deterministic run length. The function h mapping from XT11 × · · · × XTmm
to R is assumed computable given the inputs Xi ’s.
Our premise is that each Pi is unknown but ni i.i.d. data Xi,1 , . . . , Xi,ni are available. The true value of(1) is therefore unknown even under abundant simulation runs. Our goal is to find a (1 − α) level CI for
the true performance measure.
Our approach is the following. For model i, we consider the probability simplex of the weights
wi = (wi,1 , . . . , wi,ni ) over the support set {Xi,1 , . . . , Xi,ni }. We consider the pair of optimization programs
Lα /Uα := min / max Z(w1 , . . . , wm )msubject toni

2− 2 ∑ ∑ log ni wi, j ≤ X1,1−αi=1 j=1ni(2)∑ wi, j = 1, for all ij=1

wi, j ≥ 0, for all i, j
2where X1,1−αis the 1 − α quantile of the chi-square distribution with degree of freedom one, and each wiin Z(w1 , . . . , wm ) should be interpreted as the distribution putting weight wi, j on Xi, j , for j = 1, . . . , ni .
The formulation (2) can be interpreted as the worst-case optimizations over the m independent inputidistributions subject to a weighted-averaged divergence. Note that the quantity −(1/ni ) ∑nj=1
log ni wi, j isLam and Qian
the Burg-entropy divergence (Ben-Tal et al. 2013) between the probability weights wi and the uniform
weights. Thus, letting n = ∑m
i=1 ni be the total sample size, we have!
m
1 ni
1 m nini− ∑ log ni wi, j− ∑ ∑ log ni wi, j = ∑n i=1 j=1
ni j=1i=1 nwhich can be viewed as a weighted average of the individual Burg-entropy divergences imposed on different
2
input models, each weight being ni /n. The first constraint in (2) imposes a bound X1,1−α/(2n) on thisaveraged divergence.
2.1 Empirical Likelihood Theory for Sum of Means
We justify the formulation (2) using the EL method. First proposed by Owen (1988), this method can beviewed as a nonparametric counterpart of maximum likelihood theory. Analogous to Wilks’ Theorem (Coxand Hinkley 1979) in maximum likelihood theory that states the convergence of the so-called logarithmic
likelihood ratio to a χ 2 -distribution, the nonparametric profile likelihood ratio in EL converges similarly.
This will be the key to obtaining our formulation (2).
We describe the EL method. Given m independent sets of data {Xi,1 , . . . , Xi,ni }, we define the nonparaiwi, j . The
metric likelihood, in terms of the probability weights wi for the i-th input model, to be ∏nj=1nmimulti-sample likelihood is ∏i=1 ∏ j=1 wi, j . It can be shown, by a simple convexity argument, that uniform
niweights wi, j = 1/ni for each model maximize ∏mi=1 (1/ni ) . Moreover, uniform weights still maximizeieven if one allows putting weights outside the support of data, in which case ∑nj=1
wi, j < 1 for some i,
ni
m
m
n
i
making ∏i=1 ∏ j=1 wi, j even smaller. Therefore, ∏i=1 (1/ni ) can be viewed as the nonparametric maximum
likelihood estimate.
To proceed, we need to define a parameter of interest that is determined by the input models. Inference
about this parameter is carried out based on the so-called profile nonparametric likelihood ratio, which is
the maximum likelihood ratio between all weights that empirically produce a given value of the parameter,
and the uniform weights (i.e. the MLE). In our case the parameter of interest is the performance measure
Z(P1 , . . . , Pm ), which is however possibly nonlinear in Pi ’s and hence not easy to work with. Fortunately, itturns out that there isn’t much loss to deal with the special case h(X1 , . . . , Xm ) = ∑mi=1 hi (Xi (1)) for somehi : R → R and Ti = 1 for all i, i.e. the parameter of interest is simply the sum of means of random variables
hi (Xi (1)). We will focus on this case for now. To simplify further, let us consider estimating ∑m
i=1 EXi ,where for convenience we write Xi in place of Xi (1). The profile nonparametric likelihood ratio is defined
as
(
)
m

R(µ) = max

ni

m

ni

ni

∏ ∏ ni wi, j ∑ ∑ wi, j Xi, j = µ, ∑ wi, j = 1, wi, j ≥ 0, for all i, j ,i=1 j=1i=1 j=1(3)

j=1and is defined to be 0 if the optimization problem in (3) is infeasible.
Before getting into details of the asymptotic theory, we note that usually the logarithmic profile
nonparametric likelihood ratio at the true value, i.e. −2 log R(∑m
i=1 EXi (1)), asymptotically follows some chisquare distribution whose degree of freedom is equal to the number of non-probability-simplex constraints,
nie.g. ∑mi=1 ∑ j=1 wi, j Xi, j = µ in (3). As another interpretation, treating the weights wi as parameters, thedegree of freedom is the difference in dimensions of the full and constrained parameter space, and typicallyd constraints result in a loss of d dimensions of the parameter space.
Now we state our theorem.
Theorem 1 Let Xi be a random variable distributed under Pi . Assume 0 < ∑m
i=1 Var(Xi ) < ∞, and mini∈I ni ≥c ∑i∈I ni always holds for some constant c > 0, where I = {i|Var(Xi ) > 0}. Then −2 log R(∑m
i=1 EXi ) converges
in distribution to X12 , the chi-squared distribution with degree of freedom one, as ni → ∞ for i ∈ I.
In this theorem only sample sizes of those having positive variance are required to grow to infinity. To
see the reason for this, note that data of inputs with zero variance are always equal to the true mean so weLam and Qian
can always assign them the uniform weights in (3). Other than these zero-variance inputs, the conditionmini∈I ni ≥ c ∑i∈I ni basically forces the sample sizes to grow at the same rate. Theorem 1 is a multi-sample
generalization of the well known empirical likelihood theorem for single-sample mean, which can be foundin Chapter 2 of Owen (2001). We state it as a special case where m = 1.
Theorem 2 Let Y1 ,Y2 , . . . ,Yn be i.i.d. random variables distributed under some distribution P, 0 < Var(Y1 ) <∞. Then −2 log R(EY1 ) converges in distribution to X12 , as n → ∞. The function R(·) here is the same asthat in (3) with m = 1, n1 = n, X1, j = Y j .
2.2 Empirical-likelihood-based Confidence Interval
We discuss how to construct confidence interval for the quantity Z(P1 , . . . , Pm ) based on the EL theory
presented in the last section, and thus justify the validity of the formulation (2). We will first study thelinear output case, and then discuss the general performance measure in (1).
As pointed out in the last section, linear output takes the form h(X1 , . . . , Xm ) = ∑mi=1 hi (Xi ), where Xi
is distributed under Pi . For output of this particular form, Theorem 1 implies that the optimization pair (2)gives an asymptotically correct confidence interval with coverage probability 1 − α, namely
Theorem 3 Assume 0 < ∑m
i=1 Var(hi (Xi )) < ∞, and mini∈I ni ≥ c ∑i∈I ni for some constant c > 0, whereI = {i|Var(hi (Xi )) > 0}. Then!
m

P Lα ≤ ∑ Ehi (Xi ) ≤ Uα

→ 1 − α, as ni → ∞ for i ∈ I,

i=1

where
Lα /Uα := min / max
subject tomniwi, j hi (Xi, j )i=1 j=1m ni∑∑

2− 2 ∑ ∑ log ni wi, j ≤ X1,1−αi=1 j=1(4)

ni∑ wi, j = 1, for all ij=1

wi, j ≥ 0, for all i, j
2
and X1,1−αis the 1 − α quantile of X12 .
2
Proof. By applying Theorem 1 to hi (Xi, j ), we know that the set {µ ∈ R| − 2 log R(µ) ≤ X1,1−α} containsmthe true sum ∑i=1 Ehi (Xi ) with probability 1 − α asymptotically. Note that this set can be identified as
)
(
m

U=

ni

ni

ni

i=1 j=1

j=1m

2, ∑ wi, j = 1, wi, j ≥ 0 .
∑ ∑ wi, j hi (Xi, j ) − 2 ∑ ∑ log ni wi, j ≤ X1,1−αi=1 j=1

It is obvious that Lα /Uα = min / max{µ|µ ∈ U}. So if the set U is convex, then U = [Lα , Uα ], and we
conclude the theorem. To show convexity, it’s enough to notice that the feasible set of (4) is convex, andthe objective is linear in wi, j .
We now extend our result to the general performance measure (1). We show that, by decomposing
Z(w1 , . . . , wm ) into linear and nonlinear parts via a Taylor-like expansion, applying Theorem 3 to the linear
part, and controlling the magnitude of the nonlinear part, we can construct a CI that has an asymptotic
guarantee similar to Theorem 3. We assume that
Ti
Assumption 1 ∑m
i=1 Var(Gi (Xi )) > 0, where Gi (x) = ∑ j=1 EP1 ,...,Pm [h(X1 , . . . , Xm )|Xi ( j) = x].

Lam and Qian
Assumption 2 Denote the index Ii = (Ii (1), . . . , Ii (Ti )) ∈ {1, 2, . . . , Ti }Ti , and Xi,Ii = (Xi (Ii (1)), . . . , Xi (Ii (Ti ))).
Assume that there exists some positive integer k such that h(X1,I1 , . . . , Xm,Im ) has finite 2k-th moment for
all possible choices of Ii ’s.
It is not difficult to see that Assumption 1 is the counterpart of the condition ∑m
i=1 Var(hi (Xi )) > 0 in
Theorem 3, which is needed in dealing with the linear part of Z(w1 , . . . , wm ), and Assumption 2 is used
for controlling the nonlinear part. We state our result for the general performance measure below.
Theorem 4 If Assumptions 1, 2 hold, and mini ni ≥ c ∑m
i=1 ni always holds for some constant c > 0, then





1
1
≤
Z(P
,
.
.
.
,
P
)
≤
U
+
O
→ 1 − α, as ni → ∞,P Lα + O p
1
m
α
p
k−1
k−1
n k
n k
where n = ∑m
i=1 ni and Lα , Uα are defined as in (2).
This theorem shows that if at least sixth moments (k ≥ 3) of the function h are finite, then the difference√k−1
O p (1/n k ) is asymptotically negligible compared with O p (1/ n), the length of the CI. This holds trivially
if, for instance, h is a bounded function.
3

MIRROR DESCENT STOCHASTIC APPROXIMATIONThis section focuses on solving optimization problem (2), using the mirror descent (MD) algorithm. The
MD algorithm was first proposed by Nemirovsky et al. (1982) for deterministic convex optimization in
general normed space, motivated by the challenge that the standard gradient descent in the Euclidean
space may not generally make sense because the primal space, where the solution lies on, can be different
from its dual, where the gradient is defined on. It resolves the issue by mapping the solution to the dual
space and conducts gradient descent therein. For optimizations in Rn , one does not necessarily need touse MD since the primal and the dual space are the same. However, with a given set of constraints, ajudicious choice of a primal-dual map can result in iteration subroutine that is computationally efficient.
Such observations extend to the setting where gradient can only be simulated. In this case, the algorithm
becomes a constrained SA procedure, and leads to MDSA.
Here we outline the MDSA procedure for the following generic version of problem (2) with η > 0
min

Z(w1 , . . . , wm )subject to− 2 ∑ ∑ log ni wi, j ≤ ηnimi=1 j=1ni(5)∑ wi, j = 1, for all ij=1

wi, j ≥ 0, for all i, jninni
The feasible set here, denoted by A , is convex and satisfies A ⊂ ∏m
i=1 P ⊂ R , where P denotes themnprobability simplex associated with a support set of size ni , and n = ∑i=1 ni . The norm on R is chosen toni
be the 1-norm k(w1 , . . . , wm )k1 = ∑mi=1 ∑ j=1 |wi, j |. First we take the entropy distance generating functionmnω on ∏i=1 P i (used in the so-called entropic descent algorithm (Beck and Teboulle 2003))mni

ω(w1 , . . . , wm ) = ∑ ∑ wi, j log wi, j .
i=1 j=1

Strong convexity of ω is guaranteed by the following proposition, whose proof is similar to that of thespecial case m = 1 in Nemirovski et al. (2009).

Lam and Qian
Proposition 1 ω(w1 , . . . , wm ) is strongly convex with parameter 1/m w.r.t. 1-norm in the relative interiorniof ∏m
i=1 P , i.e.
!2
m
m ni
1
∂ω
(ui, j − wi, j ) +ω(u1 , . . . , um ) − ω(w1 , . . . , wm ) ≥ ∑ ∑∑ kui − wi k1 .
2m i=1i=1 j=1 ∂ wi, j (w1 ,...,wm )Corresponding to ω is the prox-function responsible for projecting the solution back to the primal
decision space
m

ni

∂ω
(ui, j − wi, j )i=1 j=1 ∂ wi, j (w1 ,...,wm )V (u1 , . . . , um ; w1 , . . . , wm ) = ω(u1 , . . . , um ) − ω(w1 , . . . , wm ) − ∑ ∑mni= ∑ ∑ wi, j log
i=1 j=1

wi, j
.
ui, j

In each step of MDSA, the current solution (wk1 , . . . , wkm ) is updated via the prox-mappingd∂Z(wi, j − wki, j ) +V (wk1 , . . . , wkm ; w1 , . . . , wm ),
∂
w
i, j
i, j

prox(wk1 , . . . , wkm ) = argmin(w1 ,...,wm )∈A γk ∑(6)\where γk is the step size, and ∂ Z/∂
wi, j is an estimate for each component of the gradient.
3.1 Gradient Estimation
Here we discuss how to estimate the gradient of the objective Z(w1 , . . . , wm ), which is needed in theprox-mapping (6). It is not hard to see that Z(w1 , . . . , wm ) is a multivariate polynomial in wi, j , and hence
its gradient is also a polynomial. However, a closer examination reveals that the number of terms in thispolynomial grows as fast as nT , where n is the order of the sample size and T the time horizon. Thus
it is generally impossible to compute the gradient by naively summing up all the terms, and simulationis needed. We adopt the proposal by Ghosh and Lam (2015a) and Ghosh and Lam (2015b) of using theddirectional derivative ψi, j = dεZ(w1 , . . . , (1 − ε)wi + εei, j , . . . , wm )|ε=0 instead of the standard derivative∂ Z/∂ wi, j , where ei, j is the j-th coordinate vector of Rni . It’s shown that substituting ψi, j in place of∂ Z/∂ wi, j retains all properties needed in the prox-mapping (6), and appealingly ψi, j is simulable. Here isan extension of their result.
dni
Proposition 2 For any (w1 , . . . , wm ) ∈ ∏m
i=1 P with each wi, j > 0, let ψi, j = dε Z(w1 , . . . , (1 − ε)wi +εei, j , . . . , wm )|ε=0 , and ∂ Z/∂ wi, j be the standard derivative. Then
m

ni

m

ni

∂Z

m

∑ ∑ ψi, j (ui, j − wi, j ) = ∑ ∑ ∂ wi, j (ui, j − wi, j ), for any (u1 , . . . , um ) ∈ ∏ P n

i=1 j=1ii=1i=1 j=1

andψi, j = Ew1 ,...,wm [h(X1 , . . . , Xm )Si, j (Xi )] ,where
1 Xi (t) = Xi, j
Si, j (Xi ) = ∑− Ti .
wi, j
t=1
Ti

This proposition suggests the following unbiased estimator for ψi, jψ̂i, j =1 R∑ h(Xr1 , . . . , Xrm )Si, j (Xri ),
R r=1

(7)

where Xr1 , . . . , Xrm , r = 1, . . . , R are R independent replications of the m input processes generated underweights (w1 , . . . , wm ), and are used simultaneously in all ψ̂i, j ’s.

Lam and Qian
3.2 Computing the Prox-mapping\
We discuss how to solve the minimization problem in prox-mapping (6) with ∂ Z/∂
wi, j replaced by ψ̂i, j .
Let’s work with the generic form
m

ni

min

ξi, j wi, j +V (u1 , . . . , um ; w1 , . . . , wm )i=1 j=1
m ni

subject to− 2 ∑ ∑ log ni wi, j ≤ η

∑∑

(8)

i=1 j=1ni∑ wi, j = 1, for all ij=1

wi, j ≥ 0, for all i, j.
where ξi, j ’s are put in place of the gradient estimators. Seeing that the decision space of (8) has dimension∑m
i=1 ni , which can be significantly larger than the number of functional constraints m + 1, our strategy is to
first solve the Lagrangian dual problem, then use the dual optimal solution to recover the primal optimal
solution.
Consider the LagrangianL(w1 , . . . , wm , λ , β )
" 
m

=∑

i=1

ni

wi, j
ξi, j wi, j + wi, j log
ui, j

∑

j=1

ni


+ λi

!#

∑ wi, j − 1

j=1

m

−β

!

ni

2 ∑ ∑ log ni wi, j + η
i=1 j=1

The dual function defined for λ ∈ Rm , β ≥ 0 is
λ ,β) =g(λL(w1 , . . . , wm , λ , β )
 m

wi, j
− 2β log ni wi, j − ∑ λi − β η= ∑ ∑ min (ξi, j + λi )wi, j + wi, j log
w ≥0
ui, j
i=1
i=1 j=1 i, j mm ni 
w̃i, j= ∑ ∑ (ξi, j + λi )w̃i, j + w̃i, j log− 2β log ni w̃i, j − ∑ λi − β η,ui, j
i=1 j=1i=1min

w1 ,...,wm ≥0
m ni

(9)
(10)

where w̃i, j is the minimizer of the inner optimization in (9). This inner optimization can be solved efficiently
λ , β ) = L(w̃1 , . . . , w̃m , λ , β ) can be easily evaluated. Moreover,using Newton’s method. In other words, g(λit can be shown that w̃i, j is always strictly positive, hence ∂ L/∂ wi, j wi, j =w̃i, j = 0. By this observation andthe chain rule, the first derivatives of the dual function are
ni
m ni
∂g
∂g
= ∑ w̃i, j − 1,
= −2 ∑ ∑ log ni w̃i, j − η.
∂ λi j=1∂β
i=1 j=1(11)w̃

To obtain the second derivatives, note that w̃i, j is determined by ξi, j + λi + log µi,i, jj + 1 − w̃2βi, j = 0. Implicit
differentiation gives−w̃2i, j∂ w̃i, j∂ w̃i, j∂ w̃i, j2w̃i, j=,= 0, for s 6= i,=.
∂ λiw̃i, j + 2β ∂ λs∂βw̃i, j + 2βLam and Qian
So the second derivatives areni−w̃2i, j∂ 2g=,∑∂ λi2 j=1 w̃i, j + 2β

∂ 2g
= 0 for s 6= i,∂ λi ∂ λs

m ni
∂ 2g
−4
=
,
∑
∑
2
∂β
i=1 j=1 w̃i, j + 2βni
2w̃i, j∂ 2g=∑.
∂ λi ∂ βw̃
j=1 i, j + 2β

(12)

λ , β ), its gradient and its Hessian
Therefore we have shown in (10), (11) and (12) that the dual function g(λcan all be efficiently evaluated. This allows us to solve the dual problem
λ ,β)
g(λ

max
λ ,β

(13)

subject to λ = (λ1 , . . . , λm ) ∈ Rm , β ≥ 0efficiently using, e.g., the Interior Point method. The following proposition helps justify our proposal ofsolving the dual (13) instead.
Proposition 3 There is a unique optimal solution (λ1∗ , . . . , λm∗ , β ∗ ) to the dual (13). Moreover, the primal
optimal solution w∗i, j to (8) can be obtained by solving
ξi, j + λi∗ + log

w∗i, j2β ∗+ 1 − ∗ = 0, i = 1, . . . , m, j = 1, . . . , ni .
ui, j
wi, j

The full MDSA algorithm is described in Algorithm 1, which possesses the following convergenceguarantee.
Theorem 5 Suppose there exists a unique optimal solution (w∗1 , . . . , w∗m ) for (5) such that for any (w1 , . . . , wm )in the feasible set
m

ni

∂Z

∑ ∑ ∂ wi, j (w1 , . . . , wm )(wi, j − w∗i, j ) = 0 if and only if wi = w∗i for all i.

(14)

i=1 j=1

Then the sequence wk1 , . . . , wkm generated by Algorithm 1 with step size γk such that
∞

∞

k=1

k=1

∑ γk = ∞, ∑ γk2 < ∞,converges to (w∗1 , . . . , w∗m ) almost surely.
It is easy to verify that the condition (14) must hold if the objective is convex, provided the optimal
solution is unique. Besides convexity, condition (14) also applies to objectives that along any ray startingfrom the optimal solution is a strictly increasing function.
4

NUMERICAL EXPERIMENTIn this section we present some simulation study on the validity of the proposed method, and compare ourmethod with bootstrap resampling (Barton and Schruben 1993, Barton and Schruben 2001). We will showthat our method produces statistically valid CIs that have similar coverage, and in some sense are morestable compared with the bootstrap.
We consider the canonical M/M/1 queue with arrival rate 0.8 and service rate 1. The system is empty
when the first customer comes in. We are interested in the probability that the 20-th customer waits forLam and Qian
Algorithm 1 MDSA for solving (5)Input: a parameter η > 0, an initial feasible solution (w11 , . . . , w1m ), a step size sequence γk , and numberof replications per iteration R.
Iteration: set k = 1
1. Estimate ψ̂i,k j = ψ̂i, j (wk1 , . . . , wkm ), i = 1, . . . , m, j = 1, . . . , ni usingψ̂i,k j =1 R
∑ h(Xr1 , . . . , Xrm )Si, j (Xri )R r=1

where Xr1 , . . . , Xrm , r = 1, . . . , R are R independent replications of the m input processes under distributions(wk1 , . . . , wkm )2. Compute the optimal solution (λ1k+1 , . . . , λmk+1 , β k+1 ) to
λ ,β)
g(λ
subject to λ = (λ1 , . . . , λm ) ∈ Rm , β ≥ 0max

as discussed in Section 3.2, with ξi, j = γk ψ̂i,k j , ui, j = wki, j .
3. Solve the equations
wi,k+12β k+1jkk+1
γk ψ̂i, j + λi + log k + 1 − k+1 = 0
wi, j
wi, j
k+1
for (wk+1
1 , . . . , wm ), then go back to 1.

longer than 2 units of time to get served. To put it in the form of (1), let At be the inter-arrival time between
the t-th and (t + 1)-th customers, St be the service time of the t-th customer, andh(A1 , A2 , . . . , A19 , S1 , S2 , . . . , S19 ) = 1{W20 >2} ,(15)

where the waiting time W20 is calculated via the Lindley recursion
W1 = 0,Wt+1 = max{Wt + St − At , 0}, for t = 1, . . . , 19.

(16)

To test the method, we pretend that both the inter-arrival time distribution and service time distribution are
unknown, but data for both inputs are accessible. Specifically, we generate two samples of size n1 and n2
from exponential distribution with rate 0.8 and 1 respectively, and plug in the data into (2) to compute a
95% CI using Algorithm 1.
The parameters of Algorithm 1 are empirically chosen based on several test runs. In all the following√
experiments, the step size is set to be γk = 1/(2 n1 k), and the number of replications for gradient estimation
R = 30. The algorithm is terminated if the difference measured in 1-norm between the average of the last
50 iterates and that of the last 51 to 100 iterates is less than some small ε > 0, which shall be specifiedin each of the experiments below. The average of the last 50 iterates also serves as the final output of thealgorithm.
Apart from the stochastic nature of the algorithm, another source of uncertainty that affects the CI isthe final evaluation of the objective at the confidence bounds, i.e. computing Lα , Uα . The bounds arecomputed by taking the average of a number of replications, called Re , of (15). Since the length of the CI√
with only input uncertainty is of order 1/ n1 , to make the stochastic error negligible, Re is chosen to bemoderately larger than n1 , n2 . The value of Re for each experiments is also to be specified below.

Lam and Qian
4.1 Accuracy of the EL confidence interval
We draw 100 samples of sizes n1 = n2 = 30, n1 = n2 = 50, n1 = n2 = 100, respectively, and construct one
CI based on each sample. Table 1 shows that the CIs have accurate coverage probabilities among all theconsidered sample sizes, and can be computed in a reasonably short time. Here Re = 1000 is large enoughbecause the sample size is no more than 100. To verify that the algorithm has indeed converged underthese stopping criteria ε, trace plots of components of the weight vector are extracted. See Figure 1 forsome of them.
Table 1: Performance of EL method under various sample sizes. # of replications for final evaluation
Re = 1000 for all, stopping criterion parameter ε are specified in each case.
Sample size
Run time per CI(seconds)
Coverage probability estimate
95% CI for coverage probability# of replications of h per CI# of iterations per CI(min+max)

30, ε = 0.0073
15.1
0.94
[0.893, 0.987]
2.7 × 104837

50, ε = 0.0057
18.1
0.92
[0.867, 0.973]
3.4 × 104
1058

0.024

100, ε = 0.0040
22.7
0.94
[0.893, 0.987]
4.5 × 104
1430

0.026

0.023
0.025
0.022

0.021
0.024

w2,1

w1,1

0.02

0.019

0.023

0.018
0.022
0.017

0.016
0.021
0.015

0.014

0.02
0

100

200

300

400

500

600

700

800

iteration(a) w1,1

0

100

200

300

400

500

600

700

800

iteration(b) w2,1

Figure 1: Trace plot of the first component of the weight vectors for inter-arrival time and service time.

4.2 Comparison with bootstrap resampling
The bootstrap has been widely used to approximate sampling distribution of a statistic from which CIs can beconstructed. One method of constructing CIs using the bootstrap is the percentile bootstrap used in Barton
and Schruben (1993) and Barton and Schruben (2001). Given a pair of samples, A1 , A2 , . . . , An1 for theinter-arrival time and S1 , S2 , . . . , Sn2 for the service time, the percentile bootstrap proceeds as follows. Firstchoose B, the number of bootstrap samples and N, the number of simulation replications for each bootstrap
sample. Then for each b = 1, 2, . . . , B, (1) draw a simple random sample of size n1 with replacementfrom {A1 , . . . , An1 } and a sample of size n2 from {S1 , . . . , Sn2 }, denoted by {A1b , . . . , Anb1 }, and {Sb1 , . . . , Sbn2 }respectively; (2) generate N replications of h with the inter-arrival distribution being uniform over A1b , . . . , Anb1
and the service time distribution being uniform over Sb1 , . . . , Sbn2 , and take the average Zb as the estimate of
Eh. Finally output the 0.025(B + 1)-th and 0.975(B + 1)-th order statistics of {Zb }Bb=1 , Z(b0.025(B+1)c) andZ(b0.975(B+1)c) , as the lower and upper limits of the CI.

Lam and Qian
To compare our method with bootstrap resampling, we study two cases. In the first case (Table 2), we
draw 100 samples of size n1 = n2 = 50, and compute CIs for each of them. In the second case (Table 3),
we draw only one sample of size n1 = n2 = 50, and repeat computing 50 CIs based on a single sample.
To make the comparison fair, we appropriately set the parameters B, N for the bootstrap method, and ε, Re
for the EL method so that the run times are comparable. The result shows that the EL method generates
comparable but slightly narrower CIs than the bootstrap, while still keeps similar coverage probabilities. In
the first case, most of the uncertainty in either method comes from the input data, so the CIs they generate
exhibit similar variation. In the second case, the input data is fixed, and hence the fluctuations of the CIs aresolely due to the resampling and the stochastic noises. Table 3 shows that the EL method gives more stable
CIs than the bootstrap, meaning that the EL is less vulnerable to these noises. This can be attributed to thefact that EL is an optimization-based method and hence does not succumb to the additional resampling
noise introduced in the bootstrap.
Table 2: EL method versus bootstrap, 100 confidence intervals for 100 samples.

Run time per CI(seconds)
Coverage probability estimate
95% CI for coverage probability
Mean CI length
Standard deviation of CI length
Mean lower limit
Standard deviation of lower limit
Mean upper limit
Standard deviation of upper limit# replications of h per CI

Bootstrap
B = 1000, N = 500 B = 500, N = 1000
17.4
17.4
0.90
0.94
[0.841, 0.959][0.893, 0.987]
0.595
0.606
0.094
0.086
0.166
0.155
0.121
0.093
0.761
0.760
0.157
0.142
5
5 × 10
5 × 105

EL method
ε = 0.0057, Re = 500
18.1
0.93
[0.880, 0.980]
0.536
0.080
0.196
0.105
0.733
0.133
3.3 × 104

Table 3: EL method versus bootstrap, 50 replications of confidence intervals for a single sample.
Run time per CI(seconds)
Standard deviation of CI length
Standard deviation of lower limit
Standard deviation of upper limit# replications of h per CI
5

Bootstrap,B = 500, N = 1000
17.2
0.02330.01060.0213
5 × 105

EL,ε = 0.0064, Re = 2 × 105
18.9
0.00530.00380.0030
4.3 × 105

CONCLUSION

We have proposed an optimization-based method to quantify nonparametric input uncertainty for simulation
performance measures. The method computes CIs that account for input errors by positing a pair of
optimizations subject to a weighted average of Burg-entropy divergence constraints on the collection of
empirical input models. We have argued the statistical accuracy of these optimizations via the EL method.
We have also investigated an MDSA algorithm to locally solve the optimizations, and compared its numericalperformances with bootstrap resampling. In future work, we will investigate the generalization of thisapproach to other types of performance measures, and more efficient algorithms to solve the involved classof optimizations.

Lam and Qian
ACKNOWLEDGMENTS
We gratefully acknowledge support from the National Science Foundation under grants CMMI-1400391/
1542020 and CMMI-1436247/1523453.
REFERENCES
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. Evans, M. Mollaghasemi,
E. Russell, and W. Biles, 503–508. ACM.
Barton, R. R., and L. W. Schruben. 2001. “Resampling Methods for Input Modeling”. In Proceedings ofthe 2001 Winter Simulation Conference, edited by B. A. Peters, J. S. Smith, D. J. Medeiros, and M. W.
Rohrer, Volume 1, 372–378. Piscataway, New Jersey: Institute of Electrical and Electronics Engineers,Inc.
Beck, A., and M. Teboulle. 2003. “Mirror Descent and Nonlinear Projected Subgradient Methods forConvex Optimization”. Operations Research Letters 31 (3): 167–175.
Ben-Tal, A., D. Den Hertog, A. De Waegenaere, B. Melenberg, and G. Rennen. 2013. “Robust Solutions of
Optimization Problems Affected by Uncertain Probabilities”. Management Science 59 (2): 341–357.
Cox, D. R., and D. V. Hinkley. 1979. Theoretical statistics. CRC Press.
Delage, E., and Y. Ye. 2010. “Distributionally Robust Optimization Under Moment Uncertainty with
Application to Data-Driven Problems”. Operations Research 58 (3): 595–612.
Ghosh, S., and H. Lam. 2015a. “Computing Worst-Case Input Models in Stochastic Simulation”. arXiv
preprint arXiv:1507.05609.
Ghosh, S., and H. Lam. 2015b. “Mirror Descent Stochastic Approximation for Computing Worst-Case
Stochastic Input Models”. In Proceedings of the 2015 Winter Simulation Conference, edited by S. Jain,
R. R. Creasey, J. Himmelspach, K. P. White, and M. Fu, 425–436. Piscataway, New Jersey: IEEE
Press: Institute of Electrical and Electronics Engineers, Inc.
Goh, J., and M. Sim. 2010. “Distributionally Robust Optimization and its Tractable Approximations”.
Operations Research 58 (4-part-1): 902–917.
Nemirovski, A., A. Juditsky, G. Lan, and A. Shapiro. 2009. “Robust Stochastic Approximation Approachto Stochastic Programming”. SIAM Journal on Optimization 19 (4): 1574–1609.
Nemirovsky, A.-S., D.-B. Yudin, and E.-R. Dawson. 1982. Problem Complexity and Method Efficiency in
Optimization. John Wiley & Sons, Inc., Panstwowe Wydawnictwo Naukowe (PWN).
Owen, A. B. 1988. “Empirical Likelihood Ratio Confidence Intervals for a Single Functional”. Biometrika 75(2): 237–249.
Owen, A. B. 2001. Empirical Likelihood. CRC Press.
Pardo, L. 2005. Statistical Inference Based on Divergence Measures. CRC Press.
AUTHOR BIOGRAPHIESHENRY LAM is an Assistant Professor in the Department of Industrial and Operations Engineering at
the University of Michigan, Ann Arbor. His research focuses on stochastic simulation, risk analysis, andsimulation optimization. His email address is khlam@umich.edu.
HUAJIE QIAN is a Ph.D. student in Applied and Interdisciplinary Mathematics at the University ofMichigan, Ann Arbor. His research interest lies in applied probability and simulation optimization. His
email address is hqian@umich.edu.