Proceedings of the 2013 Winter Simulation Conference
R. Pasupathy, S.-H. Kim, A. Tolk, R. Hill, and M. E. Kuhl, eds.

ITERATIVE METHODS FOR ROBUST ESTIMATION UNDER BIVARIATE
DISTRIBUTIONAL UNCERTAINTYHenry Lam

Soumyadip Ghosh

Department of Mathematics and Statistics
Boston University
Boston, MA 02215, USA

Business Analytics and Math Sciences
T.J. Watson IBM Research Center
Yorktown Heights, NY 10598, USA

ABSTRACT
We propose an iterative algorithm to approximate the solution to an optimization problem that arises
in estimating the value of a performance metric in a distributionally robust manner. The optimization
formulation seeks to find a bivariate distribution that provides the worst-case estimate within a specified
statistical distance from a nominal distribution and satisfies certain independence condition. This formulation
is in general non-convex and no closed-form solution is known. We use recent results that characterize
the local “sensitivity” of the estimation to the distribution used, and propose an iterative procedure on thespace of probability distributions. We establish that the iterations of solutions are always feasible and thatthe sequence is provably improving the estimate. We describe conditions under which this sequence canbe shown to converge to a locally optimal solution. Numerical experiments illustrate the effectiveness of
this approach for a variety of nominal distributions.
1

INTRODUCTION

In many computational settings in engineering and management operations, the model assumed on theunderlying system components is based on some convenient approximation of the truth. When the actual
probability distribution that governs the model is different from the assumption, the decision output canpotentially deviate to a substantial degree. This is commonly referred to as model uncertainty or ambiguity,and there is a sizeable literature devoted to studying this issue. For example, model selection methods basedon classical and Bayesian statistics (Draper 1995; Chick 2001) and construction of input-adjusted confidence
intervals have been proposed in large simulation studies. The latter can be done via delta-method approach(Zouaoui and Wilson 2004) or more sophisticated techniques, such as metamodel-assisted bootstrapping(Barton, Nelson, and Xie 2010; Barton 2012). In decision analysis, various forms of robust optimization
have been introduced. The main streams of work include deterministic robust optimization that represents
model ambiguity as deterministic uncertainty sets (Bertsimas, Brown, and Caramanis 2011; Ben-Tal and
Nemirovski 2002), and distributionally robust optimization, where the underlying probability distributions
are assumed to be ambiguous but come from a convenient family (Lim, Shanthikumar, and Shen 2006; Hu
and Hong 2013). The latter approach is also particularly popular among the control community (Hansen
and Sargent 2011), and recently, it has been applied to derive so-called robust Monte Carlo and has beenused to compute robust bounds for performance measures that arise in the context of finance (Glasserman
and Xu 2012).
In this paper, we consider a general form of robust estimation problems in the spirit of Glasserman
and Xu (2012), which aims to generate bounds for performance measures that are subject to model
uncertainty. Typically, these problems can be formulated as constrained maximization and minimization of
the performance measure, say E[h(X)], where h(·) is the cost function that is known (but not necessarily inclosed-form) and X is the random component. The optimization is over the probability distribution P of X

Ghosh and Lamthat generates the expectation E[·]. Constraints are placed on the probability P to lie within an uncertainty
set; one common example of such uncertainty set is a ball around a benchmark distribution, say, P0 , measured
via a certain statistical distance, such as Chi-Square or Kullback-Leibler divergence. For instance, the costfunction h can be the indicator of whether a portfolio with underlying stock price X hits a certain levelof loss. While the process driving the stock is not known completely, portfolio managers typically havecertain degree of confidence about the right distribution. The idea of the maximization program under the
distribution constraint is to find the worst possible tail probability under such confidence. Similarly, theminimization program can be viewed as the best possible scenario. Depending on the application, these
programs can have different interpretations and variations.
In many situations, a careful set-up of the constraints and the program formulation leads to tractable
solution. For example, consider maximizing E f [h(X)] subject to a Kullback-Leibler constraint D(Pf kP0 ) ≤ η,
where D denotes the distance between the distribution Pf on X and the benchmark P0 , and η captures
the level of confidence around the benchmark. This formulation has a neat explicit characterization of the
optimal solution, which can be represented as an exponential change of measure on the cost h(X) and an
exponential parameter that is related to the confidence level η (Hansen and Sargent 2011). Analogous
results have been shown to hold under other statistical distances, such as χ 2 -distance (Chen and Lam 2013),
with different forms of representation for the optimal changes of measure.
Our main focus of this paper is a natural scenario, especially in the context of simulation, where suchtractable formulation breaks down. In particular, we consider restrictions that are placed on the marginal
distributions of the random variate X, defined as lower-dimensional (typically univariate) distributions that
describe the probability law for sub-sets of components of X. Previous work, such as Delage and Ye(2010) and Bertsimas, Popescu, and Sethuraman (2000), has considered marginal moment constraints inrobust estimation formulation. In this article, we shall study marginal constraints based on χ 2 -distanceon a bivariate random vector that generates the objective function E f [h(X1 , X2 )], together with a natural
constraint that the individual components Xi are i.i.d.. In other words, we restrict the space of feasible
distributions P to those that are product form bivariate:max
E f [h(X,Y )]subject to X,Y ∼ Pf , i.i.d.
χ 2 (Pf , P0 ) ≤ η
Pf ∈ P0(1)(the min formulation is analogous) where the constraint is on the χ 2 -distance between Pf and P0 , given byχ 2 (Pf , P0 ) = E0 (dPf /dP0 − 1)2 , with dPf /dP0 being the Radon-Nikodym derivative of Pf with respect toP0 . The set P0 represents all Pf that are absolutely continuous with respect to P0 , a condition needed forχ 2 -distance to be defined. The important feature here is that X and Y are assumed to be i.i.d. and they
are both generated according to Pf , which is uncertain. The difficulty with the optimization problem (1)is that it is usually non-convex, and that no explicit form of the optimal solution is available.
Such formulation, in its more general form, appears commonly in stochastic systems, such as queueing
networks, when many customers in the system each contribute as individual random sources and theassociated vector of inputs (inter-arrival times and service-times) across the customers are i.i.d. The max
and the min optimal values for the formulation type in (1) will together give a bounding interval for the
performance measure subject to the specified constraints. As another motivating example, consider thefollowing inventory model that will be studied in the numerical experiments section:
Example 1 A retailer has to decide on the size of order to place with a manufacturer that operates abatch-processing facility. The retailer requires that the order quantity q be sufficient to meet demand in two
successive time periods, where each period’s demand di , i = 1, 2 is i.i.d. The distribution of di is determined
to be P0 by running the Chi-Square test statistic over available data, which returns an uncertainty estimateof χ 2 (P, P0 ) ≤ η. A profit of p units is realized on every sale, but the manufacturer charges a cost of cunits for unsold inventory carried over at the end of each period. The net expected profit the retailer stands

Ghosh and Lam
to make under the nominal distribution P0 is
E0 h(d1 , d2 , q) = pE0 [ min{d1 + d2 , q}] − cE0 [(q − d1 )+ ] − cE0 [(q − d1 − d2 )+ ].

(2)

Note that this decision problem is distinct from the classical newsvendor problem because of the presenceof the extra carry-over cost charged at the end of the first period, represented in (2) by the middle term.
Such formulations are typical of batch-produced goods such as steel slabs and hot-rolled coils from steel
manufacturing plants.
We view (1) as an initial attempt to tackle more general robust estimation problems. With more
involved analysis, our framework can be extended to the case where the random vector X has more than two
i.i.d. components. It is also plausible that other statistical distances, such as Kullback Leibler divergence,
can be used, as long as they possess a similar type of “local expansion” behavior that we shall discuss
momentarily. Another generalization one might consider is the inclusion of moment constraints consideredin, e.g., Delage and Ye (2010); we believe our method can be extended for this problem instance too.
The main goal of this paper is to derive a descent-type iterative procedure to find local optimum of(1). A tool that we shall utilize is the optimality characterization of (1) when η is chosen small enough,as a local expansion in terms of η, which was recently studied in Lam (2013) and Chen and Lam (2013).
The result states that, for a bivariate cost function h, as η → 0,√
(3)
max E f [h(X,Y )] = E0 [h(X,Y )] + sd0 (H0 (X)) η + o(η)
χ 2 (Pf ,P0 )≤η

where H0 (x) is a “symmetrization” form of the cost function h, defined as
H0 (x) = E0 [h(X,Y )|X = x] + E0 [h(X,Y )|Y = x]and sd0 (·) denotes the standard deviation under P0 . Note that sd0 (·) only acts on a single X in (3). This“symmetrization” arises in a “product rule” when differentiating the Lagrangian formulation in (1) (Lam2013); alternatively, it can also be derived using a variation of Hoeffding decomposition (Serfling 2009). The
key observation, in either method of derivation, is that on a local level, the optimal probability distribution
within χ 2 (Pf , P0 ) ≤ η can be characterized by the change of measureL∗ (x) :=dPf
H0 (x) − E0 H0 (X) √(x) = 1 +η
dP0
sd0 (H0 (X))(4)

assuming that sd0 (H0 (X)) > 0 and η is small enough, under certain regularity conditions on h (such as
boundedness, which we shall state in the sequel). In other words, the approximation (3) can be achievedby√
E0 [h(X,Y )L∗ (X)L∗ (Y )] = E0 [h(X,Y )] + sd0 (H0 (X)) η + o(η)(5)
We propose an iterative scheme whose local move is based on (5); Section 2 provides an outline of the
proposed iterative procedure. In Section 3, we will prove some theoretical guarantees on the procedure.
We show that our scheme is always feasible and ascending for the maximization formulation (1). Notethat our proof does not need to know (3); it merely provides a guidance to find an ascending local move.
Under mild boundedness conditions on the function h, we establish in Sub-section 3.2 that this sequencemust then converge. Sub-section 3.3 provides a local optimality guarantee on the limit under stronger
assumptions on the function h.
The iterative procedure described in Section 2 can be computed exactly for any benchmark distribution
P0 that has a discrete support. Section 4 describes a heuristic scheme to apply this procedure using a particle
approach to the case where the support of the nominal distribution is continuous. Section 5 applies thisheuristic to the inventory model in Example 1 under two assumptions of the benchmark P0 and illustrates
the efficiency of this approach.

Ghosh and Lam
2

THE ALGORITHM

We state an algorithm to approximate the solution of (1). Note that we can rewrite (1) in terms of thelikelihood ratio L := dPf /dP0 as
max
E0 [h(X,Y )L(X)L(Y )]subject to E0 (L − 1)2 ≤ η
L∈L

(6)

where the χ 2 -distance is expressed as E0 (L − 1)2 , and L = {L : E0 [L] = 1, L ≥ 0 a.s.}. It is easy to seethat in general the program (6) is non-convex. As such, it is impossible to find a scheme that is guaranteed
to approach the global optimum. We shall instead focus on finding the local optimum of (6). One heuristic
to find the global optimum is to, for instance, run our algorithm at different initial points in the feasible
set, in the space L . Also, throughout our analysis we assume that h is bounded.
Suppose we start from the initial distribution P0 . At each step k, suppose that the current distribution is
updated to Pk , and the distance from the benchmark is Dk := χ 2 (Pk , P0 ). Moreover, let Lk := dPk /dP0 and
Hk (x) = Ek [h(X,Y )|X = x] + Ek [h(X,Y )|Y = x]be the “symmetrization” of the system under Pk , with Ek [·] being the expectation taken with respect to the
distribution Pk . We also let H̄k (x) = Hk (x) − Ek Hk (X) as the centered version of Hk (x). Then, the following
quantities are needed to define the next movement for Pk+1 :1. DefineAk := −
2. Definesdk (Hk (X))32(Ek [h(X,Y )H̄k (X)H̄k (Y )])
sdk (Hk (X))
Bk := min −: x ∈ supp(Pk ) with H̄k (x) < 0
H̄k (x)3. Definep(E0 [Lk (X)2 H̄k (X)])2 + E0 [Lk (X)2 H̄(X)2 ](η − Dk ) − E0 [Lk (X)2 H̄k (X)]
Ck := sdk (Hk (X))
E0 [Lk (X)2 H̄k (X)2 ]where supp(Pk ) denotes the support of Pk . We introduce a step size parameter δk+1 , which defines theχ 2 -distance that one should move in the next (k + 1)-th step. To compute this parameter, we divide into
two cases. Suppose that Hk (X) is non-degenerate under Pk ,1. If Ek [h(X,Y )H̄k (X)H̄k (Y )] < 0, thenp
δk+1 = min{Ak , Bk ,Ck }2. If Ek [h(X,Y )H̄k (X)H̄k (Y )] ≥ 0, thenp
δk+1 = min{Bk ,Ck }Next, we let
νk (x) =H̄k (x).
sdk (Hk (X))

With δk+1 and νk (·), we define our change of measure at the (k + 1)-th step asp
Lk+1 (x) = Lk (x) 1 + νk (x) δk+1

(7)

Ghosh and Lam
i.e. the probability distribution at the (k + 1)-th step ispPk+1 (x) = Pk (x) 1 + νk (x) δk+1

(8)This updates the probability distribution. After each iteration, the values of Dk+1 and Hk+1 (x) are also
updated. Repeat the procedure until either
1. Dk = η or
2. Hk (x) is degenerate under Pk .
We shall provide a few remarks on the algorithm. First, we explain some intuition on the purpose of
the quantities Ak , Bk and Ck :1. Ak ensures that the algorithm is strictly ascending from step k to k + 1. As we shall see, the form of
Ak comes from the analysis of a quadratic form that represents the increment from step k to k + 1.
2. Bk ensures that the next probability distribution Pk+1 is a valid distribution, i.e. the likelihood ratiothat represents the change of measure at each step has to be non-negative.
3. Ck ensures that the next probability distribution Pk+1 is feasible, i.e. χ 2 (Pk+1 , P0 ) ≤ η.
We note that the form of the change of measure (7) resembles the best local movement as in (5). The
parameter δk captures the confidence level such that the linear approximation in (3) holds. If δk is chosen
small enough, we can then argue that the movement is ascending in the objective value. The quantities Ak ,Bk and Ck combine to provide such guarantee.
In the next section, we will show rigorously that the algorithm is always strictly ascending, until one of
the two stopping criteria holds. The first criterion means that a boundary point of the set χ 2 (Pf , P0 ) ≤ η isreached. The second criterion, roughly speaking, is a condition for reaching a “stationary point”. We will
discuss in Sub-section 3.3 a semi-definite condition that guarantees that this “stationary point” is a local
maximum in certain sense (the definition of local maximum is not entirely obvious, since χ 2 -distance is
not a metric).
3

THEORETICAL PROPERTIESIn this section we discuss a few properties of our algorithm. These include the strict ascendency of thealgorithm, an optimality condition on the end value of the algorithm, and a convergence property of the“derivative” term sdk (Hk (X)) in (5) in a sense that we shall discuss.
3.1 Ascendency Guarantee
The first property can be summarized as follows:
Theorem 1 Suppose the cost function h is bounded. The algorithm in Section 2 is always feasible, i.e.
χ 2 (Pk , P0 ) ≤ η, and strictly ascending, i.e. Ek+1 [h(X,Y )] > Ek [h(X,Y )], until one of the two stopping
criteria is met.
p
Proof.
Note that δk+1 is the minimum of either the three terms Ak , Bk and Ck or the two terms Bk andCk , under the corresponding conditions in Section 2. Related to the discussion before, we shall show that
1. Ak ensures that the algorithm is strictly ascending from step k to k + 1.
2. Bk ensures that the next probability distribution Pk+1 is a valid distribution.
3. Ck ensures that the next probability distribution Pk+1 is feasible.
We shall show these one by one:Ghosh and LampAscendency. First, by choosing δk+1 ≤ Ak , we have




H̄k (Y ) p
H̄k (X) p
δk+1
1+
δk+1
Ek+1 [h(X,Y )] = Ek h(X,Y ) 1 +sdk (Hk (X))
sdk (Hk (Y ))Ek [h(X,Y )(H̄k (X) + H̄k (Y ))] pδk+1= Ek [h(X,Y )] +sdk (Hk (X))
Ek [h(X,Y )H̄k (X)H̄k (Y )]+δk+1
Vark (Hk (X))
Ek [Ek [h(X,Y )|X]H̄k (X)] + Ek [Ek [h(X,Y )|Y ]H̄k (Y )] p
δk+1
= Ek [h(X,Y )] +sdk (Hk (X))
Ek [h(X,Y )H̄k (X)H̄k (Y )]+δk+1
Vark (Hk (X))
p
Ek [h(X,Y )H̄k (X)H̄k (Y )]= Ek [h(X,Y )] + sdk (Hk (X)) δk+1 +δk+1Vark (Hk (X))(9)

by using the i.i.d. assumption on X and Y in the last equality. Suppose Ek [h(X,Y )H̄k (X)H̄k (Yp)] ≥ 0, thenpicking a positive δt+1 will make the algorithm ascending; if Ek [h(X,Y )H̄k (X)H̄k (Y )] < 0, then δk+1 ≤ Ak
guarantees thatp
Ek [h(X,Y )H̄k (X)H̄k (Y )]sdk (Hk (X)) δk+1 +δk+1 ≥ 0.
Vark (Hk (X))
In fact, pickingpsdk (Hk (X))3
δk+1 = −2Ek [h(X,Y )H̄k (X)H̄k (Y )]maximizes the increase in Ek [h(X,Y )]. Hence the algorithm is strictly ascending.
Validity of probability distribution. Note that we have defined in (8) thatpdPk+1
= 1 + νk (x) δk+1 .
dPk
To guaranteepthat Pk+1 is valid, we pmust have Ek [dPk+1 /dPk ] = 1, which is easily seen to hold since
Ek [1 + νk (X) δk+1 ] = 1 + Ek [νk (X)] δk+1 = 1. Secondly, we also need δk+1 to satisfy the constraintp
1 + νk (x) δk+1 ≥ 0
for all x ∈ supp (Pk ). We need only scrutinize the x that have νk (x) < 0, or H̄k (x) < 0, and make sure thatp1
δk+1 ≤ −νk (x)
for all such x. This gives Bk .
Feasibility. Finally, we need to ensure that Dk+1 = χ 2 (Pk+1 , P0 ) ≤ η. Note thatp
Dk+1 = χ 2 (Pk+1 , P0 ) = E0 [Lk (X)2 (1 + νk (x) δk+1 )2 ] − 1p= E0 [Lk (X)2 + 2Lk (X)2 νk (X) δk+1 + Lk (X)2 νk (X)2 δk+1 ] − 1p= E0 [Lk (X)2 ] − 1 + 2E0 [Lk (X)2 νk (X)] δk+1 + E0 [Lk (X)2 νk (X)2 ]δk+1
p= Dk + 2E0 [Lk (X)2 νk (X)] δk+1 + E0 [Lk (X)2 νk (X)2 ]δk+1 .

Ghosh and Lam
Hence we needp
Dk + 2E0 [Lk (X)2 νk (X)] δk+1 + E0 [Lk (X)2 νk (X)2 ]δk+1 ≤ η.
p
Since this is quadratic and convex in δk+1 , we need thatp
δk+1 ≤ ξ

where ξ is the positive root of the equationp
Dk + 2E0 [Lk (X)2 νk (X)] δk+1 + E0 [Lk (X)2 νk (X)2 ]δk+1 = η.
This givespδk+1 ≤

p(E0 [Lk (X)2 νk (X)])2 + E0 [Lk (X)2 νk (X)2 ](η − Dk ) − E0 [Lk (X)2 νk (X)].
E0 [Lk (X)2 νk (X)2 ]Note that the right hand side above is Ck . This will ensure that Dk+1 ≤ η.

3.2 Convergence
We shall also discuss a convergence property of the algorithm. Let us focus on how the algorithm behavesin the interior points of the set χ 2 (Pf , P0 ) ≤ η. In particular, we have the following:
Proposition 1 Consider a bounded cost function h. Suppose η = ∞, i.e. the constraint χ 2 (Pf , P0 ) ≤ η
is removed from (1). Then our algorithm will converge to a probability distribution P∗ that satisfies
Var∗ (H ∗ (X)) = 0 (where Var∗ (·) and H ∗ (·) are defined as the variance and symmetrization of h under P∗ ).
Proof.p Note that in this case, the quantity
p Ck in our algorithm is not needed. If Ek [h(X,Y )H̄k (X)H̄k (Y )] < 0,then δk+1 = min{Ak , Bk }, otherwise δk+1 = Bk .
We shall prove by contradiction. Suppose that the algorithm does not lead the quantity Vark (Hk (X))
to 0, then there must exist a subsequence {ki }i=1,2,... such that Varki (Hki (X)) > ε for all i for some ε > 0.
Now, from the calculation in (9), we know the increment Ek+1 [h(X,Y )] − Ek [h(X,Y )] isp
Ek [h(X,Y )H̄k (X)H̄k (Y )]δk+1 .
sdk (Hk (X)) δk+1 +Vark (Hk (X))(10)

Consider the step ki , and we shall analyze (10) explicitly. We shall divide into three cases:p
1. Eki [h(X,Y )H̄ki (X)H̄ki (Y )] < 0 and δki +1 = Aki : We have (10) equal to
Ek [h(X,Y )H̄ki (X)H̄ki (Y )] 2
sdki (Hki (X))Aki + i
Aki
Varki (Hki (X))
sdki (Hki (X))3= sdki (Hki (X)) −2Eki [h(X,Y )H̄ki (X)H̄ki (Y )]2Eki [h(X,Y )H̄ki (X)H̄ki (Y )]sdki (Hki (X))3+−
Varki (Hk (X))
2Eki [h(X,Y )H̄ki (X)H̄ki (Y )]2
Varki (Hki (X))= −4Eki [h(X,Y )H̄ki (X)H̄ki (Y )](11)

But since h is bounded, so is −Eki [h(X,Y )H̄ki (X)H̄ki (Y )], and since we assume Varki (Hki (X)) > ε,we have (11) bounded away from 0 (and positive).

Ghosh and Lam
2. Eki [h(X,Y )H̄ki (X)H̄ki (Y )] < 0 andthis case, since Aki ≥ Bki , we have−p
δki +1 = Bki : Let H̄k∗ := ess infx∈supp(Pk ):H̄k (x)<0 H̄k (x) < 0. In

sdki (Hki (X))3
sdk (Hk (X))≥ − i ∗i
2Eki [h(X,Y )H̄ki (X)H̄ki (Y )]H̄kiwhich givesH̄k∗i Varki (Hki (X))− Eki [h(X,Y )H̄ki (X)H̄ki (Y )] ≤ −2Note that (10) is equal to(12)

Eki [h(X,Y )H̄ki (X)H̄ki (Y )] 2
Bki
Varki (Hki (X))!
Ek [h(X,Y )H̄ki (X)H̄ki (Y )] Varki (Hki (X))
sdki (Hki (X))= sdki (Hki (X)) −+ i∗
Varki (Hki (X))
H̄ki(H̄k∗i )2
sdki (Hki (X))Bki += −Varki (Hki (X)) Eki [h(X,Y )H̄ki (X)H̄ki (Y )]+H̄k∗i(H̄k∗i )2≥ −Varki (Hki (X))
2H̄k∗i(13)where the last inequality follows from (12). Since Varki (Hki ) > ε and H̄k∗i is bounded from above
since h is bounded, we have (13) bounded away fromp0.
3. Eki [h(X,Y )H̄ki (X)H̄ki (Y )] ≥ 0: In this case we have δki +1 = Bki . Similar to above, (10) is equaltoEki [h(X,Y )H̄ki (X)H̄ki (Y )] 2
B ki
Varki (Hki (X))
Varki (Hki (X)) Eki [h(X,Y )H̄ki (X)H̄ki (Y )]
+
= −
H̄k∗i
(H̄k∗i )2
sdki (Hki (X))Bki +(14)

Since in this case Eki [h(X,Y )H̄ki (X)H̄ki (Y )] ≥ 0, it follows similarly that (14) is bounded away from0.
In conclusion, for every step in the sequence ki , the increment Eki +1 [h(X,Y )] − Eki [h(X,Y )] is always
bounded away from 0. Hence Eki [h(X,Y )] grows in an unbounded manner. But we have assumed that his bounded, hence a contradiction.
3.3 Nesting and Optimality
In this sub-section we shall discuss a simple optimality condition of our algorithm. Before so, we emphasizethat the use of χ 2 -distance as a measurement of disparateness between distributions, our decision variables,
poses some non-standard features to the definition of optimality. At the foremost, χ 2 -distance is not a
proper metric: it does not satisfy triangle inequality nor commutativity. This itself does not pose direct
issue for using iterative schemes. However, it leads to a “one-way” update of Pk , or a nesting behavior interms of absolute continuity. Namely, the sequence of update {Pk }k=0,1,2,... satisfies Pk+1 ≺ Pk , where ≺denotes “absolutely continuous with respect to”. In other words, the support of Pk shrinks at each step.
We shall call P∗ a nested locally optimal distribution if for any P ≺ P∗ , we have E[h(X,Y )] ≤ E ∗ [h(X,Y )].
We call P∗ a nested strict locally optimal distribution if the inequality is strict. The next proposition gives
a condition for local optimality, if the algorithm stops at an interior point of the feasible region. It requires
a semi-definite property on the cost function h.

Ghosh and Lam
Proposition 2 Consider a probability distribution P∗ that satisfies Var∗ (H ∗ (X)) = 0. If h(x, y) is negative
semi-definite in supp(P∗ ), in the sense thatZh(x, y)dε(x)dε(y) ≤ 0(15)

for any signed measure ε on supp(P∗ ) with bounded variation, then P∗ is a nested local maximum of the
program (1). Supposing that h(x, y) is negative definite, i.e. the inequality (15) is strict for any non-zeromeasure, then P∗ is a nested strict local maximum.
Proof.
Suppose P∗ satisfies Var∗ (H ∗ (X)) = 0. Consider any P ≺ P∗ in a neighborhood of P∗ . The proof
follows from the following decomposition (with some abuse of notation)Z

E[h(X,Y )] =h(x, y)dP(x)dP(y)Zh(x, y)d(P∗ (x) + (P − P∗ )(x))d(P∗ (y) + (P − P∗ )(y))

Zh(x, y)dP∗ (x)dP∗ (y) +

=
=

Z

+

Z

h(x, y)dP∗ (x)d(P − P∗ )(y) +Zh(x, y)dP∗ (y)d(P − P∗ )(x)h(x, y)d(P − P∗ )(x)d(P − P∗ )(y)= E ∗ [h(X,Y )] +Z

H ∗ (x)d(P − P∗ )(x) +Zh(x, y)d(P − P∗ )(x)d(P − P∗ )(y)where the last equality follows since x and y are dummies for the variable X. Hence if Var∗ (H ∗ (X)) = 0,then H ∗ (X) is degenerate and the second term above vanishes. If h is negative semi-definite in supp(P∗ ),then we have E[h(X,Y )] ≤ E ∗ [h(X,Y )], and hence P∗ is a nested local maximum. Similarly, if h is negative
definite, then P∗ is a nested strict local maximum.
Note that the proof above resembles the optimality conditions in Euclidean space, the main difference
now being that we are working on the space of distributions. The symmetrization H ∗ (x) can be interpreted
as the “derivative” in such space (see also Lam (2013)), and the semi-definite condition is also analogous.
4

SIMULATION HEURISTIC FOR CONTINUOUS P0

This section describes an algorithm for efficiently calculating the distributions implied by the iterative
procedure described above. The calculations are exact for benchmark distributions that have finite discrete
support. For the case when the support of the benchmark distribution P0 is continuous, we suggest a naı̈ve
approximation procedure of approximating the original benchmark P0 with a discrete distribution P̂0 thatequi-weighs a set of N0 samples generated from P0 in an i.i.d. fashion.
H EURISTIC 1: ROBUST E STIMATION A PPROXIMATION A LGORITHM

1.

2.

Given: benchmark distribution P0 , distribution χ 2 -discrepancy target η, oracle to measure performance metric h(x1 , x2 )If benchmark P0 has support X that is discrete, set Xˆ = X . Else, set Xˆ = {xi , i = 1, . . . , N0 },where each xi is sampled from the benchmark distribution P0 . Let P̂0 = {1/N0 , . . . , 1/N0 } be an
N0 -dimensional row-vector that represents the approximate benchmark distribution over support-set
Xˆ . L0 = {1, . . . , 1}, a N0 −dimensional row-vector that records the likelihood ratio P̂k /P̂0 . Record
h0 = {h(x1 , x j ), ∀xi , xi ∈ Xˆ }, the matrix of values taken by the performance oracle h over Xˆ × Xˆ .
For k = 0, 1, 2, . . .

Ghosh and Lam(a) Calculate Ek h = P̂k h0 P̂kt .
(b) Calculate Dk = χ 2 (P̂k , P̂0 ) = ∑Ni 0 (Lk (i) − 1)2 ∗ P̂k (i). If Dk = η, proceed to Step 3.
(c) Calculate Hk = (h0 + ht0 )P̂k , the symmetric measure under the kth density iterate.
(d) Calculate Ek Hk and sdk (Hk ), the mean and standard deviation of the symmetric measure underP̂k . If sdk (Hk ) = 0, then proceed to Step 3.
(e) Calculate the three limits Ak , Bk and Ck on the distributional distance δk as described in Section 2.
Set δk = min{Ak , Bk ,Cproceed to Step 3.
k }. If δk = 0, p
Hk (i)−Ek Hkδk , ∀i = 1, . . . , N0 .
(f) Set Lk+1 (i) = Lk (i) ∗ 1 + sdk (Hk )3.

5(g) Set P̂k+1 (i) = P̂0 (i) ∗ Lk+1 (i), ∀i = 1, . . . , N0 .
Return P̂k and Ek h as estimates of the desired worst-case density and the performance metric.

NUMERICAL RESULTS

This section analyzes Example 1 that models a retailer of batch-produced goods. Recall that the expected
profit the retailer stands to make under the nominal distribution P0 is given by (2) to be:
E0 h(d1 , d2 , q) = pE0 [ min{d1 + d2 , q}] − cE0 [(q − d1 )+ ] − cE0 [(q − d1 − d2 )+ ].

(16)

A retailer has to decide on the size of the order q to place with a manufacturer that operates a batch-processing
facility. The classical approach to this problem is to first determine the benchmark distribution P0 fromavailable data, and then set the order quantity to be the optimizer of the stochastic optimization problem
that maximizes the realized profit:q∗ = argminq E0 h(d1 , d2 , q).

(17)

Differentiating (16) with respect to q gives us the following optimization criterion:∞∂ E0 h(d1 , d2 , q)
= p − cP0 (q) − (p + c)
P0 (q − x)dP0 (x) = 0.
∂q
−∞

Z

where we abuse notation to denote P0 also as the (identical) distribution function of d1 and d2 . Suppose
the nominal distribution P0 is exp(λ ) and p = rc, r > 1. Then, the optimal q∗ is the unique solution to theequation e−λ q (r + 2 + λ (r + 1)q) = 2. Taking the second derivative of (16) with respect to q and settingq = q∗ confirms that this is a maximizer. For instance, if λ = 1 and r = 3, then q∗ = 1.812. As remarked
earlier, the problem (17) is distinct from the newsvendor model because of the presence of the additional
first-period cost term in (16). In its absence, the optimal order quantity from the newsvendor model would
be the p/(p + c)-th quantile of the distribution of the sum (d1 + d2 ). In the case where p/c = r = 3 andP0 ∼ exp(1), (d1 + d2 ) is an Erlang distribution, and the optimal newsvendor order quantity is 2.6926.
Thus, the optimal order quantity falls to 1.812 because of the additional carryover cost in the first period.
We apply our algorithm in Section 2. Since the exponential distribution P0 is continuous, we use
Heuristic 1 in Section 4 to encode the distribution. For illustrative purpose, we keep track of the discrepancy
distance χ 2 (Pk , P0 ) instead of putting a hard constraint that it is bounded by η (when the constraint is
imposed, the algorithm will stop whenever the η-boundary is reached). By Proposition 1, we know thateventually the algorithm should stop when sdk (Hk (X)) = 0. Figure 1(a) plots the densities of the Pk iterates
generated by applying Heuristic 1 until sdk (Hk (X)) is within a small tolerance around 0. The initial sample
size was set to N0 = 10, 000. The performance measure (16) rises from E0 [h(X,Y )] = 5.941 to the optimal
E ∗ [h(X,Y )] = 9.076, where χ 2 (P∗ , P0 ) = 0.33. (Figure 1(a) plots a smoothed version of the empiricaldensity maintained by Heuristic 1, where the smoothing is obtained by binning the N0 support points
into 50 equi-sized bins and then constructing a Gaussian kernel density with an appropriate smoothingGhosh and Lam
0.5

0.25

0.45
0.4

0.2

0.35
0.3

0.15

0.25
0.2

0.1

0.15
0.1

0.05

0.05
0

0

1

2

3

4

5

6

7

8

9

0

10

0

0.1

0.2

0.3

(a) P0 ∼ exp(1)

0.4

0.5

0.6

0.7

0.8

0.9

1

(b) P0 ∼ U(0, 1)

Figure 1: Densities Generated by Iterative Heuristic 1, Smoothed Using Gaussian Kernels.
parameter. This smoothing is meant as a visual aid and has no bearing on the calculations in Heuristic 1.)
Figure 1(b) plots the densities of the iterates for the same problem but with P0 ∼ U[0, 1], q∗ = 1.524 andr = 8.
Of perhaps more interest is Figure 2 that plots the objective function Ek [h(X,Y )] in (2), the discrepancy
measure χ 2 (Pk , P0 ) and the “gradient” sdk (Hk (X)) associated with the k-th iteration of Heuristic 1 for boththe exponential and the uniform benchmark cases. The results show that Heuristic 1 increasingly tends toconcentrate the mass in a subset of the original sample-set of N0 points that reduces the value of sdk (Hk (X)).
In the exponential-P0 case, the method was able to identify a distribution with sdk (Hk (X)) = 0 within k = 22
iterations. The identified solution is however only an approximation of the true optimal density (when nodiscrepancy constraint is imposed), since the starting P0 used in Heuristic 1 is a sample set generated fromthe original benchmark distribution.
12
Exponential
Uniform

1

10

Exponential
Uniform

11
10
9
0

10

8
7
0

10

6
−1

10

5
4
3

0

5

10

15

(a) χ 2 (Pk , P0 )

20

25

30

2

Exponential
Uniform
0

5

10

15

20

25

(b) Performance Metric Ek h(X,Y )

−2

30

10

0

5

10

15

20

25

30

(c) sdk (Hk (X))

Figure 2: Performance of Heuristic 1 as Iteration Count Grows, in the Exponential and Uniform Benchmark
Cases.
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. M.
Uhrmacher, 1–12. Piscataway, New Jersey: Institute of Electrical and Electronics Engineers, Inc.

Ghosh and Lam
Barton, R. R., B. L. Nelson, and W. Xie. 2010. “A framework for input uncertainty analysis”. In Proceedingsof the 2010 Winter Simulation Conference, edited by B. Johansson, S. Jain, J. Montoya-Torres, J. Hugan,and E. Yucesan, 1189–1198. Piscataway, New Jersey: Institute of Electrical and Electronics Engineers,
Inc.
Ben-Tal, A., and A. Nemirovski. 2002. “Robust optimization–methodology and applications”. Mathematical
Programming 92 (3): 453–480.
Bertsimas, D., D. B. Brown, and C. Caramanis. 2011. “Theory and applications of robust optimization”.
SIAM Review 53 (3): 464–501.
Bertsimas, D., I. Popescu, and J. Sethuraman. 2000. “Moment problems and semidefinite optimization”.
In Handbook of semidefinite programming, 469–509. Springer.
Chen, X., and H. Lam. 2013. “Robust Assessment of Model Uncertainty for Steady-State Estimators”.
Working Paper.
Chick, S. E. 2001. “Input distribution selection for simulation experiments: accounting for input uncertainty”.
Operations Research 49 (5): 744–758.
Delage, E., and Y. Ye. 2010. “Distributionally robust optimization under moment uncertainty with application
to data-driven problems”. Operations Research 58 (3): 595–612.
Draper, D. 1995. “Assessment and propagation of model uncertainty”. Journal of the Royal StatisticalSociety. Series B (Methodological):45–97.
Glasserman, P., and X. Xu. 2012. “Robust Risk Measurement and Model Risk”. Available at SSRN 2167765.
Hansen, L. P., and T. J. Sargent. 2011. Robustness. Princeton university press.
Hu, Z., and J. Hong. 2013. “Kullback-Leibler Divergence Constrained Distributionally Robust Optimization”.
Working Paper.
Lam, H. 2013. “Robust Sensitivity Analysis for Stochastic Systems”. arXiv preprint arXiv:1303.0326.
Lim, A. E., J. G. Shanthikumar, and Z. M. Shen. 2006. “Model uncertainty, robust optimization, andlearning”. Tutorials in Operations Research: Models, Methods, and Applications for Innovative Decision
Making:66–94.
Serfling, R. J. 2009. Approximation theorems of mathematical statistics, Volume 162. Wiley-Interscience.
Zouaoui, F., and J. R. Wilson. 2004. “Accounting for input-model and input-parameter uncertainties insimulation”. IIE Transactions 36 (11): 1135–1151.
AUTHOR BIOGRAPHIESHENRY LAM is an Assistant Professor in the Department of Mathematics and Statistics at Boston University. He graduated from Harvard University with a Ph.D. degree in statistics in 2011. His research interests
lie in applied probability and Monte Carlo methods with applications in queueing, operations management
and insurance modeling. His email is khlam@bu.edu and his web page is at http://math.bu.edu/people/khlam/.
SOUMYADIP GHOSH is a Research Staff Member in the Business Analytics and Mathematical Sciences
Department at the IBM T.J. Watson Research Center. His current research interests lie in simulation based
optimization techniques for stochastic optimization problems, with a focus on applications in Energy andPower systems and supply chain management. His email is ghoshs@us.ibm.com and his web page is athttps://researcher.ibm.com/researcher/view.php?person=us-ghoshs.