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

SIMULATING TAIL EVENTS WITH UNSPECIFIED TAIL MODELSHenry Lam

Clementine Mottet

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

Department of Mathematics and Statistics
Boston University
111 Cummington Mall
Boston, MA 02215, USA

ABSTRACT
Reliable simulation estimation builds on accurately specified input models. In the context of simulating
tail events, knowledge on the tail of the input model is especially important, yet is often hard to obtain
due to a lack of data. In this paper, we consider tail event estimation without any knowledge on the input
tail, but rather only making a general assumption that it is convex. We focus on the standard problemof estimating the probability for i.i.d. sum, and set out goal as to compute its worst-case bound among
all summand distributions that have convex tails. Our main procedure relies on a stochastic, and in asense infinite-dimensional, version of the Frank-Wolfe method in nonlinear programming. We demonstrate
through a numerical example how the level of knowledge on the tail of the summands relates to theconservativeness in computing bounds for the aggregate tail quantity.
1

INTRODUCTION

This paper explores simulation-based estimation of probabilities related to tail events without full knowledgeof the input model, especially its tail. The motivation is that, in most practical scenarios, information on
the tail of an input model is difficult to obtain because of the inavailability of data. Specifically, we shallconsider the standard problem of simulating the tail probability of i.i.d. sum, but each summand randomvariable has a distribution that is known only up to a certain threshold, but nothing beyond that. Our
goal is to construct a methodology that can output reasonable estimate for the tail probability under such
uncertainty.
The problem we post is generally ill-posed when no additional assumption is imposed, in the sensethat the tail probability estimate of the i.i.d. sum can be anywhere from 0 to 1. Our first step to avoid
such issue is to make a geometric assumption: the tail density of the i.i.d. random variables is convex.
This assumption is clearly satisfied by any known parametric distributions, including normal, exponential,
lognormal, Weibull etc., as well as distributions that are more pertinent to tail events like the generalized
extreme value distribution, and therefore we believe it is a reasonable assumption to make.
Our approach is a worst-case analysis: imposing convexity as a constraint, we aim to find the mostextreme tail probability of i.i.d. sum. This requires solving an optimization program that has decisionvariable being a probability density. This optimization has a nonlinear objective function and is infinitedimensional. Our main contribution of the paper is to demonstrate a stochastic version of the Frank-Wolfe(FW) method that is well-suited for the problem, including being able to run in the infinite-dimensional
space by exploiting the problem structure. In fact, due to this structure, our algorithm will operate on a
space of what we call “augmented” probability densities, which is designed to handle sequences of densities
that are not tight and have masses escaping to infinity. In overall, the implementability of the algorithm
will provide a tractable way to obtain robust estimate of tail probabilities when the input model’s tail is
not fully known.

Lam and Mottet
We close this introduction by discussing some related work. Our worst-case framework is inspired by
the literature on distributionally robust optimization (Delage and Ye 2010, Ben-Tal et al. 2013), which
considers decision-making in stochastic environment where the underlying probability is not fully known
and is assumed to lie in an uncertainty set. There is also literature on so-called moment problems andChebyshev’s inequalities, applied to stochastic programming and decision analysis (Birge and Wets 1987,
Smith 1995). In contrast to these literature, we impose tail convexity constraint, studied in Lam and Mottet(2015) and is also related to the general framework of moment problems in distributional class in Popescu(2005). References on the FW method can be found in Bertsekas (1999), and our stochastic version of FW
is related to Ghosh and Lam (2015) and Goeva et al. (2014), which consider optimization over discretized
input models as means for measuring model uncertainty or carrying out model calibration.
The paper is organized as follows. We present our formulation and notation in Section 2, and provide
an overview of our procedure in Section 3. Then we discuss several aspects of the procedure, including
the gradient form in Section 4, the solution to linearized subprograms in Section 5, and the constructionof augmented probability densities in Section 6. We lay out the algorithmic details of our FW procedurein Section 7, followed by some numerics in Section 8. We discuss future work in Section 9.
2

FORMULATION

We are interested in estimating P(Sn > b) where Sn = X1 + · · · + Xn , Xi ’s are i.i.d. random variables and b issome large number. We assume that Xi has density f (x) for x ∈ R but is known only up to some thresholda. In other words, f (x) = g(x) for some known function g for x ≤ a.
Beyond a, we assume that f is convex. We denote
A = { f (x), x ∈ R : f (x) = g(x) for x ≤ a, f (x) is convex for x ≥ a}(1)

as the set of densities that are equal to g for x ≤ a and is convex beyond a. Our goal is to compute theworst-case tail probability
max Pf (Sn > b)
(2)
f ∈A

where Pf denotes the probability measure generated by the i.i.d. random variables each with density f .
It is easy to check that the feasible region A is convex. However, the objective function in (2), which
can be written asZ
Z
Pf (Sn > b) = · · ·f (x1 ) · · · f (xn )dx1 · · · dxn(3)x1 +···+xn >b

is not linear, nor convex with respect to f in general. In view of this, we shall set our target as to obtain alocal optimum for (2).
3MAIN ELEMENTS OF OUR PROCEDURE

The strategy that we shall undertake borrows from the FW method (Bertsekas 1999; also known as theconditional gradient method) in nonlinear programming. This method applies when the feasible region
is convex. It is an iterative scheme that, at each iteration, linearizes the objective function and solves a
subprogram with linear objective and the given feasible region in order to find a direction to move along.
To explain further, let us denote Z( f ) = Pf (Sn > b). At each iteration k, given the current solution fk , one
would need to solve
maxh∇Z( fk ), f − fk i
(4)
f ∈A

where ∇Z( f ) denotes some notion of the gradient of Z(·) with respect to f , and h·, ·i is some suitable innerproduct. Say the optimal solution of (4) is rk . Then one would update the solution fk+1 = fk + εk (rk − fk ) =(1 − εk ) fk + εk rk for some step size εk .

Lam and Mottet
For our particular problem, it is difficult to evaluate the objective function or its gradient directly
because this would involve high-dimensional convolution when n is large. So one has to resort to simulation\
estimation. Hence one would need to replace ∇Z( fk ) by some estimate ∇Z(fk ), and the corresponding
stepwise subprogram will be\maxh∇Z(fk ), f − fk i
(5)
f ∈A

with the update being fk+1 = (1 − εk ) fk + εk r̂k where r̂k is the optimal solution of (5)In the following sections, we will present the above notions and steps in detail.
4

LINEARIZATION OF OBJECTIVE FUNCTION

We first discuss the notion of gradient for Pf (Sn > b) with respect to the density f . This is summarized as:
Lemma 1 For any probability densities f1 and f2 , we haveZ ∞

Pf2 (Sn > b) = Pf1 (Sn > b) +−∞ζ (x; f1 )( f2 (x) − f1 (x))dx + O(k f1 − f2 k21 )(6)

as k f1 − f2 k1 → 0, whereζ (x; f1 ) = nPf1 (Sn−1 > b − x)R∞and k f1 − f2 k1 = −∞ | f1 (x) − f2 (x)|dx denotes the L1 -distance between f1 and f2 .
The proof follows from an expansion

Proof of Lemma 1.

Pf2 (Sn > b)
Z

=
Z

=
Z

=

···
···
···

n

Z

∏ f2 (xi )dxi

∑ni=1 xi >b i=1
Zn∏( f1 (xi ) + ( f2 (xi ) − f1 (xi )))dxi∑ni=1 xi >b i=1nnZ∑ni=1 xi >b∏ f1 (xi ) + ∑i=1∏ f1 (x j )( f2 (xi ) − f1 (xi ))i=1 j=1,...,n
j6=i

n

!

r

+∑

∑

r=2 (i1 ,...,ir )∈P(n,r)∏

j=1,...,nj6=iu ,u=1,...,rf1 (x j ) ∏( f2 (xil ) − f1 (xil ))
l=1

n

∏ dxi
i=1

where P(n, r) denotes the set of all permutations of r items from {1, . . . , n}. Now, integrating the firstterm in the integrand gives Pf1 (Sn > b). Integrating the second term gives


n Z

∑

i=1

Z

= n

Z




···

Z

∏ f1 (x j )dx j  ( f2 (xi ) − f1 (xi ))dxi∑l=1,...,n xl >b−xi j=1,...,n
j6=i
l6=i

Pf1 (Sn−1 > b − x)( f2 (x) − f1 (x))dx

by noting that xi ’s are exchangeable. Note that each summand in the third term contains two or morefactors of f2 (xil ) − f1 (xil ), and the integral of the term is absolutely bounded by
Z Z

C| f2 (x) − f1 (x)|| f2 (y) − f1 (y)|dxdy

for some C > 0. Combining the above gives (6).

Lam and Mottet
The function ζ (·; f1 ) can be interpreted as the first order gradient of Pf (Sn > b) at f1 . When f1 andf2 are very close in terms of the L1 -distance, the second order term in (6) becomes relatively negligible.
The function ζ (·; f1 ) is essentially the Gateaux derivative of Z( f1 ) and is related to the so-called influencefunction (Huber 2011) in robust statistics, which is an important tool in measuring the infinitesimal effecton a given statistics due to changes in the distribution of data.
Therefore, from (6), given a density f1 , the first order approximation of Pf2 (Sn > b) centered at f1 isgiven by
Pf2 (Sn > b) ≈ Pf1 (Sn > b) + E f2 [ζ (X; f1 )] − E f1 [ζ (X; f1 )].
This suggests that at each iteration k in our optimization procedure, given fk , one would solve max f ∈A E f [ζ (X; fk )]as the subprogram. When ζ (X; fk ) is not computable exactly, an estimated version ζ\(X; fk ) will be substituted.
5

WORST-CASE CHARACTERIZATION FOR LINEAR OBJECTIVES

We discuss how to solve max f ∈A E f [ζ (X)] for a given function ζ (·). This borrows results from Lam andMottet (2015). To formulate the optimization more explicitly, let us pin down the necessary informationfrom the known function g. We define:
1. η = g(a) as the value of g at a.
2. −ν =R g0− (a) as the left derivative of g at a.
Rag(x)dx.
3. β = a∞ g(x)dx as the tail distribution at a, which is equal to 1 − −∞We assume that η, ν, β are all positive. Then the formulation max f ∈A E f [ζ (X)] can be written asamax fζ (x)g(x)dx + a∞ ζ (x) f (x)dxR−∞∞subject to a f (x)dx = β
f (a) = η
f+0 (a) ≥ −ν
f convex for x ≥ af (x) ≥ 0 for x ≥ a.

R

R

(7)

In this optimization, the first constraint makes sure that f is a valid probability density by equating thetail probability at a to β . The second constraint states that f has to be a continuous extrapolation from g,since f is convex. The third constraint specifies that the right derivative of f at a needs to be at least −ν,
because of the convexity assumption again. The whole set of constraints therefore ensures that f is a valid
density that is a convex extrapolation from g.
The following result is adopted from Lam and Mottet (2015):
Theorem 1 Consider the optimization max f ∈A E f [ζ (X)]. Suppose that ζ (·) is bounded, non-negativeand non-decreasing, and η, ν, β > 0. If η 2 > 2β ν, there is no feasible solution. Otherwise, optimality is
characterized by either one of the two scenarios:
1. There is an optimal density that is piecewise linear with at most two line segments for x > a. The
first line segment has slope −ν and the second line segment hits the x-axis. (The first line segment
can be degenerate, i.e. zero length.)
2. There is no optimal density. Instead, there is a sequence of densities f (s) , s = 1, 2, . . . whose objective
value converges to the optimum. Each of these f (s) has at most three line segments. The first line
segment has slope −ν and the third line segment hits the x-axis. The third line segment becomescloser and more parallel to the x-axis as s → ∞. (The first and the second line segments can bedegenerate.)Lam and Mottet
Hence, solving (7) reduces to finding the kinks of piecewise linear functions. Algorithm 1 depicted
below can classify between the above two scenarios and also solve for the optimal density or sequence ofdensities explicitly.
Let us explain the output of Algorithm 1. The output f (x) encodes either an optimal density, in thefirst case of Theorem 1, or the pointwise limit of a sequence of densities f (s) that converges to optimality,
in the second case of Theorem 1. We emphasize the fact that in the second case, the algorithm does not
return the sequence f (s) but only its pointwise limit f (x). Note that the limit density f (x) is not a valid
density itself, in the sense that it does not integrate to 1. The quantity q∗ encodes the amount of probability
mass that has “escaped” to infinity.
The variable x1 represents the location of the first kink in the optimal piecewise linear density (or limit
density). The function W (x1 ) is introduced to find this first kink and, in the second case of Theorem 1,
an additional function V (x1 , ρ) is needed. When x̃1 = argmaxx1 ∈[0,µ−ε]W (x1 ) 6= µ − ε, we have q∗ = 0,
which signifies that the optimality characterization falls into the first case in Theorem 1. On the other
hand, if x̃1 = µ − ε, then q∗ > 0 which brings to the second case. The tolerance parameter ε is introduced
to detect whether the maximum of W (x1 ) or V (x1 , ρ) occurs at x1 = µ. Note that the definitions of W (x1 )and V (x1 , ρ) both have singularities at x1 = µ (it can be checked that the values of W (x1 ) and V (x1 , ρ)remain bounded as x1 → µ). Introducing ε and checking x̃1 = µ − ε is an implementable way to detectthat the solution is at the singularity µ.
In the last scenario in the algorithm, when x1∗ = µ and ρ ∗ = µ 2 , we define p∗1 = 1 and x2∗ = x1∗ , or inother words output f (x) = η − ν(x − a) for a ≤ x ≤ x1∗ + a and 0 for x ≥ x1∗ + a.

Lam and Mottet
Algorithm 1 Procedure for solving (7)Input: η, ν, β , ζ (x), and a small “tolerance” parameter ε > 0
Definition: Let
Z xZ uH(x) =

ζ (v + a)dvdu
0

λ

0

H(x)
x2
x→∞


σ − µ2σ − µx1µ − x12H(x1 ) +H
σ − 2µx1 + x1 2
µ − x1σ − 2µx1 + x12



2 !
2

ρ
−
µx
ρ
−
µx
ρ − µ2
µ
−
x
1
1
1
H
H(x1 ) − λ x12 +−λ+λσµ − x1µ − x1
ρ − 2µx1 + x12
ρ − 2µx1 + x12= lim supW (x1 ) =V (x1 , ρ) =Initialization:µ ← ηνσ ← 2βνq∗ ← 0. By default, the escaping mass is null.

Procedure:if σ < µ 2 then
STOP. There is no feasible solution.
else if σ = µ 2 thenx̃1 ← µρ∗ ← µ2elsex̃1 ← argmaxx1 ∈[0,µ−ε]W (x1 )if x̃1 6= µ − ε thenρ∗ ← σelse(x̃1 , ρ ∗ ) ← argmaxx1 ∈[0,µ−ε],ρ∈[µ 2 ,σ ]V (x1 , ρ)
if x̃1 = µ − ε thenx̃1 ← µρ∗ ← µ2
end ifq∗ ← ν2 (σ − ρ ∗ )end ifExclusion of trivial scenarios.




















Case 1 : Light tail. 















Treatment of non-trivial scenarios.











Case 2 : Heavy tail.























x1∗ ← x̃1ρ ∗ −µ 2p∗1 ← ρ ∗ −2µx
∗ +x∗ 2x2∗ ←

ρ ∗ −µx1∗
µ−x1∗

1

1


 η − ν(x − a)
η − νx1∗ − ν(1 − p∗1 )(x − a − x1∗ )f (x) ←
0end ifreturn ( f (x), q∗ )










for a ≤ x ≤ x1∗ + a
for x1∗ + a ≤ x ≤ x2∗ + a
for x ≥ x2∗ + aLam and Mottet
6

AUGMENTED PROBABILITY DENSITIES

6.1 The Calculus of Augmented Probability Densities
The type of outputs in Algorithm 1 forms a natural space for running our FW algorithm, and here we shalldescribe it in moreprobability density as ( f (x), q), where f (x) ≥ 0 for
R ∞ detail. We define an “augmented”
R∞−∞ < x < ∞, −∞f (x)dx ≤ 1, and q = 1 − −∞f (x)dx.
We define the expectation of a bounded function h(x), with limit h(∞) = limx→∞ h(x), under anaugmented density ( f (x), q) asZ ∞E f ,q [h(X)] =h(x) f (x)dx + h(∞)q.

(8)

−∞

More complicated functionals of an augmented density can be defined by viewing q as an “escaping”
probability mass to positive infinity. For instance, for X1 , X2 i.i.d. with augmented density ( f (x), q), wewould define
Z ∞Z ∞

Pf ,q (X1 + X2 > b) =−∞ −∞I(x1 + x2 > b) f (x1 ) f (x2 )dx1 dx2 + (1 − (1 − q)2 ).

To interpret the above, one can think that if say a particular realization of X1 has “escaped”, then X1 mustbe larger than b. The first term above corresponds to the part that both X1 and X2 do not escape. The
second part corresponds to the scenario that at least one of them escape.
To simulate quantities dictated by an augmented density ( f (x), q), one can use a mixture of f˜(x) =f (x)/(1 − q) and the escaping mass q, by noting that f˜(x) is a valid density. For instance, for X1 , X2
i.i.d. with augmented density ( f (x), q),f (x1 ) f (x2 )dx1 dx2 + (1 − (1 − q)2 )1−q 1−q−∞ −∞= (1 − q)2 Pf˜(X1 + X2 > b) + (1 − (1 − q)2 )Pf ,q (X1 + X2 > b) = (1 − q)2

Z ∞Z ∞

I(x1 + x2 > b)(9)where Pf˜(·) is the probability measure generated by i.i.d. X1 and X2 each under the density f˜, and(1 − (1 − q)2 ) above can be evaluated exactly.
For another example, say Pf ,q (X1 + X2 ∈ [b, c]) for X1 , X2 i.i.d. with augmented density ( f (x), q), we
can express as
Pf ,q (X1 + X2 ∈ [b, c]) = (1 − q)2 Pf˜(x1 + x2 > b).
The term corresponding to the escaping mass disappears now because if say X1 escapes, then the event
X1 + X2 ∈ [b, c] cannot happen.
6.2 Consistency between Optimal Density Sequence and Augmented Density
One key observation when building the above rules for augmented densities is that they are consistent withthe limiting behavior of the optimal sequence of densities in the second case of Theorem 1. It can beshown that for an optimal density sequence for problem (7) given by f (s) that has pointwise limit f ,E f ,q [h(X)] = lim E f (s) [h(X)]s→∞where E f ,q [h(X)] is defined by rule (8). Moreover, consistency is also preserved for the tail probability ofi.i.d. sum, namely
Pf ,q (Sn > b) = lim Pf (s) (Sn > b) = (1 − q)n Pf˜(Sn > b) + (1 − (1 − q)n )s→∞where Pf ,q is defined by rule (9) and f˜(x) = f (x)/(1 − q). The proof of these results will be shownelsewhere.

Lam and Mottet
7

ITERATIVE SCHEME VIA FRANK-WOLFE STOCHASTIC APPROXIMATION

Our FW algorithm will operate in the augmented density space. Since the linearization at each iteration
has to be simulated, it can be viewed as a constrained stochastic approximation (SA) method. It is statedas:
Algorithm 2 Iterative procedure for solving (2)Input: the function g(x), x ≤ a, the level b, step size εk , sample size mk for each step k.
Initialization: an augmented density ( f1 (x), q1 ) where f1 (x) meets the constraints of problem (7). For
instance, one can take ( f1 (x), q1 ) to be the augmented density representing argmax f ∈A Pf (X > b), i.e. theoutput from Algorithm Algorithm 1 by putting ζ (x) = I(x > b) where I(·) denotes the indicator function.
(Other initializations will also work.)
Procedure: For each iteration k = 1, 2, . . ., given an augmented density ( fk (x), qk ):1. Compute1 mkζˆk (x) =∑ I(Sn−1 > b − x)
mk i=1
using mk sample paths of Sn−1 , where Sn−1 = X1 + · · · + Xn−1 for i.i.d. Xi ’s generated under f˜k (x) =fk (x)/(1 − qk ).
2. Run Algorithm 1 with input function ζˆk (x) to get an output (rk (x), uk ).
3. Update ( fk+1 (x), qk+1 ) = (1 − εk )( fk (x), qk ) + εk (rk (x), uk ).
We explain some details related to the subprogram step. By Lemma 1, one should post the subprogramat iteration k as max f ∈A E f [nPfk ,qk (Sn−1 > b − X)], where ( fk (x), qk ) is the current augmented density. Thenby the rules in Section 6.1 we have
Pfk ,qk (Sn−1 > b − x) = (1 − qk )n−1 Pf˜k (Sn−1 > b − x) + (1 − (1 − qk )n−1 ).
Hence optimizing E f [Pfk ,qk (Sn−1 > b − X)] is equivalent to optimizing E f [Pf˜k (Sn−1 > b − x)]. The quantityζˆk (x) is an unbiased estimator for P ˜ (Sn−1 > b − x).
fk

Note that the updating step ( fk+1 (x), qk+1 ) = (1 − εk )( fk (x), qk ) + εk (rk (x), uk ) represents a mixture of( fk (x), qk ) and (rk (x), uk ), which is defined as the mixture over both the density part and the escaping mass
part. In other words, we have fk+1 (x) = (1 − εk ) fk (x) + εk rk (x) and qk+1 = (1 − εk )qk + εk uk .
8

NUMERICAL EXAMPLE

Consider the setting where f is known to be an exponential distribution with rate λ , up to a. We apply
Algorithm 2 to estimate a local optimum for max f ∈A Pf (Sn > b). We set n = 8, b = 10, and λ = 1. Wevary a from the 70-th percentile of the exponential distribution to the 99-th percentile. For each a, weapply our Algorithm to find the locally optimal augmented density, and we simulate P(Sn > b) using thatto get an estimate of the maximum objective value.
Note that even though the gradient estimator ζˆk (x) at each iteration is unbiased, the estimated best
feasible direction obtained through solving the optimization subprogram is in general biased. Consequently,
in our implementation we use a growing sample size mk = 10k1.1 along the iterations. We use a standard
SA step size specification εk = 1/k. Some theoretical guarantees regarding this specification for FWSA,
though in a separate context, are reported in Ghosh and Lam (2015).
Since g is taken to be Exp(1), we know that the actual distribution of Sn is a Gamma distribution withshape and rate parameters equal to n and 1 respectively. The actual value of P(Sn > b) is therefore the taildistribution of the Gamma distribution evaluated at the point b. This provides a benchmark for measuring
the level of conservativeness of the output from Algorithm 2.

Lam and Mottet
Figure 1 shows the estimated value of Pf100 (Sn > b), i.e. we ran 100 iterations, for different valuesof the threshold a. In terms of convergence, we found that in this example Algorithm 2 almost reachedequilibrium in one iteration, and so 100 iterations was likely more than enough. On the other hand, we also
tried the same set of experiments for the corresponding minimization problem, and found that the algorithm
converged in around 20 iterations. In Figure 1, it can be seen that, when a is at the 70-th percentile, theestimate of P(Sn > b) of about 0.7 is very conservative, compared to the true value of around 0.2 dictatedby the Gamma distribution. Nevertheless, the estimate gets progressively less conservative as a gets larger,
until it drops to almost the true value at the 99-th percentile.

Figure 1: Estimated optimal upper bound of Pf100 (S8 > 10) for Xi ∼ Exp(1) known up to different percentiles.
Lastly, we provide some discussion on how to simulate f˜k at each iteration, which is used to generatecopies of Sn−1 and subsequently ζˆk (x) in Algorithm 2. To simulate a copy from f˜k , we can use the methodof inverse distribution function. First, since fk is a mixture of augmented densities, we start with generating
one of the components with the mixture probabilities. This component will have density λ e−λ x up to a,and then followed by a straight line with some slope −s1 up to x1 , and slope −s2 up to x2 which touches
the x-axis. To simulate from this component conditional on being not escaping, we can generate a uniform

Lam and Mottet
random variable U over the interval [0, 1 − qk ]. Then, letting F(a) = 1 − e−λ x and f (a) = λ e−λ a , we output
 1− λ log(1 −U)





for U ≤ F(a)


√ 2


f
(a)+
f
(a)
−2s
(U−F(a))

1
 a+
s1
2

for F(a) ≤ U ≤ F(a) + f (a)(x1 − a) − s1 (x12−a)


q


s (x −a)2
f (a)−s1 (x1 −a)+ ( f (a)−s1 (x1 −a))2 −2s2 (U−F(a)− f (a)(x1 −a)+ 1 12
)


x
+

1

s
2

2

for F(a) + f (a)(x1 − a) − s1 (x12−a) ≤ U ≤ x2 .

This can be seen by solving for the inverse at each segment (the exponential, and the two quadratic segments)
in the distribution function. For the first quadratic segment, one needs to solve, for a given U,
F(a) + f (a)(x − a) −s1 (x − a)2=U2

in x, and for the second quadratic segment, one solves
F(a) + f (a)(x1 − a) −9s1 (x1 − a)2s2 (x − x1 )2+ ( f (a) − s1 (x1 − a))(x − x1 ) −
= U.
2
2

DISCUSSION AND FUTURE WORKThis paper demonstrates an iterative scheme, based on a stochastic version of the FW method, to findworst-case estimates of the tail probability of i.i.d. sum assuming the distribution of the summands is known
only up to a specific threshold. We have described how to handle the infinite-dimensional nature of the
problem, and also the possibility of non-existence of an optimal solution in any stepwise subprogram byintroducing the notion of augmented probability density. Several future investigation directions are in line.
The first to obtain some theoretical convergence guarantees of the proposed algorithm. Second, for a highthreshold b, naive Monte Carlo for estimating ζˆk (x) in each iteration, like what we have used, will be too
slow. The use of importance sampling will be investigated. Third, we also plan to generalize the result
to other types of objectives in addition to the tail probability of i.i.d. sum. Lastly, all results in this paper
assume complete information on the non-tail part of the distribution, whereas in most cases it is partiallyknown through data. In the future, we plan to extend our analysis to joint uncertainty on both the tail andthe non-tail regions.
ACKNOWLEDGMENTS
We gratefully acknowledge support from the National Science Foundation under grants CMMI-1400391
and CMMI-1436247.
REFERENCES
Ben-Tal, A., D. Den Hertog, A. De Waegenaere, B. Melenberg, and G. Rennen. 2013. “Robust Solutions of
Optimization Problems Affected by Uncertain Probabilities”. Management Science 59 (2): 341–357.
Bertsekas, D. P. 1999. Nonlinear Programming. Athena scientific Belmont.
Birge, J. R., and R. J.-B. Wets. 1987. “Computing Bounds for Stochastic Programming Problems by Meansof a Generalized Moment Problem”. Mathematics of Operations Research 12 (1): 149–162.
Delage, E., and Y. Ye. 2010. “Distributionally Robust Optimization under Moment Uncertainty with
Application to Data-Driven Problems”. Operations research 58 (3): 595–612.
Ghosh, S., and H. Lam. 2015. “Computing Worst-Case Input Models in Stochastic Simulation”. submitted,
available at: http://www-personal.umich.edu/ khlam/files/robustSAOR3.pdf .

Lam and Mottet
Goeva, A., H. Lam, and B. Zhang. 2014. “Reconstructing Input Models via Simulation Optimization”. In
Proceedings of the 2014 Winter Simulation Conference, 698–709. IEEE Press.
Huber, P. J. 2011. Robust Statistics. Springer.
Lam, H., and C. Mottet. 2015. “Tail Analysis Without Tail Information: A Worst-Case Perspective”.
submitted, available at: http://www-personal.umich.edu/ khlam/files/tailOR7.pdf .
Popescu, I. 2005. “A Semidefinite Programming Approach to Optimal-Moment Bounds for Convex Classesof Distributions”. Mathematics of Operations Research 30 (3): 632–657.
Smith, J. E. 1995. “Generalized Chebychev Inequalities: Theory and Applications in Decision Analysis”.
Operations Research 43 (5): 807–825.
AUTHOR BIOGRAPHIESHENRY LAM is an Assistant Professor in the Department of Industrial and Operations Engineering at
the University of Michigan, Ann Arbor. He graduated from Harvard University with a Ph.D. degree instatistics in 2011, and has been an Assistant Professor in the Department of Mathematics and Statistics
at Boston University until 2014. His research focuses on stochastic simulation, rare-event analysis, andsimulation optimization, with application interests in service systems and risk management. His email
address is khlam@umich.edu.
CLEMENTINE MOTTET is a Ph.D. student in the Department of Mathematics and Statistics at Boston
University. She joined their PhD program in 2013 after graduating from the the National Institute of Applied
Sciences of Toulouse with a Master’s Degree in Applied Mathematics. She is currently conducting researchunder the supervision of Professor Lam. Their collaborative research has been focusing on robust estimationof tail quantities under limited information, which is relevant in fields related to risk management. She
can be reached at cmottet@bu.edu.