Proceedings of the 2018 Winter Simulation Conference
M. Rabe, A. A. Juan, N. Mustafee, A. Skoogh, S. Jain, and B. Johansson, eds.

SAMPLING UNCERTAIN CONSTRAINTS UNDER PARAMETRIC DISTRIBUTIONSHenry Lam
Fengpei Li
Department of Industrial Engineering and Operations Research
Columbia University
500 W. 120th Street
New York, NY 10027, USA

ABSTRACT
We consider optimization problems with uncertain constraints that need to be satisfied probabilistically.
When data are available, a common method to obtain feasible solutions for such problems is to impose
sampled constraints, following the so-called scenario generation (SG) approach. However, when the data
size is small, the sampled constraints may not support a guarantee on the feasibility of the obtained solution.
This paper studies how to leverage parametric information and the power of Monte Carlo simulation to
obtain feasible solutions even when the data are not sufficient to support the use of SG. Our approachmakes use of a distributionally robust optimization formulation that informs the Monte Carlo sample size
needed to achieve our guarantee.
1

INTRODUCTION

We consider optimization problems in the formminx∈X ⊆RdcT x,subject to P(x ∈ Xξ ) ≥ 1 − ε,(1)where P is a probability measure governing the random variable ξ , and Xξ ⊂ X is a set depending on ξ .
Problem (1) enforces a solution x to satisfy x ∈ Xξ with high probability, namely at least 1 − ε. This problem
is often known as a chance-constrained or probabilistically constrained optimization (e.g., Prékopa 2003).
It provides a natural framework for decision-making under stochastic resource capacity or risk tolerance,and has been applied in various domains such as production planning (Murr and Prékopa 2000), inventory
management (Lejeune and Ruszczynski 2007), reservoir design (Prékopa and Szántai 1978; Prékopa et al.
1978), communications (Shi et al. 2015), and ranking and selection (Hong et al. 2015).
We focus on the situations where P is unknown, but some data, say ξ1 , . . . , ξn , are available. One
common approach to handle (1) in these situations is to use the so-called scenario generation (SG) or
constraint sampling. This replaces the unknown constraint in (1) with x ∈ Xξi , i = 1, . . . , n, namely, by
considering
min
cT x,
x∈X ⊆Rd(2)
subject to x ∈ Xξi , i = 1, . . . , n.
Note that the chance-constrained problem (1) is generally difficult to solve, even when the set Xξ isconvex and tractable for any given ξ (Prékopa 2003). Thus, in the latter case, the sampled problem (2)offers an additional benefit in approximating the intractable problem with a more tractable one.
Our goal is to find a good feasible solution for (1) in the described data-driven context above. Note
that, because of the statistical noise from the data, we can at best find a solution that is feasible with a highLam and Li
confidence. Define V (x, P) = P(x ∈/ Xξ ) to be the violation probability of a solution x for the conditionx ∈ Xξ under probability measure P. For any x̂ that is obtained from finite data, we want to make sure that
Pdata (V (x̂, P) ≤ ε) ≥ 1 − α,(3)for a given confidence level 1 − α (e.g., α = 5%), where Pdata denotes the measure generating the dataξi , i = 1, . . . , n (and hence x̂).
As n → ∞, one would expect that the space of ξ is sufficiently populated by the sample and V (x̂, P) → 0.
The question is then how many observations are enough to see this behavior. The seminal work of Campi
and Garatti (2008) shows a tight bound on the number of observations, or sampled constraints,
needed

N i
ε
(1
−
ε)N−i
to guarantee (3), for given values of ε and α. In particular, define γ(N, ε, d) = ∑d−1i=0 i
for positive integers N ≥ d ≥ 1 and 0 < ε < 1. Under the convexity of Xξ and some additional mild
assumptions, Campi and Garatti (2008) proves that an optimal solution obtained from solving (2) satisfies
Pdata (V (x̂, P) > ε) ≤ γ(n, ε, d).

(4)Thus, if we can find n such that γ(n, ε, d) ≤ α, then we achieve (3).
In this paper, we focus on the small-sample situation, so that our data size n is not enough to support
γ(n, ε, d) ≤ α. Note that, in using this result from Campi and Garatti (2008), n can be seen to be linear in
the decision dimension d, and so for high-dimensional problems this small-sample situation can happen
frequently. This dimensional dependence also appears in other sample size bounds (e.g., De Farias and
Van Roy 2004; Luedtke and Ahmed 2008). To overcome this challenge, several recent methods have beensuggested, such as the use of support rank and solution-dependent support constraints (Schildbach et al.
2013; Campi and Garatti 2018), regularization (Campi and Carè 2013), and sequential approaches (Carèet al. 2014; Calafiore et al. 2011; Chamanbaz et al. 2016; Calafiore 2017). They aim to alleviate thedependence on d, and thus substantially extend the scope of applicability of SG.
Our main contribution in this paper is to study a different path in obtaining guarantee (3) in small-sample
situation, in the settings where P is assumed a parametric structure. The unknown and estimable quantity is
the set of parameters in P. We also assume that one can simulate from the parametric model P using Monte
Carlo (an assumption applied for all common parametric models). We will see how such a capability, which
can be viewed as a way of generating additional synthetic data, can be combined with the relation (4) toobtain a scheme applicable when n does not support using (4) directly. Unlike some other techniques that
reduce the data size needed for SG, our procedure resembles closely the standard SG (2). The differences
are in the algorithmic parameters and the distribution we sample from in the Monte Carlo scheme.
Our derivation relies on casting the uncertain-parameter problem as a distributionally robust optimization(DRO) (Delage and Ye 2010; Wiesemann et al. 2014). This approach considers the worst-case situation
among all parameters that lie in a so-called uncertainty set or ambiguity set. In the chance constraint
framework, this entails replacing the chance constraint with unknown distribution with a worst-case chanceconstraint within the uncertainty set (Hanasusanto et al. 2015; Zymler et al. 2013; Hanasusanto et al.
2017; Li et al. ; Jiang and Guan 2016; Zhang et al. 2016). The DRO approach has been used in stochastic
simulation, bearing names such as the robust Monte Carlo (Hu et al. 2012; Glasserman and Xu 2014;
Lam 2016; Hu and Hong 2015; Lam 2018; Ghosh and Lam 2018). Our procedure will rely on a suitable
DRO formulation using statistical distances (Petersen et al. 2000; Ben-Tal and Nemirovski 2000; Limet al. 2006; Love and Bayraksan 2015), and the particular SG we use has a similar favor as the robust
Monte Carlo considered in Hu et al. (2012) and Glasserman and Xu (2014), as we also utilize a changeof-measure argument in arriving at our scheme. Lastly, Erdoğan and Iyengar (2006) considers an SG for
robust chance-constrained problems that is related to our proposal. Erdoğan and Iyengar (2006) considers
uncertainty set based on the Prohorov distance and derives bounds for the required sample size. On theother hand, our approach uses the class of φ - or f -divergences, which is readily estimatable especially inthe parametric setting. Moreover, our focus is on utilizing the convexity-based sample size estimate inCampi and Garatti (2008), in contrast to the Vapnik-Chervonenkis dimension used in Erdoğan and Iyengar

Lam and Li(2006). We will study the required Monte Carlo size in relation to the data size, and the choice of theMonte Carlo distribution, which are quite different from Erdoğan and Iyengar (2006).
2

OUTLINE OF THE METHOD

Recall that we are interested in finding a solution x̂ such that (3) holds, in the case that P is observable
only through data. Suppose that P ∈ P, the class of all possible probability distributions for ξ (which weshall exemplify later). Suppose that given our data, we can find an uncertainty set Udata ∈ P such that
Pdata (P ∈ Udata ) ≥ 1 − α.

(5)We then proceed to consider a distributionally robust chance-constrained problemminx∈X ⊆Rd

subject to

cT x,
(6)

inf Q(x ∈ Xξ ) ≥ 1 − ε,Q∈Udatawhere Q is the probability measure for ξ . If we can find a solution x̂ feasible for (6), then this x̂ is also
feasible for (1) with confidence at least 1 − α. This is because if P ∈ Udata , then for any x feasible for (6),P(x ∈ Xξ ) ≥ infQ∈Udata Q(x ∈ Xξ ) ≥ 1 − ε, and therefore Pdata (V (x̂, P) ≤ ε) ≥ Pdata (P ∈ Udata ) ≥ 1 − α.
We now consider a Monte Carlo scheme from some “baseline” distribution P0 (which can depend on
the data), to obtain a solution that is feasible for (6) with a confidence of, say, 1 − β . This is achievable byusing a certain Monte Carlo size N. To obtain this number, we find a bound on sup V (x̂(P0 ), Q), where
Q∈Udatax̂(P0 ) is obtained from solving (2) using ξi ’s generated from P0 . Call this bound M(P0 , Udata ,V (x̂(P0 ), P0 )).
In other words, it satisfies
sup V (x̂(P0 ), Q) ≤ M(P0 , Udata ,V (x̂(P0 ), P0 )).

(7)

Q∈UdataMoreover, suppose we also have that M(P0 , Udata , v) is non-decreasing in v > 0. We then find a δ > 0such that
M(P0 , Udata , δ ) ≤ ε.
(8)Then, by using the result in Campi and Garatti (2008) discussed in the introduction, we know that
PMC,0 (V (x̂(P0 ), P0 ) > δ ) ≤ γ(N, δ , d) where PMC,0 is the measure generating N Monte Carlo samples from
P0 to obtain x̂(P0 ), which holds independent of the choice of P0 . We find N such thatγ(N, δ , d) ≤ β .

(9)This choice of N then gives us, with confidence 1−β , thatV (x̂(P0 ), P0 ) ≤ δ and hence sup V (x̂(P0 ), Q) ≤Q∈Udataε by (7), (8) and the monotonicity of M. This in turn leads to the conclusion that x̂(P0 ) is feasible for (6)with confidence at least 1 − β .
Thus, overall, if we choose Udata to satisfy (5) and N to satisfy (8) and (9), our obtained solution x̂(P0 )is feasible for (1) with confidence at least 1 − α − β . Note that our argument relies crucially on generating
a bound M that translates the ambiguous violation probability sup V (x̂(P0 ), Q) into a quantity depending
Q∈Udatainstead on the baseline violation probability V (x̂(P0 ), P0 ). Next section will address how to do so.
2.1 Bounding Ambiguous Violation Probability
We study the bound M(P0 , Udata ,V (x̂(P0 ), P0 )). In fact, we will consider a more general result. Consider
any measurable set A ⊂ Y and set U ⊂ P. We will find an M such thatsup Q(ξ ∈ A ) ≤ M(P0 , U , P0 (ξ ∈ A )).
Q∈ULam and Li
For two probability measures P1 and P2 that are both dominated by a common measure ν on Y , withdP2
2
1
Radon-Nikodym derivatives dPdν and dν respectively, we define the χ -distance between P1 and P2 asχ 2 (P1 , P2 ) =

Z

dP1 2
2
( dP
dν − dν )
dP1
dν

Y

Z

ν(dy) =Y(

dP2 dP1
dP1/− 1)2
ν(dy).
dν dνdνWhen P2 is absolutely continuous with respect to P1 , we have in particularχ 2 (P1 , P2 ) =

Z
Y

(

dP2
− 1)2 P1 (dy) =dP1

dP2
P2 (dy) − 1.
Y dP1

Z

Suppose that Q is absolutely continuous with respect to P0 . We have
sup Q(ξ ∈ A ) = P0 (ξ ∈ A ) + sup Q(ξ ∈ A ) − P0 (ξ ∈ A )
Q∈U

Q∈U

= P0 (ξ ∈ A ) + supZ

1{y ∈ A }Q(dy) −Z

1{y ∈ A }P0 (dy)

Q∈U




dQ
= P0 (ξ ∈ A ) + sup 1{y ∈ A }− 1 P0 (dy)
dP0
Q∈U
Z
1/2
1/2  Z 2
dQ≤ P0 (ξ ∈ A ) + sup− 1 P0 (dy)1{y ∈ A }P0 (dy)
Y
Y dP0
Q∈U
Z

= P0 (ξ ∈ A ) + P0 (ξ ∈ A )1/2 · ( sup χ 2 (P0 , Q))1/2 ,

(10)

Q∈U

where the inequality follows from the Cauchy-Schwarz inequality. Note that the bound (10) holds if A
and U are random sets, potentially dependent on each others, but independent from ξ generated from Q
or P0 above. Thus, by plugging in x̂(P0 ) ∈/ Xξ as ξ ∈ A and Udata as U , we havesup V (x̂(P0 ), Q) ≤ V (x̂(P0 ), P0 ) +V (x̂(P0 ), P0 )1/2 · ( sup χ 2 (P0 , Q))1/2 .
Q∈UdataQ∈UdataThus, we have identified
M(P0 , Udata , v) = v + v1/2 · ( sup χ 2 (P0 , Q))1/2 ,(11)Q∈Udatawhich is non-decreasing in v.
2.2 Choices of Baseline and Uncertainty Set
Our next step as outlined in the beginning of Section 2 is to find δ such that M(P0 , Udata , δ ) ≤ ε, andmoreover a Udata that satisfies (5). For these, we now make an assumption that P, the true distribution ofξ , is known to lie in a parametric family P = {Pθ }θ ∈Θ⊂R p indexed by θ , which has dimension p.
We make a convenient choice for P0 and Udata . Namely, we choose P0 to be Pθ̂ , where θ̂ is the
maximum likelihood estimator for θ from the data ξi , i = 1, . . . , n. Then we set
(
)
2
χ
1−α,p
Udata = Q ∈ P : χ 2 (Pθ̂ , Q) ≤
,
(12)
n
2
where χ1−α,pis the 1 − α-quantile of the χ 2 -distribution with degree of freedom p. By divergence-based
inference (e.g., Pardo 2005), we have that (12) satisfies

lim Pdata (P ∈ Udata ) = 1 − α.

n→∞Lam and Li
Moreover, (11) becomes
v+v

Thus, we choose δ such that δ + δ 1/2



1/2

2
χ1−α,pn2χ1−α,p

n

!1/2
.

1/2
≤ ε, or equivalently

v

2
u
u
2
2
2
χ1−α,pt χ1−α,pχ1−α,p
,
δ =ε+
− ε·
+
2n
n
4n2

(13)

and the Monte Carlo size N such that γ(N, δ , d) ≤ β . Using this N number of samples generated fromPθ̂ to construct (2) then guarantees that the obtained solution x̂(Pθ̂ ) is feasible for (1) with a confidence
1 − α − β + o(1) as n → ∞. We note that formula (13) appears in a related context in Proposition 2 in Tsenget al. (2016), but here we are motivated by the use of SG in very limited data situations and investigate
the joint construction of the uncertainty set and SG with overall statistical guarantees, which is different
from the work of Tseng et al. (2016).
Note that the χ 2 -distance we used above is one of many distances between probability distributions
that can be categorized under the framework of φ - or f -divergences (Vajda 1972). There are reasons tobelieve that χ 2 is close to the best under our discussed framework. First is that absolute continuity between
the true distribution P and the baseline distribution P0 to generate Monte Carlo samples is critical forour approach. Suppose they are not absolutely continuous. One could attempt to form an uncertainty set
Udata that contains both P and P0 with high probability (this can be done using, e.g., Wasserstein distance;
Esfahani and Kuhn 2015; Blanchet and Kang 2016; Gao and Kleywegt 2016). However, the difference
sup Q(ξ ∈ A ) − P0 (ξ ∈ A ) becomes more difficult to control as it can depend intricately on the set
Q∈Udata

A . On the other hand, one can instead use Udata that contains distributions absolutely continuous withrespect to P or vice versa. But if P0 lies outside this set, the difference sup Q(ξ ∈ A ) − P0 (ξ ∈ A ) is
Q∈Udataagain hard to control. Therefore, to apply our framework, choosing a set Udata that contains distributions
absolutely continuous with respect to P0 and also contains P facilitates the error control tremendously.
Next, many other φ -divergences can work under our framework, since one can find Udata in much thesame way using general divergence-based inference tools. This comes from the fact that these divergences
between two distributions indexed by θ1 and θ2 have the same expansion up to the second order (Nielsen
and Nock 2014). Then, as long as we can bound the difference sup Q(ξ ∈ A ) − P0 (ξ ∈ A ), the same
Q∈Udataargument to look for a suitable N will apply. However, it appears that χ 2 -distance can perform better
than other common candidates, including the Kullback-Leibler (KL) divergence, precisely because of thisdifference bound. In the KL case, Pinkser’s inequality gives rise to an analog of (10) as
sup Q(ξ ∈ A ) ≤ P0 (ξ ∈ A ) + 0.51/2 ( sup KL(P0 , Q))1/2 .
Q∈UdataQ∈Udatawhere KL(P0 , Q) denotes the KL divergence between P0 and Q, which, using our machinery describedabove, is not as tight as using (10) when the ultimate value of δ is less than 0.5.
Lastly, note that we have constructed Udata that satisfies (5) asymptotically. To provide finite-sample
bound, we can use concentration inequalities regarding the MLE estimators (see Korostelev and Korosteleva2011). Though we do not go into the details here, we note that our previous discussion is exact for Gaussian
distributions with unknown means.

Lam and Li
3

GENERALIZATIONS AND MIXTURE BASELINE DISTRIBUTIONSNote that an important element to determine the sample size requirement of our scheme is the quantity
sup χ 2 (P0 , Q) in (11). We have chosen a convenient choice in the last section for this quantity, which
Q∈Udata

involves choosing P0 and Udata . Within the class of φ - or f -divergences, the choice Udata is asymptoticallythe same. To express this more precisely, we re-parametrize the uncertainty set U˜data to be over the spaceof the parameter θ ∈ Θ, and work with sup χ 2 (P0 , Pθ ) instead. We letθ ∈U˜data(
U˜data ,Tθ ∈ Θ : (θ − θ̂n ) I(θ̂n )(θ − θ̂n ) ≤2
χ1−α,p

)

n

,

(14)

where I(θ̂n ) is the estimated Fisher information, and θ̂n is the MLE in which we highlight the dependence
on n. All divergence-based uncertainty sets Udata has an equivalent asymptotic form as (14).
The natural question to ask is which baseline measure P0 we should choose to generate the Monte Carlosample. For convenience, we denote Ddata (P0 ) = sup χ 2 (P0 , Pθ ). Ideally, we would like to choose P0 soθ ∈U˜datathat Ddata (P0 ) is computable and also minimized, which then leads to a computable and minimized required
number of sample N. To proceed, suppose further that Pθ has a density p(y; θ ) on Y . Furthermore, we
focus on P0 that has a density p0 (y) in the mixture formZ

p0 (y) =p(y; θ )µ(dθ ),U˜datawhere µ is a distribution on θ ∈ U˜data , and for convenience we call this associated class of distributionsP(U˜data ). This choice of density is handy to simulate from because we can simply sample θ ∼ µ(dθ )and then ξ ∼ Pθ .
To minimize Ddata (P0 ), we write
Ddata (P0 ) = supZ p(y; θ )p0 (y)

θ ∈U˜data Y

2− 1 p0 (y)dy(p(y; θ ))2dy − 1p0 (y)
θ ∈U˜data Y
Z

= sup

Z

= sup

θ ∈U˜data Y

We defineZ

L(µ, θ ) ,Y(p(y; θ ))2dy − 1.
U˜data p(y; θ )µ(dθ )R(p(y; θ ))2
dy andU˜data p(y; θ )µ(dθ )

R

(15)

φ (µ) , sup L(µ, θ ).
θ ∈U˜data

It follows from (15) that minimizing Ddata (P0 ) is equivalent to solving:min

max L(µ, θ ) =µ∈P(U˜data ) θ ∈U˜dataminµ∈P(U˜data )φ (µ).

(16)

Note that the right hand side of (16) is a convex minimization problem, thanks to the fact that 1/x isa convex function for x > 0.
However, even though the minimization in (16) is convex, solving it exactly is difficult due to theinner maximization. Thus we will aim for some heuristic methods to improve the choice of the baseline
distribution. To gain some insights, we first consider the value of φ (µ) in the case where we pick

Lam and Liµ = δθ̂n ∈ P(U˜data ), the point mass at the MLE estimator at δθ̂n , so that P0 = Pθ̂n . Now we want to see if
we can find some other µ ∈ P(U˜data ) with a smaller value of φ (µ). In this way, even if we cannot exactly
solve (16), we can still decrease the value of Ddata (P0 ). Here, we propose a mixture µ prop (dθ ) ∈ P(U˜data )on the boundary of U˜data that is easy to simulate and shows promising practical performance. Specifically,
we define the support of measure µ prop (dθ ) to beΘ(θ̂n ) = {θ ∈ U˜data : (θ − θ̂n )T I(θ̂n )(θ − θ̂n ) =2
χ1−α,pn}.

pNow, to sample θ ∼ µ prop (dθq), we first sample a random vector η ∈ R uniformly on the surface of2
χ1−α,pbe achieved by sampling from p number ofn . In particular, this canq2χ1−α,punit of length (see Muller (1959)).
independent standard normal variables and scale them to haven
−1/2Then, we set θ = θ̂n + (I(θ̂n ))η. In other words, we have chosen µ prop such that if θ ∼ µ prop (dθq), 2thenχ1−α,p(I(θ̂n ))1/2 (θ − θ̂n ) is uniformly distributed on the surface of the p-dimensional ball with radiusn .

the p dimension ball with radius

We want to simplify problem (16) by finding 0 < t < 1 that minimizes
φ (t) = φ ((1 − t)δθ̂n + tµ prop ),using line search or the bisection method. In the case where we have checked φ (1) < φ (0), we are certain
to make improvement over using µ as the point mass at θ̂Rn .
In general, if we cannot directly compute p0 (y) = U˜data p(y; θ )µ prop (dθ ) and solve for φ (µ prop ),S
an
R alternate is to use Monte Carlo samples of {θi }i≤S and approximate pˆ0 S (y) = ∑i=1 p(y, θi )/N =U˜data p(y; θ )µemprical (dθ ), where µempirical is the empirical distribution of these θi ’s. Then we try to
calculate φ ((1 − t)δθ̂n + tµempricial ) for different 0 < t < 1 to see if there is improvement.
4

NUMERICAL EXPERIMENTS

We present some numerical examples to illustrate the performance of our method. We first focus on thecomputational costs for different level of accuracy ε and dimension d. Next, we compare with the result
from standard SG. Finally, we demonstrate how choosing an effective baseline distribution P0 can reduce
the required sample size and decrease the computational cost.
We will perform experiments on multivariate Gaussian random variables and exponential randomvariables on single linear chance constrained programs (CCP) and joint linear CCP.
•••

For every experiment in each problem, we obtain an optimal solution x̂ by solving SG, and evaluate
the violation probability V (x̂, P) under the true probability measure Pθ0 (where θ0 denotes the true
parameter value) either through exact calculation by CDF or numerical expeiments of 10000 times.
Then we calculate the value ε̂ as the average of violation probability V (x̂, P) across 1000 cases. In
addition, we compute the 95th quantile of the violation probability and also the avarage of objective
value as “Ave.Obj.Val” in these 1000 cases.
In all examples we consider α = 0.05 and β = 0.05 for convenience and alternate between different
values of ε and d.
For each ε and d, we let Next to be the sample size we need if we can sample from P directly, andNamb be the Monte Carlo size using our method given the data size n. This number varies under
different choices of baseline measure P0 . Note that we must have n < Next < Namb .

Lam and Li
4.1 Single Linear Chance Constraints
We consider a single linear CCP
cT x,

min

x∈X ⊆Rd

subject to P((a + ξ )T x ≤ b) ≥ 1 − ε,where x ∈ Rd is the decision variable, a, c ∈ Rd and b ∈ R are constants. The random vector ξ ∈ Rd follows
N (ϑ , Id ) for some unknown ϑ . We suitably choose the parameters so that our problem is always feasible.
The true underlying measure has ϑ = 0. We choose a baseline measure P0 using ϑ̂n , the MLE estimator.
The χ 2 -distance in Udata can be explicitly computed in this case. For different levels of ε and d, we
compute Next and set n < N. Then we compare the values of Next and Namb to gain insight on the differenceof the computational cost in our method with that of the standard SG. Finally, we report “Ave.Obj.Val”
and ε̂ and its 95% quantiles Q95 from 1000 experiments. Table 1 shows our results.
Table 1: A Single Linear CCP for Gaussian with unknown mean.
Different levels of ε and d

n
Next
Namb
Ave.Ob j.Valε̂Q95

ε = 0.1
d=5

ε = 0.1
d = 10ε = 0.1
d = 20ε = 0.05
d=5

ε = 0.05
d = 10ε = 0.05
d = 2060
89342
0.80860.00940.0203100
154
585
0.86860.01910.0341

180
275
1010
0.90540.02210.0333100
181
748
0.80370.00710.0137

200
311
1144
0.85660.00910.0149

300
554
2209
0.89440.00970.0145Note that the data size is not enough to support a standard SG in all our considered settings. For
comparison, we show the Monte Carlo Namb needed in our method, which are quite big compared to Next ,though one should keep in mind that the Monte Carlo samples are easy to generate. As we can see, the ε̂from our method are all below the tolerance level, exhibiting the validity of our method.
Next, we explore the case where ξ ∼ exp(λ ) with some unknown λ . We set the true underlying measure
P to have λ = 1 and we choose the baseline measure P0 to use λ̂n , the MLE estimator. We consider thecase where ε = 0.01, d = 1 and summarize the results in Table 2. Again, the ε̂ we obtain is well below
the tolerance level.
Table 2: A Single Linear CCP for Exponential with unknown mean.

Value

n

NexactNamb

Ave.Ob j.Val

ε̂

Q95

100

299

1761

0.36440.000360.0014

Finally, we consider the case where ξ ∈ N (ϑ , Σ) with both parameters unknown. In such cases, we use
the Wilks’ theorem and construct a joint confidence region for both ϑ and Σ based on the χ 2 -distribution(Greene 2008; Arnold and Shavelle 1998). Specifically, We show the results for d = 1, ε = 0.01 fordemonstration. We set the true underlying measure P to have ϑ = 0, σ 2 = 1 and we choose the baseline
measure P0 to use (ϑ̂n , σ̂n ), the MLE estimator. We summarize the results in Table 3, which show similar
patterns as the previous cases.

Lam and Li
Table 3: A Single Linear CCP for Gaussian with unknown mean and unknown variance.

Value

n

NexactNamb

Ave.Ob j.Val

ε̂

Q95

200

299

1426

0.52070.00110.00444.2 Joint Linear Chance Constraints
We consider a joint linear CCP
minx∈X ⊆RdcT x,subject to P((A + ξ )T x ≤ b) ≥ 1 − ε,where x ∈ Rd is the decision variable and c ∈ Rd ,b ∈ Rl and A ∈ Rd×l are constants. The random vectorξ ∈ Rd×l satisfies vec(ξ ) ∼ N (0, Σ) where Σ ∈ Rdl×dl is some non-identity positive definite matrix. Again,
we suitably choose the parameters so that our problem is always feasible. Also, we choose different valuesof d, l and ε to demonstrate the results. We summarize the results in Table 4.
Table 4: A Joint Linear CCP.
Different levels of d, l and ε

n
Next
Namb
Ave.Ob j.Valε̂Q95

ε = 0.1
d=5
l = 10ε = 0.1
d = 10d = 15ε = 0.05
d=5
l = 10ε = 0.05
d = 10l = 1560
89342
0.63880.00170.0041100
154
585
0.65630.00240.0052100
181
748
0.63870.00170.0037

200
311
1144
0.66690.00440.0091In this joint CCP case, the number of samples Namb can grow quickly when ε decreases. Our ε̂ arestill well below the tolerance level, thus showing the validity of our approach. However, to reduce thecomputational costs, we should consider choosing a more efficient P0 . We will address this question in
the following subsection.
4.3 Different Choices of the Baseline Distribution
We investigate different choices of P0 , motivated from results in the last subsection that using P0 at theMLE estimator θ̂n may lead to a large Namb . For demonstration, we consider the case where the parametricfamily is Gaussian with unknown mean
R and unit variance. We propose different choices of the mixture
distribution µ(dθ ), compute p0 (y) = U˜data p(y; θ )µdθ and solve directly for Ddata (P0 ) to compute δ andNamb . Note that, in the case ξ ∼ N (ϑ , 1) and α = 0.05, the uncertainty set U˜data is simply the interval√ , ϑ̂n + 1.96
√ ].
[ϑ̂n − 1.96
nn
We compare 4 types of mixture on U˜data . We denote µ1 to be the point mass δθ̂n , µ2 to be the mixtureof two equal point mass at the two points that split the interval U˜data into three equal-distance pieces.
Then, we set µ3 to be the uniform distribution on the entire interval U˜data . Finally, µ prop is the uniform
distribution on the boundary of U˜data , which is just two equal point masses at two end points. We set

Lam and Li
ε = 0.05. We summarize the results in Table 5. As, we can see, in this example, µ prop has the best
performance as it gives the lowest number of required Monte Carlo samples.
Table 5: A Comparison among four choices of P0 for single CCP.

δθ̂n
µ2
µ3
µ propDdata (P0 )

δ

Namb

0.46837
0.426110.367020.28765

0.01530.01640.01820.0215

195
182
164
138

We also demonstrate the case of multivariate Gaussian with unknown mean. In this case, we can still
explicitly compute the p0 (y) under µ prop using the hypergeometric function defined in Nath (1951). We
show the comparisons in Table 6. We again set ε = 0.05, d = 2 and n = 80. In this case Next = 93. Our
choice of µ prop gives a much smaller Namb than using µ set as the point mass at θ̂n . Moreover, the smallvalue of ε̂ shows that our solution statistically satisfies feasibility for the CCP, which demonstrates that
this choice of µ prop is valid for our method.
Table 6: A Comparison among two choices of P0 for joint CCP.

δθ̂nµ propDdata (P0 )

δ

Namb

Ave.Obj.Val

ε̂

Q95

0.01550.005430.01530.0360305
1300.9126
0.92110.00160.01530.00800.0402ACKNOWLEDGMENTS
We gratefully acknowledge support from the National Science Foundation under grants CMMI-1542020
and CAREER CMMI-1653339/1834710.

REFERENCES
Arnold, B. C., and R. M. Shavelle. 1998. “Joint Confidence Sets for the Mean and Variance of a Normal
Distribution”. The American Statistician 52(2):133–140.
Ben-Tal, A., and A. Nemirovski. 2000. “Robust solutions of Linear Programming problems contaminatedwith uncertain data”. Mathematical Programming 88(3):411–424.
Blanchet, J., and Y. Kang. 2016. “Sample Out-Of-Sample Inference based on Wasserstein Distance”. arXiv
preprint arXiv:1605.01340.
Calafiore, G. C. 2017. “ Repetitive Scenario Design”. IEEE Transactions on Automatic Control 62(3):1125–
1137.
Calafiore, G. C., F. Dabbene, and R. Tempo. 2011. “ Research on Probabilistic Methods for Control System
Design”. Automatica 47(7):1279–1293.
Campi, M. C., and A. Carè. 2013. “Random Convex Programs with L1 -Regularization: Sparsity andGeneralization”. SIAM Journal on Control and Optimization 51(5):3532–3557.
Campi, M. C., and S. Garatti. 2008. “The Exact Feasibility of Randomized Solutions of Uncertain Convex
Programs”. SIAM Journal on Optimization 19(3):1211–1230.
Campi, M. C., and S. Garatti. 2018. “Wait-and-Judge Scenario Optimization”. Mathematical Programming 167(1):155–189.

Lam and Li
Carè, A., S. Garatti, and M. C. Campi. 2014. “FAST–Fast Algorithm for the Scenario Technique”. Operations
Research 62(3):662–671.
Chamanbaz, M., F. Dabbene, R. Tempo, V. Venkataramanan, and Q.-G. Wang. 2016. “Sequential Randomized
Algorithms for Convex Optimization in the Presence of Uncertainty”. IEEE Transactions on AutomaticControl 61(9):2565–2571.
De Farias, D. P., and B. Van Roy. 2004. “On Constraint Sampling in the Linear Programming Approach
to Approximate Dynamic Programming”. Mathematics of Operations Research 29(3):462–478.
Delage, E., and Y. Ye. 2010. “Distributionally Robust Optimization Under Moment Uncertainty with
Application to Data-Driven Problems”. Operations Research 58(3):595–612.
Erdoğan, E., and G. Iyengar. 2006. “Ambiguous Chance Constrained Problems and Robust Optimization”.
Mathematical Programming 107(1-2):37–61.
Esfahani, P. M., and D. Kuhn. 2015. “Data-Driven Distributionally Robust Optimization using the Wasserstein
Metric: Performance Guarantees and Tractable Reformulations”. Mathematical Programming:1–52.
Gao, R., and A. J. Kleywegt. 2016. “Distributionally Robust Stochastic Optimization with Wasserstein
Distance”. arXiv preprint arXiv:1604.02199.
Ghosh, S., and H. Lam. 2018. “Robust Analysis in Stochastic Simulation: Computation and Performance
Guarantees”. To appear in Operations Research, available as arxiv preprint at arXiv:1507.05609.
Glasserman, P., and X. Xu. 2014. “Robust Risk Measurement and Model Risk”. Quantitative Finance 14(1):29–58.
Greene, W. 2008. Econometric Analysis. 6th ed. Upper Saddle River, New Jersey: Pearson/Prentice Hall.
Hanasusanto, G. A., V. Roitch, D. Kuhn, and W. Wiesemann. 2015. “A Distributionally Robust Perspective on Uncertainty Quantification and Chance Constrained Programming”. Mathematical Programming 151(1):35–62.
Hanasusanto, G. A., V. Roitch, D. Kuhn, and W. Wiesemann. 2017. “Ambiguous Joint Chance Constraints
under Mean and Dispersion Information”. Operations Research 65(3):751–767.
Hong, L. J., J. Luo, and B. L. Nelson. 2015. “Chance Constrained Selection of the Best”. INFORMS
Journal on Computing 27(2):317–334.
Hu, Z., J. Cao, and L. J. Hong. 2012. “Robust Simulation of Global Warming Policies using the DICE
Model”. Management Science 58(12):2190–2206.
Hu, Z., and L. J. Hong. 2015. “Robust Simulation of Stochastic Systems with Input Uncertainties Modeledby Statistical Divergences”. In Proceedings of the 2015 Winter Simulation Conference, edited by
L. Yilmaz et al., 643–654. Piscataway, New Jersey: IEEE.
Jiang, R., and Y. Guan. 2016. “Data-Driven Chance Constrained Stochastic Program”. Mathematical
Programming 158(1-2):291–327.
Korostelev, A. P., and O. Korosteleva. 2011. Mathematical Statistics: Asymptotic Minimax Theory. Providence, Rhode Island: American Mathematical Society.
Lam, H. 2016. “Robust Sensitivity Analysis for Stochastic Systems”. Mathematics of Operations Research 41(4):1248–1275.
Lam, H. 2018. “Sensitivity to Serial Dependency of Input Processes: A Robust Approach”. Management
Science 64(3):1311–1327.
Lejeune, M. A., and A. Ruszczynski. 2007. “An Efficient Trajectory Method for Probabilistic ProductionInventory-Distribution Problems”. Operations Research 55(2):378–394.
Li, B., R. Jiang, and J. L. Mathieu. “Ambiguous Risk Constraints with Moment and Unimodality Information”.
Mathematical Programming:1–42.
Lim, A. E., J. G. Shanthikumar, and Z. M. Shen. 2006. “Model Uncertainty, Robust Optimization, andLearning”. In Tutorials in Operations Research: Models, Methods, and Applications for Innovative
Decision Making, edited by M. P. Johnson et al., 66–94. Catonsville, Maryland: INFORMS.

Lam and Li
Love, D., and G. Bayraksan. 2015. “Phi-Divergence Constrained Ambiguous Stochastic Programs for DataDriven Optimization”. Technical report, Department of Integrated Systems Engineering, The Ohio State
University, Columbus, Ohio.
Luedtke, J., and S. Ahmed. 2008. “A Sample Approximation Approach for Optimization with Probabilistic
Constraints”. SIAM Journal on Optimization 19(2):674–699.
Muller, M. E. 1959. “A Note on a Method for Generating Points Uniformly on N-Dimensional Spheres”.
Communications of the ACM 2(4):19–20.
Murr, M. R., and A. Prékopa. 2000. “Solution of a Product Substitution Problem Using Stochastic Programming”. In Probabilistic Constrained Optimization, edited by U. Stanislav, 252–271. Manhattan,
New York: Springer.
Nath, P. 1951. “Confluent Hypergeometric Function”. Sankhyā: The Indian Journal of Statistics:153–166.
Nielsen, F., and R. Nock. 2014. “On the Chi Square and Higher-Order Chi Distances for Approximatingf -Divergences”. IEEE Signal Processing Letters 21(1):10–13.
Pardo, L. 2005. Statistical Inference Based on Divergence Measures. New York: Chapman and Hall/CRC.
Petersen, I. R., M. R. James, and P. Dupuis. 2000. “Minimax Optimal Control of Stochastic Uncertain
Systems with Relative Entropy Constraints”. IEEE Transactions on Automatic Control 45(3):398–412.
Prékopa, A. 2003. “Probabilistic Programming”. In Handbooks in Operations Research and Management
Science, Volume 10: Stochastic Programming, edited by A. Ruszczynski and A. Shapiro. Amsterdam,
Netherlands: Elsevier.
Prékopa, A., T. Rapcsák, and I. Zsuffa. 1978. “Serially Linked Reservoir System Design Using Stochastic
Programing”. Water Resources Research 14(4):672–678.
Prékopa, A., and T. Szántai. 1978. “Flood Control Reservoir System Design Using Stochastic Programming”.
In Mathematical Programming in Use, edited by M. Balinski and C. Lemarechal, 138–151. Manhattan,
New York: Springer.
Schildbach, G., L. Fagiano, and M. Morari. 2013. “Randomized Solutions to Convex Programs with Multiple
Chance Constraints”. SIAM Journal on Optimization 23(4):2479–2501.
Shi, Y., J. Zhang, and K. B. Letaief. 2015. “Optimal Stochastic Coordinated Beamforming for Wireless
Cooperative Networks with CSI Uncertainty”. IEEE Transactions on Signal Processing 63(4):960–973.
Tseng, S.-H., E. Bitar, and A. Tang. 2016. “Random Convex Approximations of Ambiguous Chance
Constrained Programs”. In 2016 IEEE 55th Conference on Decision and Control, edited by A. Giua etal., 6210–6215. Piscataway, New Jersey: IEEE.
Vajda, I. 1972. “On the f -Divergence and Singularity of Probability Measures”. Periodica Mathematica
Hungarica 2(1-4):223–234.
Wiesemann, W., D. Kuhn, and M. Sim. 2014. “Distributionally Robust Convex Optimization”. Operations
Research 62(6):1358–1376.
Zhang, Y., R. Jiang, and S. Shen. 2016. “Ambiguous Chance-Constrained Bin Packing under MeanCovariance Information”. arXiv preprint arXiv:1610.00035.
Zymler, S., D. Kuhn, and B. Rustem. 2013. “Distributionally Robust Joint Chance Constraints with SecondOrder Moment Information”. Mathematical Programming:1–32.
AUTHOR BIOGRAPHIESHENRY LAM is an Associate Professor in the Department of IEOR at Columbia University. He received his
Ph.D. degree in statistics from Harvard University in 2011, and was on the faculty of Boston University andUniversity of Michigan before joining Columbia in 2017. His research focuses on Monte Carlo simulation,risk and uncertainty quantification, and stochastic optimization. His email address is khl2114@columbia.edu.
FENGPEI LI is a Ph.D. candidate in the Department of IEOR at Columbia University. He receives his B.S.
degree in mathematics from UCSD. His research focuses on Monte Carlo simulation, stochastic modeling,
stochastic optimization and robust optimization. His email address is fl2412@columbia.edu.