Proceedings of the 2015 Winter Simulation Conference
L. Yilmaz, W. K V. Chan, I. Moon, T. M. K. Roeder, C. Macal, and M. D. Rossetti, eds.

QUANTIFYING UNCERTAINTY IN SAMPLE AVERAGE APPROXIMATIONHenry Lam

Enlu Zhou

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

School of Industrial & Systems Engineering
Georgia Institute of Technology
755 Ferst Drive NW
Atlanta, GA 30332, USA

ABSTRACT
We consider stochastic optimization problems in which the input probability distribution is not fully known,and can only be observed through data. Common procedures handle such problems by optimizing anempirical counterpart, namely via using an empirical distribution of the input. The optimal solutions
obtained through such procedures are hence subject to uncertainty of the data. In this paper, we explore
techniques to quantify this uncertainty that have potentially good finite-sample performance. We consider
three approaches: the empirical likelihood method, nonparametric Bayesian approach, and the bootstrap
approach. They are designed to approximate the confidence intervals or posterior distributions of the
optimal values or the optimality gaps. We present computational procedures for each of the approaches
and discuss their relative benefits. A numerical example on conditional value-at-risk is used to demonstratethese methods.
1

INTRODUCTION

We are interested in stochastic optimization problems in the form
max E[h(X; θ )],
θ ∈Θ

(1)

where θ ∈ Θ ∈ R p is the decision variable and X ∈ Rd is a random variable. In many applications, theunderlying probability distribution that controls the expectation E[·] is not fully known and can only beaccessed via limited data or Monte Carlo samples. Then it is customary to work on an empirical counterpartof the problem, namely by solving
1 n
max ∑ h(Xi ; θ ),(2)
θ ∈Θ n i=1
where X1 , . . . , Xn are the data. This is well known as sample average approximation (SAA) (Shapiro,Dentcheva, and Ruszczyński 2014) in the stochastic programming literature. Obviously, (2) is subject toapproximation errors with respect to (1), and this error is the main focus of this paper.
Our premise is that beyond the n observed samples, new samples are not easily accessible, either
because of lack of data or because of limited computational capacity in running further Monte Carlosimulation. In this setup, previous results on quantifying the approximation errors have fallen in two
directions. The first is central limit convergence of the optimality gap, i.e. difference between (1) and (2),typically to some Gaussian random variable under suitable scaling (Shapiro, Dentcheva, and Ruszczyński
2014; Kim, Pasupathy, and Henderson 2015). Such results can be used to construct confidence interval(CI) that has asymptotically correct coverage probability as sample size increases. The second approachconsists of sampling complexity results, in the form of concentration inequalities on the optimality gap.
For instance, Shapiro and Nemirovski (2005) present large deviations inequalities that depend on certain
size measures of Θ and the Lipschitz constant of h. In classification-type learning problems, where h isLam and Zhou
often taken as an indicator function, deviation bounds have been obtained based on the entropy or theso-called Vapnik-Chervonenkis dimension of h (e.g. Mendelson 2003).
The above two approaches are both valuable in evaluating the reliability of an adopted solution or
choosing a good sample size n in advance. However, they have antagonistic limitations: Central limit
convergence typically only works well for the large-sample case and has no guarantee for smaller sample
size, while sampling complexity results are typically loose because of their worst-case nature, even though
they hold for any sample size. Driven by these limitations, we are interested in exploring techniques thatcan potentially yield improvement on both ends, by giving tight quantification of approximation errors
even in situations of limited samples.
We consider three distinct approaches: the empirical likelihood (EL) method, the Dirichlet method,and nonparametric bootstrapping. These approaches differ in statistical philosophy and have distint benefits
and disadvantages in terms of both the uncertainty assessment scope and computational load. To give some
initial highlight, the EL method can be viewed as a nonparametric analog of classical likelihood inference,and can generate CI for the true optimal value or the optimality gap in the frequentist sense. On the other
hand, the Dirichlet method is utilized from a Bayesian perspective to approximate a corresponding posterior
distribution, whereas nonparametric bootstrapping can be viewed as a limit of the Dirichlet method bytaking suitable asymptotic on the prior distribution.
To demonstrate our proposed approaches, we consider conditional value-at-risk (CVaR) as a specificexample of (1). CVaR is widely used in portfolio risk management in finance and insurance, among otherapplications (Rockafellar and Uryasev 2000). The distributions of these portfolios are unknown and can
only be represented via data in any real-life contexts. We shall provide some numerics on quantifying theuncertainty of the empirical optimization of CVaR.
Our paper is organized as follows. In Section 2, we first present the three approaches for quantifying
uncertainty of the optimal value output by (2). We then adapt these approaches for optimality gap in Section
3. We discuss some comparisons among the approaches in Section 4, followed by some numerical resultson CVaR in Section 5.
2

METHODS FOR QUANTIFYING UNCERTAINTY OF THE OPTIMAL VALUE

2.1 Empirical Likelihood Method
The EL method is a nonparametric analog of maximum likelihood estimation first proposed by Owen
(1988). To introduce the method, let us first fix some notation. Given the set of data X1 , X2 , . . . , Xn , wedefine a probability simplex over {X1 , . . . , Xn }, denoted w = (w1 , . . . , wn ) where ∑ni=1 wi = 1 and wi ≥ 0 for
2 as the 1 − β quantile of a χ 2 distribution with degree of freedom q. The EL method
all i. We denote χq,βentails that the optimizations
max / minw maxθ ∈Θ ∑ni=1 wi h(Xi ; θ )2subject to−2 ∑ni=1 log(nwi ) ≤ χ p+1,β∑ni=1 wi = 1
wi ≥ 0 for all i(3)

where max / min denote a pair of maximization and minimization, contain the true value of maxθ ∈Θ E[h(X; θ )]with probability at least 1 − β asymptotically.
To get some intuition, the quantity − ∑ni=1 log(nwi ) can be interpreted as a statistical distance between
two probability measures P and Q on the support {X1 , . . . , Xn }, defined by the weights w and the uniform
weights (1/n, . . . , 1/n) respectively. In fact, this distance belongs to the class of so-called φ -divergence,
where the φ function can be identified as φ (x) = − log x here (Pardo 2005). The optimizations in (3) cantherefore be viewed from a distributionally robust optimization perspective: when the underlying probability
distribution of a stochastic problem is not fully known, one can impose an uncertainty set in which theprobability distribution is believed to lie, and calculate the worst-case scenario among all distributions

Lam and Zhou
within the set. In our case, (3) is computing the worst-case scenarios as the max or min for the quantitymaxθ ∈Θ E[h(X; θ )], where the uncertainty set can be thought as the neighborhood around the uniform
distribution on {X1 , . . . , Xn } measured by a φ -divergence with φ (x) = − log x.
2.1.1 Computational LoadWe shall discuss how to compute (3), which consists of a min-max problem and a max-max problem. Notethat the outer objective function maxθ ∈Θ ∑ni=1 wi h(Xi ; θ ) is a convex function in w. Suppose that h(X; ·) is
concave in θ . Then evaluating maxθ ∈Θ ∑ni=1 wi h(Xi ; θ ) for fixed w boils down to a convex optimization, andhence both the outer and inner optimizations in the min-max problem are convex and efficiently solvable.
Moreover, to speed up computation, by observing that the feasible region for w is compact one can consider
an exchange of min-max to max-min via Sion’s minimax theorem (Sion et al. 1958), which givesnn

min max ∑ wi h(Xi ; θ ) = max min ∑ wi h(Xi ; θ )w∈A θ ∈Θ i=1(4)θ ∈Θ w∈A i=1

where A denotes the feasible region in (3). Fixing θ , the inner minimization minw∈A ∑ni=1 wi h(Xi ; θ ) in
the maximin formulation can be solved by Lagrangian relaxation
!
!
n

n

2
max min ∑ wi h(Xi ; θ ) + ν −2 ∑ log(nwi ) − χ p+1,βν≥0,λ ∈R w≥0 i=1

n

−λ

i=1

∑ wi − 1

i=1

n

=

max

2
+λ.
∑ min{wi (h(Xi ; θ ) − λ ) − 2ν log(nwi )} − 2νn log n − ν χ p+1,βν≥0,λ ∈R i=1 wi ≥0(5)Straightforward calculation reveals that(min {wi (h(Xi ; θ ) − λ ) − 2ν log(nwi )} =

wi ≥0




2ν 1 − log h(Xi2ν;θ )−λif h(Xi ; θ ) > λ−∞otherwisewith the optimal solution of wi being 2ν/(h(Xi ; θ ) − λ ) or ∞ correspondingly. This implies that (5) isequal to


n
2ν
2
max
− ∑ 2ν log+ 2νn(1 − log n) − ν χ p+1,β+λh(Xi ; θ ) − λν≥0, λ <mini h(Xi ;θ ) i=1
where the optimal solution of wi is 2ν/(h(Xi ; θ ) − λ ). Hence the inner minimization of the right hand sideof (4) becomes a convex maximization problem with two variables (instead of having n variables in theoriginal formulation).
On the other hand, the max-max problem in (3) can be more challenging, because the outer optimization
involves maximizing the convex function maxθ ∈Θ ∑ni=1 wi h(Xi ; θ ) over w. A quick heuristic for obtaining
a solution for the max-max problem is to do alternating maximization, namely repeatedly fixing w andmaximizing over θ , and fixing θ and maximizing over w, until no improvement is observed. Csiszar andTusnady (1984) contains global convergence results for this sort of procedure under some improvementbound assumptions uniformly on the alternating steps. Similar as above, when θ is fixed, one can writemaxw∈A ∑ni=1 wi h(Xi ; θ ) via Lagrangian relaxation as


n
2ν
2
min
2ν ∑ log+ 2νn(log n − 1) + ν χ p+1,β
+λ
λ
−
h(X
;
θ
)
ν≥0, λ >maxi h(Xi ;θ )ii=1with the optimal solution of wi being 2ν/(λ − h(Xi ; θ )).

Lam and Zhou
2.1.2 Statistical Guarantees
2
can be made rigorous
The statistical guarantee for (3) and the choice of the “neighborhood size” χ p+1,βthrough the theory of the EL method. The key of the method is an analog of the celebrated Wilks’ Theorem(Cox and Hinkley 1979) in parametric likelihood inference, namely that the ratio between the maximum
likelihood and the true likelihood, i.e. the so-called likelihood ratio, converges to a chi-square distribution
in a suitable logarithmic scale.
To apply the EL method on the data set {X1 , . . . , Xn }, one would first define an empirical likelihood∏ni=1 wi . Treating w as the “parameters”, it is not difficult to see that the maximum value of ∏ni=1 wi ,among all w in the probability simplex, is ∏ni=1 (1/n), which can be viewed as a nonparametric analog
of maximum likelihood. The next step is to define a target parameter of interest, i.e. the quantity whose
statistical uncertainty is to be assessed. In our case this is maxθ ∈Θ E[h(X; θ )]. The EL method bridges to
Wilks’ Theorem via a quantity known as the profile likelihood, defined as the maximum ratio between theempirical likelihood and the nonparametric maximum likelihood, optimized over all probability weights withsupport over the data set and satisfying some given moment conditions incorporating the target parameterof interest. The crux is that profile likelihood satisfies similar asymptotic properties as the likelihood ratio
in the parametric context.
We state the above mathematically. We first make the following assumptions:
Assumption 1
1. E[h(X; θ )] is differentiable in θ with ∇θ E[h(X; θ )] = E[∇θ h(X; θ )] for θ ∈ Θ.
2. ∇θ E[h(X; θ0 )] = 0 if and only if θ0 ∈ argmaxθ ∈Θ E[h(X; θ )], and ∑ni=1 wi ∇θ h(Xi ; θ0 ) = 0 if and only
if θ0 ∈ argmaxθ ∈Θ ∑ni=1 wi h(Xi ; θ ) for any Xi and w = (w1 , . . . , wn ) is a probability vector.
3. There exists a θ0 ∈ argmaxθ ∈Θ E[h(X; θ )] such that the covariance matrix of the random vector(∇θ h(X; θ0 ), h(X; θ0 )) has rank p + 1 > 0.
Note that there is no assumption on the uniqueness of θ0 . Condition 2 is a first order condition thatis satisfied if, for instance, h(X; ·) is a coersive concave function for any X a.s. and Θ = R p . Condition
3 states that all partial derivatives of h(X; θ0 ) and also h(X; θ0 ) itself are linearly independent. Then we
have:
Theorem 1 Let the maximum and minimum values of the programs (3) be Z and Z respectively. Then,
under Assumption 1, we have
lim inf P max E[h(X; θ )] ∈ [Z, Z] ≥ 1 − β .
n→∞

θ ∈Θ

We shall prove Theorem 1. We denote θ0 as an element in argmaxθ ∈Θ E[h(X; θ )] such that thecovariance matrix of the random vector (∇θ h(X; θ0 ), h(X; θ0 )) has rank p + 1 > 0. The existence of such
θ0 is guaranteed in Condition 3 in Assumption 1. By Conditions 1 and 2 in Assumption 1, θ0 satisfiesE[∇θ h(X; θ0 )] = 0, where 0 is a length p vector of zeros.
We define a profile likelihood as()
R(θ , z) = max

n

n

n

n

∏ nwi : ∑ wi ∇θ h(Xi ; θ ) = 0, ∑ wi h(Xi ; θ ) = z, ∑ wi = 1, wi ≥ 0 for all i .
i=1

i=1

i=1

(6)

i=1

This is the maximum ratio between the empirical likelihood and the nonparametric maximum likelihood,among all probability weights that respect the condition E[∇θ h(X; θ0 )] = 0 and E[h(X; θ0 )] = z. The
parameter z is our target. We let z0 be the true optimal value of maxθ ∈Θ E[h(X; θ )], which is equal toE[h(X; θ0 )]. The following result is an immediate application from Owen (2001):
2
Theorem 2 Under Assumption 1, the profile likelihood (6) satisfies −2 log R(θ0 , z0 ) ⇒ χ p+1
as n → ∞,22where χq denotes the χ -distribution with degree of freedom q.

Lam and Zhou
2 ) → 1−β . So the set C = {(θ , z) : −2 log R(θ , z) ≤In other words, we know that P(−2 log R(θ0 , z0 ) ≤ χq,β2
} forms an asymptotically valid 1−β level confidence region for (θ0 , z0 ), i.e. P((θ0 , z0 ) ∈ C ) → 1−β .
χ p+1,βThis implies that with probability 1 − β asymptotically, there exists a probability vector w such that2
. This further implies that with
∑ni=1 wi ∇θ h(Xi ; θ0 ) = 0, ∑ni=1 wi h(Xi ; θ0 ) = z0 and −2 ∑ni=1 log(nwi ) ≤ χ p+1,βprobability 1 − β asymptotically, z0 is contained by the pair of optimizations

max / minw ∑ni=1 wi h(Xi ; θ0 )subject to∑ni=1 wi ∇θ h(Xi ; θ0 ) = 0
2−2 ∑ni=1 log(nwi ) ≤ χ p+1,βn∑i=1 wi = 1
wi ≥ 0 for all ior equivalently
max / minw maxθ ∈Θ ∑ni=1 wi h(Xi ; θ )subject toθ0 is an optimal solution for maxθ ∈Θ ∑ni=1 wi h(Xi ; θ )2−2 ∑ni=1 log(nwi ) ≤ χ p+1,β∑ni=1 wi = 1
wi ≥ 0 for all i(7)

by using Condition 2 in Assumption 1. By relaxing the first constraint in (7), we conclude that withprobability at least 1 − β asymptotically, z0 is contained in (3), which proves Theorem 1.
Theorem 2 is an asymptotic statement. For small samples, accuracy (measured in terms of the empirical
coverage probability being correctly larger than 1 − β ) can be enhanced by other calibration methods than2
using χ p+1,βin the optimization (3). Chapters 2 and 13 in Owen (2001) provide some discussion.
2.2 Dirichlet Methods
The empirical likelihood method is a frequentist approach based on asymptotics. In contrast, the Bayesian
approach takes a different perspective by characterizing the likelihood of possible alternatives based on
the available data and the chosen prior. It encodes all the information in the posterior distribution anddoes not rely on asymptotics. In this section, we develop a Bayesian, non-parametric method to quantify
the uncertainty in the optimal value of (1). To have a tractable posterior of distribution over a general
support, we will use the Dirichlet process, which was first introduced by Ferguson (1973) and is currently
one of the most popular Bayesian nonparametric models. To ease the explanation, we will first illustrate
the approach on the simple case when the input distribution has a finite support.
2.2.1 Input Distribution with Finite Support
Suppose the true input distribution of X lives on a finite support {x1 , . . . , xn }. The data we have about
the input distribution are Ni ’s — the number of times that we observe each xi . Please note that the finite
support is not merely a simplification; such scenario does happen, for example, when we want to estimate
the distribution of the customer demand for a certain product based on observations of sales (assuming
no lost sales), where the sales quantity is an integer ranging from 0 to the stock-out level. We denote by
F = [p1 , . . . , pn ] an empirical distribution on the support {x1 , . . . , xn }. The Bayesian approach views F as arandom variable, and yields a belief (posterior distribution) of F based on a chosen prior and the likelihoodof observing the data ψ , [N1 , . . . , Nn ]:l(ψ|F) = Πni=1 pNi i .
If we have no prior knowledge of the true distribution, we can assume a uniform prior, which is equivalent
to a Dirichlet distribution with the parameter being a unit vector, denoted as Dir(e) where e is the unit

Lam and Zhou
vector. We can also use an informative prior by choosing an appropriate parameter value in the Dirichletprior. Under the noninformative prior, the posterior distribution is the Dirichlet distribution with parametere + ψ, i.e., Dir(e + ψ) with the probability density functionf (p; e + ψ) =1
Πn pNi ,B(e + ψ) i=1 iwhere the normalization constant B(·) is the multinomial Beta function. Hence, the Dirichlet distribution
Dir(e + ψ) is essentially a distribution of the distribution F given the data ψ.
The uncertainty in the optimal value of (1) can be characterized by the posterior distribution on the
optimal value given the data, i.e., P(maxθ ∈Θ E[h(X; θ )]|ψ). Hence, we need to propagate the posterior on
the distribution of X to the optimal function value, and that can be done by the following two steps:
Sampling: draw F j = [p1j , . . . , pnj ], j = 1, . . . , M from Dir(e + ψ).
Optimization: for each j = 1, . . . , M, solve the optimization problem
(
)

1.
2.

n

max EF j [h(X; θ )] = ∑ pij h(xi ; θ ) ,
θ ∈Θ

i=1

and denote the optimal value by h∗ j .
The procedure above will output M sample optimal values {h∗1 , . . . , h∗M }, which form an empirical
distribution that approximates the posterior distribution of the optimal value h∗ of (1). From this empirical
distribution further information can be extracted. Specifically, we can sort the sample optimal values fromββ
the smallest to the largest such that h∗(1) ≤ . . . ≤ h∗(M) , and take [h∗(dM 2 e) , h∗(dM(1− 2 )e) ] as an estimate ofthe (1 − β ) Bayesian confidence interval (also called “credible interval”) for the optimal value h∗ . The
procedure above can be made more efficient by first evaluating and storing the function values on the
support points, i.e., h(x1 ; θ ), . . . , h(xn ; θ ), and then simply taking a linear combination of these stored values
to get the objective function EF j [h(X; θ )] according to the probabilities of each generated F j .
ββ
The following theorem shows that the estimated credible interval [h∗(dM 2 e) , h∗(dM(1− 2 )e) ] is asymptoticallythe true credible interval as the sample size M goes to infinity. First note that the posterior distribution of
the optimal value h∗ given data ψ isPh∗ (t|ψ) , Pr{max E[h(X; θ )] ≤ t|ψ}.
θ ∈Θ

Hence, the true γ-quantile of this posterior distribution is h∗γ (ψ) = inf{h : Ph∗ (t|ψ) ≥ γ}. Then we have thefollowing theorem.
Theorem 3 Assume Ph∗ (t|ψ) is a continuous distribution. Given the data ψ,(i) as M → ∞, the empirical distribution formed by {h∗1 , . . . , h∗M } provides a uniformly consistent estimatorof the posterior distribution of h∗ , i.e., Ph∗ (t|ψ);(ii) for any γ ∈ (0, 1), limM→∞ h∗(dMγe) = h∗γ (ψ) almost surely.
The proof is similar to that of Theorem 1 in Xie, Nelson, and Barton (2015). Note that {h∗1 , . . . , h∗M }is an i.i.d. sample from the posterior distribution Ph∗ (t|ψ). By the Glivenko-Cantelli theorem, the empirical
distribution formed by {h∗1 , . . . , h∗M } converges uniformly to Ph∗ (t|ψ) almost surely. Since Ph∗ (t|ψ) is
continuous, according to Lemma 21.2 in Vaart (1998), as M → ∞ the quantile estimate h∗(dMγe) converges
to the true quantile h∗γ (ψ) almost surely. By setting γ to be β /2 and (1 − β )/2, we conclude thatββ[h∗(dM 2 e) , h∗(dM(1− 2 )e) ] is a consistent estimator of the credible interval on h∗ .

Lam and Zhou
2.2.2 General Input Distributions
The approach above can be generalized to the case when we have a general input distribution. It follows
essentially the same idea but uses the Dirichlet process instead of the Dirichlet distribution to address thegeneral support. Similar to a Dirichlet distribution, a Dirichlet process can be viewed as a distribution over
distributions, i.e., each sample from a Dirichlet process is a distribution. To facilitate understanding, we
will first give a brief overview of the Dirichlet process (more introductions on Dirichlet process can befound in, for example, Teh 2010, Ghosal 2010). Specifically, let F0 be a distribution over the probability
space Θ and a be a positive number. Then G is Dirichlet process distributed with base distribution F0 andconcentration parameter a, denoted as G ∼ DP(a, F0 ), if(G(A1 ), . . . , G(Ar )) ∼ Dir(aF0 (A1 ), . . . , aF0 (Ar )),for every finite measurable partition {A1 , . . . , Ar } of Θ. Simply put, a Dirichlet process is a stochastic
process that has Dirichlet distributed finite-dimensional marginal distributions, just as a Gaussian process
has Gaussian distributed finite-dimensional marginal distributions.
Just like the Dirichlet distribution is a conjugate prior, the Dirichlet process is also conjugate for
estimating an unknown distribution based on i.i.d. observations. Let X1 , . . . , Xn be i.i.d. observations fromthe unknown input distribution F c . Denote the empirical distribution by F̂(t) = n1 ∑ni=1 I{Xi ≤ t}. If wechoose the prior over F as DP(a, F0 ), then posterior of F given the data ψ = {X1 , . . . , Xn } isF|ψ ∼ DP(a + n,

a
n
F̂).
F0 +
a+n
a+n

Note that the base measure of the posterior is a weighted average of the prior base F0 and the empirical
distribution F̂, and the weight associated with the prior base is proportional to a. Hence, it is easy to see
the role of the parameters F0 and a in the Bayesian sense: the base measure F0 is the mean of the prior, andthe concentration parameter a is the strength of the prior. As a → 0, the posterior becomes the empirical
distribution and it boils down to a non-informative prior. On the other hand, as n, the number of data
points, increases, the posterior is dominated by the empirical distribution, which in turn becomes a close
approximation of the true distribution.
The Dirichlet process is weakly consistent at the true distribution√
F c with a convergence rate 1/ n; or in other words, as the number of observations/data, n, goes to infinity,
the Dirichlet process converges to the true distribution F c in probability (Lo 1983, Ghosal 2010). Thisjustifies the use of Dirichlet process as the posterior distribution to estimate the unknown true distribution.
We are now ready to present our approach to characterizing the posterior distribution of the optimalvalue of (1) given the data. Conceptually that can be done by the following two steps:
1.
2.

n
a
F0 + a+nF̂); then, draw
Sampling: first, draw sample distributions F 1 , . . . , F M from DP(a + n, a+njjji.i.d. samples Yi , . . . ,YN from each F .
Optimization: for each j = 1, . . . , M, solve the optimization problem
(
)
1 N
j
max EF j [h(X; θ )] ≈ ∑ h(Yi ; θ ) ,
θ ∈Θ
N i=1

and denote the optimal value as h∗ j .
The procedure above outputs the sample optimal values {h∗ j }, which form an empirical distribution thatββapproximates the posterior distribution of h∗ given the data. In particular, [h∗(dM 2 e) , h∗(dM(1− 2 )e) ], whereh∗(1) ≤ . . . ≤ h∗(M) , is an estimate of the (1−β ) credible interval for the optimal value. Similar as Theorem 3,
under the assumption that the posterior distribution of h∗ given data is a continuous distribution, we canββshow that [h∗(dM 2 e) , h∗(dM(1− 2 )e) ] converges to the true credible interval almost surely as M → ∞ and N → ∞.

Lam and Zhou
The difficulty with the conceptual procedure above is in the first step, i.e., sampling from the Dirichlet
process. It is impossible to draw a distribution F j directly from the Dirichlet process, since a full
distribution requires storing an infinite amount of information. Therefore, we resort to sampling strategies
that circumvents drawing distributions from the Dirichlet process but rather draws samples that approximates
the distribution F j . Here we describe two such sampling methods. The first, Polya urn scheme (Teh 2010),
treats each F j as a latent distribution and directly draws i.i.d. samples {Yi j } from F j . It is described below(the superscript j is dropped for convenience):••na
F0 + a+nDraw Y1 from a+n
F̂.
ny
a
n
a+n
draw Yi from a+n
F0 + a+n
F̂; with probability a+n+i−1
For i = 2, . . . , N, with probability a+n+i−1
set Yi = y, where ny is the number of previous samples that take value y.

As N goes to infinity, the Polya urn scheme ensures that {Yi } becomes an exact draw from the Dirichlet
process. An alternative and simpler sampling method is to draw {Y1 , . . . ,YN } directly from the distributionana+n F0 + a+n F̂, which is equal to E[F(t)|ψ] and thus a natural Bayesian estimate for F(t). This is essentiallybootstrap sampling from the weighted average of the prior and the empirical distribution. Although thissampling strategy does not give an exact draw from the Dirichlet process, it has the same asymptotic
performance as the exact sampling method (see, e.g. Hjort 1985, Lo 1987).
2.3 Bootstrap Approaches
The Dirichlet method above has a close relation with the nonparametric bootstrap approaches proposed by
Barton and Schruben (1993) and Barton and Schruben (2001) for quantifying the uncertainty in stochastic
simulation. Specifically, they proposed the direct bootstrap approach based on Efron (1982) and the Bayesianbootstrap approach based on Rubin (1981), both resampling empirical distributions from the data set andthen using direct simulation to construct confidence intervals on the simulation output. We will now adapt
their direct bootstrap and Bayesian bootstrap approaches to our setting in quantifying the uncertainty instochastic optimization.
1.

2.

Sampling: for each j = 1, . . . , M,
– If direct bootstrap: draw i.i.d. samples Y1j , . . . ,YNj from F̂. Then F j (t) = N1 ∑Ni=1 I{Yi j ≤ t} isan empirical distribution.
– If Bayesian bootstrap (also called “uniformly randomized resampling”): draw i.i.d. samplesjjj≤ ... ≤; generate ordered samples u(1)≤ . . . ≤ Y(N)Y1j , . . . ,YNj from F̂, order the samples as Y(1)
j
j
j
j
u(N)
from a Uniform(0,1) distribution, and let u(0)= 0. Then F j (t) = ∑Ni=1 (u(i)− u(i−1))I{Y(i)j ≤t} is a smoothed empirical distribution.
Optimization: for each j = 1, . . . , M, solve the optimization problem

max EF j [h(X; θ )],θ ∈Θ

and denote the optimal value as h∗ j .
Direct bootstrap resamples with replacement from the data and assigns equal probabilities to thesamples; whereas Bayesian bootstrap also resamples with replacement but reweights the samples usingjj
uniform spacings, because the joint distribution of the probabilities F(Y(1)), . . . , F(Y(N)) corresponds to thatof the N order statistics from a uniform distribution. The Bayesian bootstrap leads to a more smoothed
empirical distribution than the direct bootstrap.
The two bootstrap resampling schemes are closely related with the Dirichlet methods. When theconcentration parameter a in the Dirichlet process goes to 0 (which is the non-informative case), the Dirichlet
process collapses to a finite-dimensional Dirichlet distribution, and a distribution from the Dirichlet processLam and Zhou
can only have probability mass on the observed data points. Hence, drawing a sample distribution from the
Dirichlet process becomes bootstrap resampling from the data set. More precisely, when a → 0, samplingfrom the Dirichlet process boils down to Bayesian bootstrap, which explains the name of “Bayesian”
bootstrap (Lo 1987); or in other words, samples from the Bayesian bootstrap correspond to discrete
distributions supported at the observed data points with Dirichlet distributed weights. On the other hand,the approximate sampling scheme for Dirichlet process in Section 2.2.2 is equivalent to direct bootstrap as
a → 0. Therefore, both direct and Bayesian bootstrap approaches can be interpreted as Dirichlet methods
with non-informative priors.
3

QUANTIFYING UNCERTAINTY OF THE OPTIMALITY GAP

The three approaches in Section 2 can all be adapted readily to quantify the uncertainty for the optimality
gap between (1) and (2). In this section, we suppose that θ̂ has been picked as our selected solution, andwe are interested in the uncertainty of the difference ∆ = maxθ ∈Θ E[h(X; θ )] − E[h(X; θ̂ )]. We present theadaptation for each of the methods:
Empirical likelihood: Using the notation in Section 2.1, we can obtain a confidence interval for ∆ asfollows:
Theorem 4 Under Assumption 1, with (∇θ h(X; θ0 ), h(X; θ0 )) in Condition 3 replaced by (∇θ h(X; θ0 ), h(X; θ0 )−h(X; θ̂ )), the programs
max / minw maxθ ∈Θ ∑ni=1 wi (h(Xi ; θ ) − h(Xi ; θ̂ ))2subject to−2 ∑ni=1 log(nwi ) ≤ χq,β∑ni=1 wi = 1
wi ≥ 0 for all i(8)where max / min denote a pair of maximization and minimization, contain the true value of ∆ with probability
at least 1 − β asymptotically as n → ∞.
The argument for (4) is similar to the proof of Theorem 1, but now with the profile likelihood definedas()
R(θ , ∆) = max

n

n

n

n

∏ nwi : ∑ wi ∇θ h(Xi ; θ ) = 0, ∑ wi (h(Xi ; θ ) − h(Xi ; θ̂ )) = ∆, ∑ wi = 1, wi ≥ 0 for all i
i=1

i=1

i=1

i=1

since ∆ is now the parameter of interest. R(θ , ∆) satisfies the same behavior as Theorem 1 under themodified version of Assumption 1 in Theorem 4, and the argument follows suit.
In terms of computation, (8) has essentially the same complexity as (3).
Dirichlet method and Bootstrap approaches: In the optimization step, for each j, instead of outputtingh∗ j = maxθ ∈Θ EF j [h(X; θ )], one should outputh∗ j − EF j [h(X; θ̂ )],where EF j [h(X; θ̂ )] is approximated by N1 ∑Ni=1 h(Yi j ; θ̂ ) in the Dirichlet method, or computed according to
the empirical distribution F j in the bootstrap approaches.
4

COMPARISON OF THE PROPOSED METHODS

We have presented three methods for quantifying uncertainty in the optimal value and the optimality gapof a stochastic optimization problem. The EL method is a frequentist approach and based on asymptotics

Lam and Zhouof large samples; it does not require sampling, but needs to solve a two-layer deterministic optimization
problem. The Dirichlet method is fully Bayesian; it requires sampling and needs to solve a significant
number (usually at least 1000 recommended) of single-layer deterministic optimization problems. The
two bootstrap approaches can be viewed as the Dirichlet method with non-informative prior under theexact and approximate sampling schemes respectively, and they also require to solve the same numberof single-layer optimization problems as the Dirichlet method. The computational efficiency of thesemethods is largely determined by the time of solving the two-layer optimization problem or solving alarge number of single-layer problems. For the numerical example that we are presenting next, we findthat solving the two-layer optimization in the EL method appears more efficient than the latter. However,as discussed before, the max-max problem is non-convex and there is no guarantee of reaching a global
optimum. Moreover, there is also a silver lining in spending more computational effort in the Dirichlet andthe bootstrap approaches: these methods output much richer information by providing an estimate of theentire posterior distribution of the optimal value. This approximate posterior distribution can be used for
estimating many quantities of interest, such as statistics of the optimal value (mean, variance, etc.), andBayesian confidence intervals under different confidence levels. In contrast, the EL method is not able toprovide estimates of such statistics, and will have to solve a different two-layer optimization problem eachtime when we change the confidence level in order to get the corresponding confidence interval.
5

NUMERICAL STUDY

We demonstrate the presented methods numerically on a simple example of estimating CVaRα,F c (X), theα-level Conditional-Value-at-Risk (CVaR) of a random variable X, which we assume follows an unknown
distribution F c . This can be rewritten as a stochastic optimization problem:
1+
E[(X − θ ) ] ,(9)min θ +
θ1−αwhere (·)+ is short for max(·, 0).
We assume that F c is a standard normal distribution, and set α = 0.9. Assuming we are given n datafrom the normal distribution, we implement the EL method (EL), Dirichlet method using Polya urn process(Dir), approximate Dirichlet method (AD), direct bootstrap (DB), and Bayesian bootstrap (BB), to obtain
the inferred 95% confidence upper and lower bounds for the optimal value of (9). Moreover, we also
compare with the confidence intervals obtained from central limit theorem and the delta method (Shapiro,Dentcheva, and Ruszczyński 2014, Theorem 5.7), given by"#σ̂ (θ̂ ∗ )∗Ẑ ± z1−β /2 √nwhere z1−β /2 is the critical value of the standard normal distribution with confidence 1−β , θ̂ ∗ is the empirical∗
2 ∗optimal solution, Ẑ ∗ is the empirical optimal value given by (1/n) ∑ni=1 h(Xq i ; θ̂ ), and σ̂ (θ̂ ) is the empiricalstandard deviation of the objective value at the optimal solution given by (1/(n − 1)) ∑ni=1 (h(Xi ; θ̂ ∗ ) − Ẑ ∗ )2 .
Regarding the parameter specifications, we pick the base distribution F0 ∼ N(0, 1) and concentration
parameter a = 0.1 for the Dirichlet process prior (note that F0 is exactly the unknown distribution F c andso we would expect a good performance of the Dirichlet method, which is confirmed below). For both
Dirichlet and bootstrap methods, we pick the number of empirical distributions M = 2000 to be drawn in
the outer layer and the number of samples N = n in the inner layer.
We consider three settings n = 10, 50 and 100. For each setting, we repeat the experiment 100 times,
and note down the empirical coverage probability, mean upper and lower bounds, and the mean and standard
deviation of the interval width for each method. The results are summarized in Table 1. Note that the true
optimal value can be accurately calculated in our setting, and is given by 1.7550.

Lam and Zhoun = 10
EL
Dir
AD
DB
BB
CLT
n = 50
EL
Dir
AD
DB
BB
CLT
n = 100
EL
Dir
AD
DB
BB
CLT

Coverageprobability
0.39
0.53
0.25
0.20
0.20
0.63

Mean lower
bound
0.80
0.73
1.40
1.38
1.37
0.94

Mean upper
bound
1.65
1.91
1.76
1.65
1.65
2.35

Mean interval
width
0.85
1.18
0.35
0.27
0.27
1.41

Standard deviationof interval width
0.56
0.38
0.16
0.20
0.20
1.02

Coverageprobability
0.90
0.88
0.72
0.67
0.67
0.86

Mean lower
bound
1.21
1.21
1.46
1.46
1.46
1.22

Mean upper
bound
2.31
2.20
1.97
1.95
1.95
2.23

Mean interval
width
1.10
0.99
0.50
0.49
0.49
1.01

Standard deviationof interval width
0.40
0.31
0.17
0.19
0.19
0.42

Coverageprobability
0.98
0.95
0.78
0.71
0.72
0.89

Mean lower
bound
1.34
1.31
1.46
1.45
1.45
1.36

Mean upper
bound
2.28
2.13
1.96
1.95
1.95
2.07

Mean interval
width
0.94
0.82
0.50
0.50
0.50
0.71

Standard deviationof interval width
0.27
0.20
0.12
0.14
0.14
0.21

Table 1: Comparison of performance among different methods, EL, Dir: Dirichet method using Polya urn,
AD: approximate Dirichlet, DB: direct bootstrap, BB: Bayesian bootstrap, CLT: classical method
We see that the methods are relatively comparable. For all three cases, EL, Dir and CLT have the highest
coverage probabilities, and among the three EL and Dir show more accurate coverage probabilities than
CLT as n becomes larger. The other methods, namely AD, DB, and BB, produce much narrower intervals,
which is a favorable characteristic that compensates their much lower coverage probabilities. However,their mean lower and upper bounds do not cover the true optimal value 1.7750 when n = 10, suggestingthese methods are not reliable when data size is small. Overall EL and Dir strike a good balance between
coverage probability and interval width when there is a reasonable amount of data. Note that when n issmall Dir is sensitive to choice of the prior (i.e., the base distribution F0 ) and the concentration parametera, but the effect of the prior fades away when n becomes larger or when a goes to zero. In the case when
we do not have much prior information, it is better to use a noninformative prior such as a very flat normal
distribution and a small a; and when we have a good prior, then we can use this prior distribution with arelatively large a.
ACKNOWLEDGMENTS
The first author is grateful to the support by the National Science Foundation under Grants CMMI-1400391
and CMMI-1436247. The second author is grateful to the support by the National Science Foundation underGrant CMMI-1413790 and Grant CAREER CMMI-1453934, and Air Force Office of Scientific Researchunder Grant YIP FA-9550-14-1-0059.
REFERENCES
Barton, R., and L. Schruben. 1993. “Uniform and Bootstrap Resampling of Input Distributions.”. In
Proceedings of the 1993 Winter Simulation Conference, 503–508.

Lam and Zhou
Barton, R., and L. Schruben. 2001. “Resampling Methods for Input Modeling”. In Proceedings of the2001 Winter Simulation Conference, edited by D. J. M. B. A. Peters, J. S. Smith and e. M. W. Rohrer,
372–378.
Cox, D. R., and D. V. Hinkley. 1979. Theoretical Statistics. CRC Press.
Csiszar, I., and G. Tusnady. 1984. “Information Geometry and Alternating Minimization Procedures”.
Statistics & Decisions Supplement Issue (1): 205–237.
Efron, B. 1982. “The Jackknife, the Bootstrap and Other Resampling Plans”. In Volume 38 of CBMS-NSF
Regional Conference Series in Applied Mathematics.
Ferguson, T. 1973. “A Bayesian Analysis of Some Nonparametric Problems”. The Annals of Statistics 1(2): 209–230.
Ghosal, S. 2010. Bayesian Nonparametrics, Chapter 2: The Dirichlet process, related priors and posterior
asymptotics, 35–79. Cambridge Series in Statistical and Probabilistic Mathematics (No. 28). Cambridge
University Press.
Hjort, N. 1985. “Bayesian Nonparametric Bootstrap Confidence Intervals”. Technical Report 240, Department of Statistics, Stanford University.
Kim, S., R. Pasupathy, and S. G. Henderson. 2015. “A Guide to Sample Average Approximation”. In
Handbook of Simulation Optimization, 207–243. Springer.
Lo, A. 1983. “Weak Convergence for Dirichlet Processes”. Sankhya: The Indian Journal of Statistics,
Series A 45 (1): 105–111.
Lo, A. 1987. “A Large Sample Study of the Bayesian Bootstrap”. The Annals of Statistics 15 (1): 360–375.
Mendelson, S. 2003. “A Few Notes on Statistical Learning Theory”. In Advanced lectures on machinelearning, 1–40. Springer.
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.
Rockafellar, R. T., and S. Uryasev. 2000. “Optimization of Conditional Value-at-Risk”. Journal of risk 2:21–
42.
Rubin, D. 1981. “The Bayesian Bootstrap”. Annals of Statistics 9:130–134.
Shapiro, A., D. Dentcheva, and A. Ruszczyński. 2014. Lectures on Stochastic Programming: Modelingand Theory, Volume 16. SIAM.
Shapiro, A., and A. Nemirovski. 2005. “On Complexity of Stochastic Programming Problems”. In Continuous
optimization, 111–146. Springer.
Sion, M. et al. 1958. “On General Minimax Theorems”. Pacific J. Math 8 (1): 171–176.
Teh, Y. 2010. “Dirichlet Process”. Encyclopedia of Machine Learning.
Vaart, A. V. D. 1998. Asymptotic Statistics. Cambridge UK: Cambridge University Press.
Xie, W., B. Nelson, and R. Barton. 2015. “A Bayesian Framework for Quantifying Uncertainty in Stochastic
Simulation”. Operations Research. Accepted.
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.
ENLU ZHOU is an Assistant Professor in the H. Milton Stewart School of Industrial and Systems Engineeringat Georgia Institute of Technology. Her research interests include stochastic control and simulation optimization, with applications towards financial engineering. Her email address is enlu.zhou@isye.gatech.edu
and her web page is http://enluzhou.gatech.edu/.