Proceedings of the 2017 Winter Simulation Conference
W. K. V. Chan, A. D’Ambrogio, G. Zacharewicz, N. Mustafee, G. Wainer, and E. Page, eds.

COMPUTING WORST-CASE EXPECTATIONS GIVEN MARGINALS VIA SIMULATIONJose Blanchet
Fei HeDepartment of Management Science and Engineering
Department of IEOR
Stanford University
Columbia University
475 Via Ortega 3rd Floor
500 W 120th St
Stanford, CA 94305, USA
New York, NY 10027, USA
Henry Lam
Department of IEOR
Columbia University
500 W 120th St
New York, NY 10027, USA

ABSTRACT
We study a direct Monte-Carlo-based approach for computing the worst-case expectation of two multidimensional random variables given a specification of their marginal distributions. This problem is motivatedby several applications in risk quantification and statistics. We show that if one of the random variables
takes finitely many values, a direct Monte Carlo approach allows to evaluate such worst case expectation
with O(n−1/2 ) convergence rate as the number of Monte Carlo samples, n, increases to infinity.
1

INTRODUCTION

We focus on the problem of computing lower and upper bounds among any dependence structure for afunction of two random vectors whose marginal distributions are assumed to be known. This problem ismotivated from several applications in risk quantification and statistics. Before discussing its applications,let us first describe it precisely.
Suppose that X ∈ Rd follows distribution µ and Y ∈ Rl follows distribution ν. We define Π (µ, ν) to
be the set of joint distributions π in Rd×l such that the marginal of the first d entries coincides with µ andthe marginal of the last l entries coincides with ν. In other words,
for any probability measure π in Rd×ll(endowed with the Borel
σ -field), if we let πX (A) = π A × R for any Borel measurable set A ∈ Rd , andπY (B) = π Rd × B for any Borel measurable set B ∈ Rl , then π ∈ Π (µ, ν) if and only if πX = µ and
πY = ν. We are interested in the quantity (focusing on minimization)
V = min{Eπ [c (X,Y )] : π ∈ Π (µ, ν)}(1)where c(·, ·) ∈ R is some cost function. Formulation (1) is well-defined as the class Π (µ, ν) is non-empty,
because the product measure π = µ × ν belongs to Π (µ, ν).
In operations research contexts, problem (1) arises as a means to obtain bounds for performancemeasures in situations where dependence information is ambiguous. Such situations occur because, inpractice, accurately estimating the marginal distributions of random variables is often relatively easy, e.g.,
by goodness-of-fit against well-chosen parametric distributions. They also occur in scenarios where datafrom different stochastic sources are collected independently (i.e., rather than in pairs), in which case nodependence information between these sources can be inferred. Indeed, special (i.e., discrete) cases of (1)have been analyzed in the distributionally robust optimization literature (e.g., Doan et al. 2015). Variantsof (1) to risk measures have also been studied, regarding both algorithmic approaches (e.g., Rüschendorf

Blanchet, He and Lam
1983, Embrechts et al. 2013) and sharp bounds over specific geometric classes of marginals (e.g., Wang
and Wang 2011, Puccetti 2013, Puccetti and Rüschendorf 2013).
In statistics and machine learning contexts, the value of (1) is the Wasserstein distance (of order 1)
between X and Y when c(·, ·) is taken as a metric. The optimization can be viewed as the classical
Kantorovich relaxation to Monge’s problem in optimal transport (e.g., Rachev and Rüschendorf 1998,
Villani 2008), where solutions based on differential properties have been extensively studied. Wasserstein
distance is of central importance in probabilistic analysis (e.g., quantifying model discrepancies in Bayesiansettings (Minsker et al. 2014) and convergence rates of ergodic processes (Boissard and Le Gouic 2014),
among many others). The estimation of the distance itself is also suggested as a tool for statistical inference,
including the use in goodness-of-fit tests (Del Barrio et al. 1999, Del Barrio et al. 2005) and in applications
such as image recognition (Sommerfeld and Munk 2016). It has also been used to quantify model uncertainty
in stochastic optimization problems (e.g., Esfahani and Kuhn 2015, Blanchet and Kang 2016, Blanchet
and Murthy 2016, Gao and Kleywegt 2016) and in the application of distributionally robust optimizationin machine learning settings (Blanchet et al. 2016). As such, there have been growing studies on theconvergence behaviors of its empirical estimation. Central limit theorems (CLTs) on the empirical estimationof (1), based on representations using quantile functions, have been investigated in the one-dimensionalcase (e.g., Bobkov and Ledoux 2014, Del Barrio et al. 1999). More generally, concentration bounds havebeen studied in the line of work including Horowitz and Karandikar (1994), Bolley et al. (2007), Boissard(2011), Sriperumbudur et al. (2012), Trillos and Slepčev (2014) and Fournier and Guillin (2015), so dolaws of large numbers in some special cases (e.g., Dobrić and Yukich 1995).
Since classical methods for solving (1), based for instance on Euler-Lagrange equations, may not yieldstraightforward computational schemes in general, we resort to Monte Carlo for an easy-to-implement
approximation. Our contribution is precisely to quantify the rate of convergence of such Monte Carlo
schemes. Our results also add to the literature of empirical Wasserstein estimation when these Monte Carlosamples are viewed as data. We focus on the setting where one of the marginals, say Y , is a finite-support
distribution, and another, say X, is a multi-dimensional distribution that can be continuous. To approximate
V , we consider the drawn samples from the continuous variable X, and replace the infinite-dimensional
linear program (LP) in (1) by its sampled counterpart, which can be solved by standard LP solvers.
Our main result shows that the error of our procedure is O(n−1/2 ) where n is the sample size, independentof the dimension d or l. We also identify the limiting distribution in the associated CLT. The closest workto our results, as far as we know, is the recent work of Sommerfeld and Munk (2016), who derive a CLT
when both marginal distributions are finitely discrete. Our result here can be viewed as a generalization totheirs when one of the distributions is continuous. We remark that our obtained rate differs from the typicalrate of O(n−1/d ) in high-dimensional empirical Wasserstein estimation where d ≥ 3 is the dimension of themarginal distributions. As we will see, the finite-support property of one of the marginals plays a crucialrole in applying classical results in sample average approximation (SAA) that maintain the standard Monte
Carlo rate in our scheme.
In the rest of this paper, we will first describe our algorithm, followed by our main results on theconvergence analysis.
2

ALGORITHMIC DESCRIPTION

Suppose that the distribution ν for Y has finite support {y1 , ...ym } ⊂ Rl . Supposing that X can be simulated,
we sample n i.i.d. observations X1 , . . . , Xn from µ, and approximate V by
Vn = min{Eπ [c (X,Y )] : π ∈ Π (µn , ν)}
where µn is the empirical distribution of X constructed from the Xi ’s, i.e.,µn (A) =1 n∑ I(Xi ∈ A)
n i=1

(2)

Blanchet, He and Lam
for any Borel measurable A.
Note that (2) is a finite-dimensional LP, which can be written more explicitly as
min∑ni=1 ∑mj=1 c(Xi , y j )pi j
subject to ∑mj=1 pi j = 1n ∀i = 1, . . . , n∑ni=1 pi j = ν{y j } ∀ j = 1, . . . , m
pi j ≥ 0 ∀i = 1, . . . , n, j = 1, . . . , m(3)where the decision variables pi j represent the probability masses on (Xi , y j ), and ν{y j } denotes the mass on y j
under ν. Problem (3) is an assignment problem, which is a special type of minimum cost problem and can besolved by, e.g., successive shortest path algorithms in polynomial time of order O(n2 m+n(n+m) log(n+m))(see, e.g., R.K.Ahuja et al. 2000 pp. 471, 500).
3

CONVERGENCE ANALYSISOur main result is a convergence analysis on Vn to V . We impose the assumptions:
Assumption 1 For each y j , c(., y j ) is non-negative and lower semicontinuous.
Assumption 2 Suppose that ν has finite support {y1 , ..., ym } ⊂ Rl . We haveEµ [c(X, y j )2 ] < ∞, ∀ j = 1, . . . , m.
DenoteV0 =max

β1 ,...,βm ∈R

h
Eµ

m

min {c (X, y j ) − β j } + ∑ β j ν{y j }j=1,...,m

i

(4)

j=1

which is the dual problem of (1) (see Lemma 1 for an explanation in the special case of finite-dimensionalsettings). Under Assumptions 1 and 2, strong duality (known as the Kantorovich duality) holds and V 0 = V ;see, e.g., Theorem 5.10 in Villani (2008).
In order to state our main result, we need to introduce a Gaussian random field G(·) : Rm → R withcovariance structure given byCov(G(β ), G(β 0 )) = Cov min c(X, y j ) − β j , min c(X, y j ) − β j0
j=1,...,m

j=1,...,mfor any β = (β j )mj=1 and β 0 = (β j0 )mj=1 . Our main result is the following.
p

Theorem 1 Under Assumption 2, Vn → V 0 as n → ∞. Moreover,
n1/2 Vn −V 0 ⇒ G∗
as n → ∞, where

G∗ =

max

G (β ) .

β =(β1 ,...,βm )∈S

Here S is the set of all optimal solutions β = (β j )mj=1 ∈ Rm for the convex optimization problem
max

β1 ,...,βm ∈R∑mj=1 β j =0

h
Eµ

i
m
min {c (X, y j ) − β j } + ∑ β j ν{y j } .

j=1,...,m

(5)

j=1

Remark 1 The significance of this result is that one can approximate worst-case expectations by samplingwith a rate of convergence (as measured by the sample size of the continuous distribution) of order O(n−1/2 ).
As we mentioned earlier, this might be somewhat surprising given that standard empirical estimators for
Wasserstein distances exhibit a degradation which becomes quite drastic in high dimensions.

Blanchet, He and Lam
3.1 Proof of Theorem 1
We first note that adding a constant to β j in the objective function of the dual does not change the objective
value. To remove this ambiguity we inroduce the next result.
Lemma 1 Define
)
(
n
m

1
(6)
Vbn :=
max∑ min c(Xi , y j ) − β j + ∑ β j ν{y j } .
β j ∈R, j=1,...,m. n i=1 j=1,...,mj=1∑mj=1 β j =0

We have Vn = Vbn .
Proof.

The dual formulation of Vn , depicted as the LP (3), is given by
m
1 n
max
n ∑i=1 αi + ∑ j=1 β j ν{y j }subject to αi + β j ≤ c(Xi , y j ) ∀i = 1, . . . , n, j = 1, . . . , m(7)m
where (αi )mi=1 , (β j ) j=1 are the dual variables. Note that the constraint in (7) can be written as αi ≤min j=1,...,m {c(Xi , y j ) − β j } ∀i = 1, . . . , n, which implies that (7) is equivalent to
(
)
m

1 n
max
(8)
∑ min c(Xi , y j ) − β j + ∑ β j ν{y j }β j ∈R, j=1,...,m n i=1 j=1,...,mj=1

Since shifting any (β j )mj=1 to (β j + λ )mj=1 by an arbitrary constant λ does not affect the objective valueof (8), we can always set λ = − m1 ∑mj=1 β j to enforce the constraint ∑mj=1 β j = 0, so that (8) is equal to (6).
Finally, since (3) is feasible by choosing an independent distribution, strong duality holds. We therefore
conclude the lemma.
Next we show that Vbn can be further reduced to a problem with compact feasible region, which willsubsequently facilitate the invocation of classical results in SAA:
Proposition 2 Define
)
(
n
m

1
cb :=
(9)
V
max
∑ min c(Xi , y j ) − β j + ∑ β j ν{y j } .
nβ j ∈R,|β j |≤b, j=1,...,m n i=1 j=1,...,mj=1∑mj=1 β j =0

There exists some large enough constant b > 0 such thatcb
Vn = Vn(10)eventually, i.e., holds for any n > N for some N < ∞ almost surely.
Proof.

By Lemma 1, we have
Vn = Vbn
(
= max

(
max

β j ∈R,|β j |≤b, j=1,...,m∑mj=1 β j =0

)
m

1 n
min c(Xi , y j ) − β j + ∑ β j ν{y j } ,∑ j=1,...,m
n i=1
j=1
(

max

β j ∈R, j=1,...,m,|β j |>b for some j∑mj=1 β j =0

m

1 n
min
c(X
,
y
)
−
β
+
i
j
j
∑ j=1,...,m∑ β j ν{y j }
n i=1
j=1

))
.

(11)

Blanchet, He and Lam
cb by our definition (9). We will show that there exists a
Note that the first term inside the outer max is Vndeterministic b > 0 such that the first term dominates the second term eventually, which will then conclude
the proposition.
To this end, consider the second term in (11)
(
)
m

1 n
max
∑ min c(Xi , y j ) − β j + ∑ β j ν{y j }β j ∈R, j=1,...,m,|β j |>b for some j n i=1 j=1,...,mj=1∑mj=1 β j =0

m


≤

max

β j ∈R, j=1,...,m,|β j |>b for some j∑mj=1 β j =0min



j=1,...,m− β j + ∑ β j ν{y j } +j=1

1 n
max c(Xi , y j ).
∑ j=1,...,m
n i=1

(12)

We analyze
m


max

β j ∈R, j=1,...,m,|β j |>b for some j∑mj=1 β j =0min

j=1,...,m− β j + ∑ β j ν{y j } .

(13)

j=1

Denote M = max j=1,...,m |β j |, so that M > b for any β inside the feasible region. There must exist either aβ j∗ = M or β j∗ = −M. In the first case, we havem
min − β j + ∑ β j ν{y j }max
β j ∈R, j=1,...,m,|β j |>b for some j∑mj=1 β j =0

j=1,...,mj=1

∑mj=1 β j ν{y j } max≤ −M + subject to β j ≤ M ∀ j = 1, . . . , m∑mj=1 β j = 0∑mj=1 β j ν{y j } max= −M + M × subject to β j ≤ 1 ∀ j = 1, . . . , m∑mj=1 β j = 0(14)where the last equality follows by a change of variable from β j to β j /M in the optimization. Note that the
optimal value ofmax∑mj=1 β j ν{y j }subject to β j ≤ 1 ∀ j = 1, . . . , m∑mj=1 β j = 0
is strictly less than 1. To see this, observe that the optimal value is at most 1 by using the first constraint.
The value of exactly 1 is attained under the first constraint by the unique solution β j = 1, j = 1, . . . , m,which is ruled out because it would violate the second constraint. With this claim, we conclude that (14)is equal to θ M for some θ < 0, which is bounded from above by θ b.
In the second case, we have β j∗ = −M. Let j˜∗ = argmax j=1,...,m {β j }. By the constraint ∑mj=1 β j = 0 in(13), we must have β j˜∗ ≥ M/(m − 1). Therefore, applying our argument for the first case gives that (13)is bounded from above by θ M/(m − 1) ≤ θ b/(m − 1) for the same θ < 0 chosen before.
Therefore, in either case (13) is bounded from above by θ b/(m − 1). Note that the first term inside
cb , satisfies V
cb ≥ (1/n) ∑n min
the outer max in (11), namely V
j=1,...,m c(Xi , y j ) by plugging in the feasiblenni=1solution given by β j = 0, j = 1, . . . , m. Thus, with the law of large numbers, by choosing b > 0 large enoughsuch that




θb
+ Eµ max c(X, y j ) < Eµ min c(X, y j )(15)
j=1,...,m
j=1,...,mm−1

Blanchet, He and Lam
the first term dominates the second term inside the outer max in (11) as n → ∞ almost surely.
We are now ready to prove Theorem 1:
Note that the function

Proof of Theorem 1.

F(X, β ) := min

j=1,...,mm

c(X, y j ) − β j + ∑ β j ν{y j }

(16)

j=1

on β = (β j )mj=1 ∈ Rm is Lipschitz continuous in the sense that
|F(X, β ) − F(X, β 0 )| ≤ (1 + kνk)kβ − β 0 k
where k · k denotes the L2 -norm, and ν is interpreted as a vector (ν{y j })mj=1 . This follows sinceminj=1,...,m
c(X, y j ) − β j − min c(X, y j ) − β j0
j=1,...,m≤ kβ − β 0 k∞

and
m

m

j=1

j=1

∑ β j ν{y j } − ∑ β j0 ν{y j } ≤ kνkkβ − β 0 k

by the Cauchy-Schwarz inequality. Since the set B := {β ∈ Rm : ∑mj=1 β j = 0, |β j | ≤ b, ∀ j = 1, . . . , m} iscompact and Eµ [F(X, β )2 ] < ∞ by Assumption 2, by using Theorem 5.7 in Shapiro et al. (2009), we have

and

p
cb →
Vb
V
n

(17)

√ cb
n(Vn −V b ) ⇒ G∗,b

(18)

where
Vb =

max

h

β j ∈R,|β j |≤b, j=1,...,m∑mj=1 β j =0

and

Eµ

G∗,b =im
min {c (X, y j ) − β j } + ∑ β j ν{y j }j=1,...,mmax

β =(β1 ,...,βm )∈Sb

(19)

j=1

G (β )

with Sb denoting the set of optimal solutions for (19) and G(·) is defined as in Theorem 1 but restricted
to the domain B.
√ cbp
By Proposition 2, we have n(Vn −Vn ) → 0 as n → ∞. Thus, together with (17), we havep

Vn → V b
and together with (18), we haveby Slutsky’s Theorem.

√n(Vn −V b ) ⇒ G∗,bBlanchet, He and Lam
To conclude the theorem, we show that V b = V 0 , and Sb = S so that G∗,b = G∗ . By using essentiallythe same argument as for Proposition 2 (with the empirical expectation replaced by Eµ [·]) and choosing
the same b as in (15), we have
V0 =
=

i
m
min {c (X, y j ) − β j } + ∑ β j ν{y j }hmax

β j ∈R, j=1,...,mEµj=1,...,mhmax

β j ∈R, j=1,...,m∑mj=1 β j =0

Eµ

j=1
m

i
min {c (X, y j ) − β j } + ∑ β j ν{y j }j=1,...,mj=1

mby shifting any (β j )mj=1 to (β j − (1/m) ∑ βk )mj=1 which does not affect the objective value andk=1menforces the constraint ∑ β j = 0
j=1

(
= max

i
m
min {c (X, y j ) − β j } + ∑ β j ν{y j } ,hmax

β j ∈R,|β j |≤b, j=1,...,m∑mj=1 β j =0Eµ

j=1,...,mhmax

β j ∈R, j=1,...,m,|β j |>b for some j∑mj=1 β j =0

Eµ

j=1

)
i
min {c (X, y j ) − β j } + ∑ β j ν{y j }mj=1,...,m

j=1

where
Vb
=

>

h

max

β j ∈R,|β j |≤b, j=1,...,m∑mj=1 β j =0

max

Eµ

m

min {c (X, y j ) − β j } + ∑ β j ν{y j }j=1,...,mβ j ∈R, j=1,...,m,|β j |>b for some j∑mj=1 β j =0

j=1

h
Eµ

i

m

min {c (X, y j ) − β j } + ∑ β j ν{y j }j=1,...,mi

j=1

so that V 0 = V b and Sb = S.
4

ADDITIONAL DISCUSSION AND EXTENSIONS

Finally, we briefly discuss the challenge in generalizing our procedure to the case when both X and Y
are continuous. Here, one may attempt to sample both variables (assuming both can be simulated) andformulate a sampled program like (2) or (3). However, the analog of its reformulation in (6) and (9) will
have a growing number of variables β j and an analogous limit in (5) that involves an infinite-dimensionalvariable, which challenges the use of standard SAA machinery. In fact, consider a special example whereX,Y ∼ U[0, 1]d and c(x, y) = kx − yk. In this case, (1) corresponds to the Wasserstein distance (of order 1)
between X and Y , which is of course 0. It is known that sampling X and keeping Y continuous will give, ford ≥ 3, an expected optimal value of (2) that is of order n−1/d , i.e., C1 n−1/d ≤ EVn ≤ C2 n−1/d for all n forsome C1 ,C2 > 0 (e.g., Problem 5.11 in van Handel 2014). Thus, the convergence rate deteriorates with thedimension and the standard Monte Carlo rate O(n−1/2 ) cannot be maintained without assuming additional
structure or infomation available to the modeler on the primal problem. It is of interest to investigate
reasonable assumptions which are useful in applications and which would mitigate such rate-of-convergence
deterioration.

Blanchet, He and Lam
ACKNOWLEDGMENTS
We gratefully acknowledge support from the National Science Foundation (NSF) under grants CMMI1542020, CMMI-1523453 and CAREER CMMI-1653339. Additional support from NSF grants CMMI1436700, DMS-1320550, and CMMI-1538217 is gratefully acknowledged.
REFERENCES
Blanchet, J., and Y. Kang. 2016. “Sample Out-Of-Sample Inference Based on Wasserstein Distance”. arXiv
preprint arXiv:1605.01340.
Blanchet, J., Y. Kang, and K. Murthy. 2016. “Robust Wasserstein Profile Inference and Applications to
Machine Learning”. arXiv preprint arXiv:1610.05627v2.
Blanchet, J., and K. Murthy. 2016. “Quantifying distributional model risk via optimal transport”.
Bobkov, S., and M. Ledoux. 2014. “One-dimensional empirical measures, order statistics and Kantorovich
transport distances”. preprint.
Boissard, E. 2011. “Simple bounds for the convergence of empirical and occupation measures in 1-Wasserstein
distance”. Electronic Journal of Probability 16:2296–2333.
Boissard, E., and T. Le Gouic. 2014. “On the mean speed of convergence of empirical and occupation
measures in Wasserstein distance”. In Annales de l’Institut Henri Poincaré, Probabilités et Statistiques,
Volume 50, 539–563. Institut Henri Poincaré.
Bolley, F., A. Guillin, and C. Villani. 2007. “Quantitative concentration inequalities for empirical measures
on non-compact spaces”. Probability Theory and Related Fields 137 (3): 541–593.
Del Barrio, E., J. A. Cuesta-Albertos, and C. Matrán. 1999. “Tests of goodness of fit based on the
L 2-Wasserstein distance”. The Annals of Statistics 27 (4): 1230–1239.
Del Barrio, E., E. Giné, and F. Utzet. 2005. “Asymptotics for L2 functionals of the empirical quantile
process, with applications to tests of fit based on weighted Wasserstein distances”. Bernoulli 11 (1):
131–189.
Doan, X. V., X. Li, and K. Natarajan. 2015. “Robustness to dependency in portfolio optimization using
overlapping marginals”. Operations Research 63 (6): 1468–1488.
Dobrić, V., and J. E. Yukich. 1995. “Asymptotics for transportation cost in high dimensions”. Journal of
Theoretical Probability 8 (1): 97–118.
Embrechts, P., G. Puccetti, and L. Rüschendorf. 2013. “Model uncertainty and VaR aggregation”. Journalof Banking & Finance 37 (8): 2750–2764.
Esfahani, P. M., and D. Kuhn. 2015. “Data-driven distributionally robust optimization using the Wasserstein
metric: Performance guarantees and tractable reformulations”. arXiv preprint arXiv:1505.05116.
Fournier, N., and A. Guillin. 2015. “On the rate of convergence in Wasserstein distance of the empiricalmeasure”. Probability Theory and Related Fields 162 (3-4): 707–738.
Gao, R., and A. J. Kleywegt. 2016. “Distributionally robust stochastic optimization with Wasserstein
distance”. arXiv preprint arXiv:1604.02199.
Horowitz, J., and R. L. Karandikar. 1994. “Mean rates of convergence of empirical measures in the
Wasserstein metric”. Journal of Computational and Applied Mathematics 55 (3): 261–273.
Minsker, S., S. Srivastava, L. Lin, and D. Dunson. 2014. “Scalable and robust Bayesian inference via themedian posterior”. In International Conference on Machine Learning, 1656–1664.
Puccetti, G. 2013. “Sharp bounds on the expected shortfall for a sum of dependent random variables”.
Statistics & Probability Letters 83 (4): 1227–1232.
Puccetti, G., and L. Rüschendorf. 2013. “Sharp bounds for sums of dependent risks”. Journal of Applied
Probability 50 (01): 42–53.
Rachev, S. T., and L. Rüschendorf. 1998. Mass Transportation Problems: Volume I: Theory, Volume 1.
Springer Science & Business Media.
R.K.Ahuja, T.L.Magnanti, and J.B.Orlin. 2000. “Network Flows”. Prentice Hall, Inc..

Blanchet, He and Lam
Rüschendorf, L. 1983. “Solution of a statistical optimization problem by rearrangement methods”. Metrika 30(1): 55–61.
Shapiro, A., D. Dentcheva, and A. Ruszczyński. 2009. Lectures on stochastic programming: modeling andtheory. SIAM.
Sommerfeld, M., and A. Munk. 2016. “Inference for Empirical Wasserstein Distances on Finite Spaces”.
arXiv preprint arXiv:1610.03287.
Sriperumbudur, B. K., K. Fukumizu, A. Gretton, B. Schölkopf, and G. R. Lanckriet. 2012. “On the empiricalestimation of integral probability metrics”. Electronic Journal of Statistics 6:1550–1599.
Trillos, N. G., and D. Slepčev. 2014. “On the rate of convergence of empirical measures in ∞-transportation
distance”. arXiv preprint arXiv:1407.1157.
van Handel, R. 2014. “Probability in high dimension”. Technical report, DTIC Document.
Villani, C. 2008. Optimal transport: old and new, Volume 338. Springer Science & Business Media.
Wang, B., and R. Wang. 2011. “The complete mixability and convex minimization problems with monotone
marginal densities”. Journal of Multivariate Analysis 102 (10): 1344–1360.
AUTHOR BIOGRAPHIESJOSE BLANCHET is a faculty member in the Department of Management Science and Engineering at
Stanford University and of IEOR and Statistics at Columbia University. He holds a Ph.D. in Management
Science and Engineering from Stanford University. From 2004 to 2008 he was a faculty member in the
Department of Statistics at Harvard University. Jose is a recipient of the 2009 Best Publication Award
given by the INFORMS Applied Probability Society and of the 2010 Erlang Prize. He also received a
PECASE award given by NSF in 2010. He worked as an analyst in Protego Financial Advisors, a leading
investment bank in Mexico. He has research interests in applied probability and Monte Carlo methods.
His email is jose.blanchet@stanford.edu.
FEI HE is a PhD candidate in the Department of Industrial Engineering and Operations Research atColumbia University. His research focuses on stochastic simulation, mathematical finance and risk management. His email address is fh2293@columbia.edu.
HENRY LAM is an Associate Professor in the Department of Industrial Engineering and Operations
Research at Columbia University. He received his Ph.D. degree in statistics at Harvard University, and wasa faculty member in the Department of Mathematics and Statistics at Boston University and the Department
of Industrial and Operations Engineering at the University of Michigan before joining Columbia. His
research focuses on stochastic simulation, risk analysis, and simulation optimization. His email address is
khl2114@columbia.edu.