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

RARE-EVENT SIMULATION WITHOUT STRUCTURAL INFORMATION: A
LEARNING-BASED APPROACH
Zhiyuan Huang

Henry Lam

Department of IOE
University of Michigan
1205 Beal Avenue
Ann Arbor, MI 48105, USA

Department of IEOR
Columbia University
500 W. 120th Street
New York, NY 10027, USA

Ding Zhao
Department of Mechanical Engineering
Carnegie Mellon University
5000 Forbes Avenue
Pittsburgh, PA 15213, USA

ABSTRACT
Importance sampling has been extensively studied as a variance reduction tool in rare-event simulation.
The design and efficiency of this method often relies on structural knowledge about the target problem inhand. In this paper, we consider a simple rare-event setting, driven by Gaussian input variates, where there
is no direct information about the system of interest, or where the system of interest can only be evaluatedthrough a “black-box” model. This setting mimics, on a basic level, situations where the system of interest
could be too complex to analyze as in the conventional importance sampling literature. We investigate an
approach based on two-stage sampling, where the first stage learns the rare-event set of interest, and thesecond stage uses a change of measure that exploits the first-stage finding. We present some guarantees
and discuss comparisons of this scheme with more naive approaches.
1

INTRODUCTION

Rare-event simulation concerns the estimation of a probability that is very small, which usually serves
as a risk measurement of the occurrence of catastrophic events. Because of the small target probability,
crude Monte Carlo typically bears a huge relative error, and consequently requires a huge number ofsamples to obtain a reasonably accurate estimate. To overcome this issue, importance sampling (IS) is
commonly suggested as an effective variance reduction method in such contexts (early references include,e.g., Siegmund 1976; Glynn and Iglehart 1989; Sadowsky and Bucklew 1990 etc.).
To design an efficient IS estimator, existing approaches often involve analyzing the structures ordynamics of the system of interest. A large literature studies the use of large deviations to invent andmathematically justify IS procedures (see, e.g., the surveys Bucklew 2013; Asmussen and Glynn 2007;
Rubinstein and Kroese 2016; Glasserman 2013; Juneja and Shahabuddin 2006; Blanchet and Lam 2012).
This line of studies requires, in particular, thorough knowledge about the system or the underlying functionat hand. Another line of studies use techniques such as Markov-chain Monte Carlo (MCMC) to sample
from the rare-event set of interest, or approximately from the conditional distribution given the occurrenceof the rare event (Botev et al. 2013; Grace et al. 2014; Chan and Kroese 2012).

Huang, Lam, and Zhao
In this paper, we study a setting where we do not have direct information about the underlying
performance function. This happens when, e.g., the function is only accessible as a “black-box”, so that it
can be revealed only at given input values chosen by the user. This situation, on a basic level, mimics thesetting where the function may be too complex to mathematically analyze for designing IS schemes. It is
motivated by several recent applications of IS in self-driving cars evaluation (Zhao et al. 2017; Huang et al.
2017a), neural network training (Bengio and Senécal 2008), and reinforcement learning (Frank, Mannor,
and Precup 2008; Shelton 2001). This is in addition to more established applications in, e.g., finance(Glasserman et al. 1999), queueing (Dupuis et al. 2007; Dieker and Mandjes 2006; Blanchet et al. 2009),
structural reliability (Au and Beck 1999; Dubourg et al. 2013; Grabaskas et al. 2016; Heidelberger 1995)
and manufacturing (Kanj et al. 2006).
As a motivating example, suppose we have a black-box prediction model ŷ(x) (built from, e.g., amachine learning model) for the response of a physical system y(x) that is unknown and needs to beobserved. Here x is the input that is random. Think, in the setting of self-driving car, that x encodes thesurroundings and y(x) is the movement of the surrounding objects in the next time unit. The action of
the self-driving car relies on the prediction model, and there would be an accident or a serious loss if theprediction deviates from the physical outcome too much. Here, we can use g(x) = d (ŷ(x), y(x)) to capture
this deviation, where d(·, ·) denotes some distance measure. The probability P(g(x) ≥ γ) captures the riskof a serious loss. Since a typical machine learning output ŷ(x) can be highly unstructured and reminiscentof a black box, using IS to speed up the estimation of such type of probability falls into our study setup.
This paper focuses on a simple example of the above setting. In particular, we assume the input variates
are Gaussian and that the rare-event set has one “dominating point” (Sadowsky and Bucklew 1990; Dieker
and Mandjes 2006; Blanchet and Li 2011; Adler et al. 2012). These assumptions could be possibly relaxed,e.g., to Gaussian mixtures that can be used to model or approximate a wider class of distributions, andto rare-event sets that have multiple dominant points. We propose a two-stage sampling scheme in theabove situations. In the first stage, we sample from the space using a space-filling design. This stage
serves to learn information about the rare-event set, in particular the dominant point. In the second stage,
we construct IS via an exponential tilting based on the estimated dominant point. From the viewpoint of
sequential learning, one can think of the first stage as the “exploration” phase and the second stage as“exploitation”.
We present some theoretical error bounds on this approach. Our analysis, in particular, shows theimpact of using exponential tilting based on a point that only approximates the dominating point in termsof their distance, which is in general exponential. However, there is a concentration-type behavior in our
first-stage learning on the dominant point, which offsets this exponential dependence. One of the insights
is that, under the particular setting we consider, the two-stage approach is advantageous against more naive
approaches, such as pure uniform sampling, if the sample space is large compared to the location of therare-event set.
Our study here is perhaps closest to the setting of cross-entropy method (Rubinstein and Kroese 2004,
De Boer et al. 2005), which can be used when there is no information about the underlying structure of
the rare event. This approach relies on a sequential use of sample average approximation to minimize theKullback-Leibler divergence between an IS distribution and the zero-variance distribution to update theIS parameters. Part of our motivation is that the cross-entropy method, though demonstrated to be very
effective, typically requires an explicit identification of a threshold parameter controlling the rarity levelthat can be sequentially adjusted. In contrast, our method does not need this parameter, and instead we rely
on a concentration-type exploration of the rare-event set in the first stage that informs our IS. Furthermore,
our study is also related to a line of IS work based on learning the rare-event set, for instance, through firstor second order approximation (Zhao and Ono 1999), quadratic response surface (Das and Zheng 2000),
neural network (Papadrakakis and Lagaros 2002), support vector machine (Hurtado 2004; Bourinet et al.
2011; Huang et al. 2018), kriging or Gaussian process (Echard et al. 2011; OHagan 2006; Dubourg et al.
2013), and the utilization of monotonic features (Huang et al. 2017b). These methods are all heuristicallyHuang, Lam, and Zhao
built on the use of data to approximate the rare-event set, and there has been no theoretical guarantees
on their performances. Our study here could provide some basic mathematical insights on the errors andtheir impacts on Monte Carlo from such type of learning. Lastly, a rare-event setting where the system ispartially known and learned is theoretically analyzed in Lam (2012), but only confined to the setting oflevel-crossing of random walks.
The remainder of this paper is as follows. Section 2 reviews the IS technique including in particular
the Gaussian input case. Section 3 presents our precise setting and the proposed procedure. Section 4shows numerical results and discusses some extensions.
2BASIC BACKGROUND

Suppose that we are interested in estimating the probability of hitting a set R, i.e, p = P(X ∈ R). IS
replaces the original probability measure P with P̃, such that the likelihood ratio P(X)/P̃(X), given by theRadon-Nikodym derivative between P(·) and P̃(·), is well defined on the set R. To adjust for the bias ofthis change, IS outputs
P(X)
.
(1)
Z = I(X ∈ R)
P̃(X)
To measure the efficiency of an IS scheme, we introduce a rarity parameter, say n, that parametrizes
the rare-event probability pn such that pn → 0 as n → ∞. Note that since the probability of interest is small,
one should focus on the relative error of the Monte Carlo estimator with respect to the magnitude of thisprobability. To this end, we call an IS estimator Zn for pn asymptotically efficient if
log Ẽ[Zn2 ]= 2.
n→∞ log Ẽ[Zn ]
lim

(2)

The notion (2) is equivalent to saying that Ẽ[Zn2 ]/Ẽ[Zn ]2 is at most polynomially growing in n. This ensuresthat the second moment, or the variance, does not explode exponentially relative to the probability ofinterest as n increases, thus preventing an exponentially large number of simulation replications to achievea given relative accuracy. For further details, see, e.g., Asmussen and Glynn (2007).
As an example that will be used in this paper, consider a standard Gaussian random vector X ∼ N(0, I)
and a rare-event set R. We recall the notion of a dominating point a∗ of R, i.e., a point a∗ ∈ R suchthat R ⊆ {x : a∗ 0 (x − a∗ ) ≥ 0}. Suppose that R is convex, then this condition is equivalent to saying thata∗ optimizes minx∈R kxk2 , as 2a∗ 0 (x − a∗ ) = ∇(kxk2 )0 (x − a∗ ) ≥ 0 is precisely the first order condition of
optimality for the convex program minx∈R kxk2 . Moreover, since we consider standard Gaussian X, wehave a∗ = argmax{φ (x) : x ∈ R}, where φ denotes the density of standard Gaussian distribution.
Suppose that a∗ is big, so that a∗ → ∞ and R becomes rare. A commonly considered IS scheme is touse exponential tilting, triggered by a mean shift of X from 0 to a∗ , i.e, we use the estimator
Z = I(X ∈ R)

φ (X).
φ (X − a∗ )(3)

The second moment is given by
"
#
2
h
i
h
i
φ (X)
2
−kXk2 +kX−a∗ k2−ka∗ k2−2a∗0 (X−a∗ )Ẽ[Z ] = Ẽ
;
X
∈
R
=
Ẽ
e
;
X
∈
R
=
e
Ẽ
e
;
X
∈
R
.
φ (X − a∗ )Since a∗ is assumed a dominating point, we have a∗ 0 (x − a∗ ) ≥ 0 for any x ∈ R, and thus the above is∗ 2bounded from above by e−ka k .
∗ 2
Since Ẽ[Z]2 = P(X ∈ R) = poly(a∗ )e−ka k , where poly(a) denotes a polynomial in a∗ , the scheme (3)is asymptotically efficient as defined by (2).

Huang, Lam, and Zhao
3

IMPORTANCE SAMPLING FOR UNKOWN RARE-EVENT SETS

Motivated by the discussion and example given in the introduction, we consider the situation where therare-event set is unknown. To do so, we refine our setting in Section 2. We assume R is unknown, but,that we know R is convex, so that the rare-event set has an unknown dominating point a∗ . Moreover,assume it is known that R ⊂ [0, B]m for some B > 0 (these could be the physical boundary of the given
problem). Given a point x ∈ Rm , we have the capability to observe whether x ∈ R, but this action would
incur a cost that is equivalent to running a simulation copy of I(X ∈ R).
The basic challenge is to attempt to learn R so that we can then apply IS effectively. Given abudget n on the number of evaluations, we investigate a two-stage IS scheme as follows. In the firststage, we use a space filling design (e.g. uniform sampling) over the design space [0, B]m and collect1
n1 observations (Xi , I(Xi ∈ R))ni=1. In the second stage, we construct an IS distribution by shifting
the mean of the standard Gaussian distribution to â = argminXi ,i=1,...,n1 {kXi k2 : Xi ∈ R}, or equivalently
argmaxXi ,i=1,...,n1 {φ (Xi ) : Xi ∈ R}. In other words, each IS output in the second stage is
Z = I(X ∈ R)

φ (X).
φ (X − â)(4)

We use n2 simulation replication in the second stage to obtain a sample mean as our final estimate ofP(X ∈ R). In the case that â does not exist, i.e., none of the first-stage samples fall into R, we perform a“back-up” procedure in the second stage that has a second moment of, say, M (this can be, e.g., a uniform
sampling on the whole space). For convenience of analysis, we use n1 = n2 . We summarize our procedurein Algorithm 1.
Algorithm 1: The two-stage IS scheme for unknown rare-event sets.
Input: Sample budget n1 , n2 .
Output: IS estimator Z.
Stage 1:
1 Use space filling design (e.g. uniform sampling) over the design space [0, B]m to collect n1
1
observations (Xi , I(Xi ∈ R))ni=1;
2 Find approximate dominating point by solving â = argminXi ,i=1,...,n {kXi k2 : Xi ∈ R};
1
Stage 2:
3 If â exists do
4

φ (X)Obtain the IS estimator as the sample average of n2 copies of Z = I(X ∈ R) φ (X−â) ;5 Else do

A back-up procedure (e.g. use space filling design to collect n2 samples and output theassociated IS estimate);
7 End;
6

3.1 Error When Using an Approximate Dominating Point
The second-stage IS is precisely the classical scheme introduced in Section 2, except that we replace theunknown a∗ with an estimator â. The following approximates the error from such an estimation:
Theorem 1 (Error from using approximate dominating point) Consider estimating P(X ∈ R) where R ∈ Rm
is a convex set and X ∼ N(0, I). Denote a∗ as a dominating point of R, and d = kâ − a∗ k for a given point∗ 22â. Suppose we use the IS scheme (4). Then its second moment satisfies Ẽ[Z 2 ] ≤ e−ka k +d .

Huang, Lam, and Zhao
Proof.

Consider the second moment of Z given by
Ẽ[Z 2 ]2

2= Ẽ[e−kXk +kX−âk ; X ∈ R]
∗ 2

0

∗ 2

2

∗ 2

∗ 2

∗ 0 ∗

∗ 2

∗ 2

∗ 0

∗

∗ 2

∗ 2

∗ 0

∗

∗ 2

∗ 2

= e−ka k Ẽ[ekâk +ka k −2X â ; X ∈ R]
∗ 0∗

= e−ka k +kâ−a k Ẽ[e−2(X−a ) a −2(X−a ) (â−a ) ; X ∈ R]
≤ e−ka k +kâ−a k Ẽ[e−2(X−a ) (â−a ) ; X ∈ R]
since (x − a∗ )0 a∗ ≥ 0 for any x ∈ R by the definition of dominating point≤ e−ka k +kâ−a k Ẽ[e−2(X−a ) (â−a ) ]= e−ka k +kâ−a k∗ 0∗where the last step follows by noting that X − a∗ ∼ N(â − a∗ , I) under P̃, and Ẽ[e−2(X−a ) (â−a ) ] is the
∗ 0
∗
∗ 2
∗ 2
moment generating function of X − a∗ at −2(â − a∗ ), so that Ẽ[e−2(â−a ) (X−a ) ] = e−2kâ−a k +4kâ−a k /2 = 1.
This concludes the theorem.
Note that Theorem 1 does not assume how â is obtained. It is applicable for any â. An implication of
the theorem is that the IS scheme (4) using â is asymptotically efficient, except for an extra multiplicative2
error ed in the ratio Ẽ[Z 2 ]/Ẽ[Z]2 that, depending on how close â is to a∗ , could potentially be significant.
3.2 Efficiency of the Two-stage Procedure
Intuitively, if we use large enough n1 in the first stage, we could get a good estimate of a∗ , i.e., d = kâ − a∗ k
is small probabilistically. To further analyze this error, let us denoteα=Area(R)
Area([0, B]m )We recall that the first-stage observations of X are sampled uniformly on [0, B]m and used to obtainâ. Consider P(kâ − a∗ k2 > y). Note that kxk2 = kx − a∗ k2 + ka∗ k2 + 2(x − a∗ )0 a∗ ≥ kx − a∗ k2 + ka∗ k2 =kx − a∗ k2 + η 2 for any x ∈ R, where we denote η = ka∗ k. We haveP(kâ − a∗ k2 > y) ≤ P(kâk2 > y + η 2 ) since â ∈ R, if exists= P(kXi k2 > y + η 2 or Xi ∈/ R)n1by our construction of â as the minimizer of kXi k2 among Xi ∈ R= (1 − αP(kXi k2 − η 2 ≤ y|Xi ∈ R))n1(5)

To proceed, we make several simplifying assumptions. We let m = 2. Moreover, we consider an
additional assumption thaty1
Area({x ∈ R : kxk2 − η 2 ≤ y1 })≥kArea({x ∈ R : kxk2 − η 2 ≤ y2 })y2for y1 < y2 running from 0 to a constant c > 0, where k is a constant such that 0 < k ≤ 1. This assumption
holds for m = 2 if the rare-event set roughly looks like a cone with rays generated from a∗ . This includes,in particular, shapes in the form of a semi-disk.
Under the above assumptions, from (5), we havey n1
P(kâ − a∗ k2 > y) ≤ 1 − αkcHuang, Lam, and Zhao
∗ 22Consider the second moment of the overall estimator of the two-stage procedure, which is e−ka k E[ed ]∗ 22where E[·] is under the first-stage uniform sampling. Here we define ed = eka k M if â and hence d does
not exist (Recall that if â does not exist, then we use a back-up procedure giving a second moment M ineach second-stage replication). Now, denote F̄(y) = P(d 2 > y), and consider
2

E[ed ] = −Z c2ey d F̄(y) + (1 − α)n1 eη M

02

≤ 1 + (1 − αk)n1 (eη M − ec ) +Z c

1 − αk0y n1 y
e dy
c

(6)

Consider the termZ c

1 − αk01
y n1 ye dy = c(1 − αkz)n1 ecz dz
c
0

Z

Z 1

=c

en1 log(1−αkz)+cz dz

0

≤c

Z 1

e−bn1 αkz+cz dz(7)0

for some constant b > 0, where the last step follows since we can find such a b such that log(1−αkz) ≤ −bαkz
for all 0 ≤ z ≤ 1.
Now consider n1 = ω(1/α), or more precisely n1 = ρ/α where ρ → ∞. Then (7) is equal toZ 1ce−bkρz+cz dz =0c(1 − e−(bkρ−c) )bkρ − cContinuing with the same choice of n1 = ρ/α, the second term in (6) becomes2(1 − αk)ρ/α (eη M − ec )Thus, if we use n2 = n1 = ρ/α replications in the second stage, we get an overall second moment of
the final estimator bounded from above by2c(1 − e−(bkρ−c) ) + (1 − αk)ρ/α (eη M − ec ) + 1−η 2 bkρ−c

e

ρ/α

(8)

3.3 Comparison with Pure Uniform SamplingTo give some insights on (8), we compare our approach with a naive uniform sampling scheme. Since
we have assumed that R ⊂ [0, B]m , using uniform sampling as IS is feasible. To compare fairly with our
two-stage procedure, suppose we use a sample size n = n1 + n2 for this uniform IS. The second moment
for each replication is given by
φ (X) 2Ẽ () ;X ∈ R
1/V
where V = Area([0, B]m ). This is equal toZ2
V Ẽ φ (X) ; X ∈ R = V
2

x∈R

if we let m = 2.

φ (x)2 dx ≈V −η 2eη2(9)Huang, Lam, and Zhao
For further comparison, we set n1 = n2 = ρ/α, so that the total replication size is n = 2ρ/α. The
second moment becomes approximately
V α −η 2e(10)2η 2 ρ
On the other hand, suppose that, in our two-stage procedure, the back-up algorithm is a pure uniform
sampling, then M is approximately replaced by (9), so that the overall second moment is approximately
bounded byc−(bkρ−c) ) + (1 − αk)ρ/α (V /η 2 − ec ) + 1
2 bkρ−c (1 − e
(11)
e−η
ρ/α
Now consider the case where α is bounded away from 1, i.e., the rare-event set does not occupy thewhole space, and that c > 0 is fixed, i.e., the volume of the rare-event set does not grow (This can becrelaxed to letting c grow slowly). Since bkρ−c(1 − e−(bkρ−c) ) is bounded, and (1 − αk)1/α ≤ e−l for someconstant l > 0, (11) is bounded by−lρ V /η 2 )2 α( f + e
(12)
e−η
ρ
for some constant f > 0.
Thus, comparing (10) and (12), we see that when V /η 2 is large or goes to ∞, the second momentprovided by the two-stage procedure is smaller, and thus it is preferable. Note that V /η 2 occurs when thespace in consideration is large compared to the position of the rare-event set. Intuitively, in this situation
most of the uniformly sampled points will miss the rare-event set. Thus a pure uniform sampling works
poorly. On the other hand, the two-stage approach makes use of a concentration-type behavior comingfrom solving the optimization based on the few hitting samples in the first stage, in order to locate thedominating point. This operation makes the overall performance of the two-stage approach better.
4

NUMERICAL EXPERIMENTS

This section uses some simple problems to illustrate the performance of our proposed approach and tocompare with uniform IS scheme. We test two types of problems. In the first type of problems, we assume
that the rare-event set is convex. In this case, the assumption for Theorem 1 holds and the proposed
approach has the analyzed efficiency. The second type of problems considers non-convex rare-event set
that violates our assumption, which serves to illustrate potential extensions of our proposed approach.
4.1 Convex Rare-event Set Example
We consider the following problem setting. The design vector x has dimension m = 2 and the design spaceis [0, 5]m , i.e. B = 5 according to the notation we used in previous sections. Note that X follows standard
Gaussian distribution. We set the rare-event set to be R = {x : (x1 − 5)2 + (x2 − 5)2 ≤ 4}. We compare theperformance of the proposed two-stage approach and the uniform IS in two sets of experiments.
In the first set of experiments, we vary the number of total sample budget from n = 1000 to n = 50, 000.
For each selection of n, we use n1 = n2 = n/2 samples in the two-stage approach. Results of theseexperiments are shown in Figures 1 and 2. Figure 1 presents the final probability estimates from thetwo approaches with different sample budgets. We observe that the two-stage approach (red solid line) isslightly more stable than the uniform IS (black dash line), as the two-stage approach has less fluctuation
when the budget is greater than 30,000. This observation is confirmed in Figure 2, where the confidence
interval half-width of the two-stage approach (red solid line) is shown to be always better than the uniform
IS (black dash line). Roughly speaking, the two-stage approach appears to require only half of the budgetof the uniform IS for achieving the same level of accuracy.
We specifically look at the case where n = 50, 000. Figure 3 shows a trajectory of the confidence
interval half-width of the two approaches in a single sample path. The blue horizontal line shows the final

Huang, Lam, and Zhao

Figure 1: Probability estimates under different Figure 2: Confidence interval half-widths undersample budgets.
different sample budgets.

Figure 3: Confidence interval half-width in the IS
stage when n = 50, 000. The blue horizontal line
shows the final confidence interval half-width of
the uniform IS.
half-width of the uniform IS for easy comparison. It indicates that at roughly 10,000 samples in the second
stage (35,000 samples in total), the two-stage approach obtains the same half-width as the uniform IS using
all 50,0000 samples. Overall, the two-stage approach performs visually better compared to the uniform IS.
Our second set of experiments serves to illustrate the analysis in Section 3.3. We consider a similar
problem setting, but with a growing design space size (so that V /η 2 increases). For the design space[0, B]2 , we have B vary from B = 4 to B = 7 with 0.1 increments. We consider the rare-event set
R = {x : (x1 − B)2 + (x2 − B)2 ≤ 4} that has a stationary size as B changes. For each value of B, we use
n = 50, 000 samples for both approaches, where the two-stage approach uses n1 = n2 = n/2. Figure 4 showsthat as the value of B increases, the probability estimate decays exponentially, and the two approaches
give very similar values. However, in Figure 5, we observe that the width of confidence interval from thetwo-stage approach (red solid line) is smaller than the uniform approach (black dash line) for all B values.
Figure 6 further presents the ratio between the width of confidence interval for the two-stage approach andthe uniform approach. The ratio decays roughly linearly, which indicates that as the size of the design
space increases, the two-stage approach has more advantages.

Huang, Lam, and Zhao

Figure 4: Probability estimates with different Figure 5: Confidence interval half-widths withdesign space sizes.
different design space sizes.

Figure 6: Ratios of the confidence interval halfwidths from the two-stage approach over theuniform IS.
These experiments indicate that the two-stage approach is efficient for problems with convex rare-event
set compared to the uniform IS. Moreover, it shows more benefits when the design space is big compared to
the location of the rare-event set, which is consistent with the findings in Section 3.3. Next, we investigate
an extension of our approach to non-convex rare-event sets.
4.2 Extension to Non-convex Rare-event Sets
We consider a heuristic approach to generalize our scheme to rare-event set R that is more complex and has
multiple dominating points, say a∗1 , ..., a∗p , such that R ⊆ ∪i {x : a∗i 0 (x − a∗i ) ≥ 0} and ai ∈ R, i = 1, ..., p. In
the first stage, we attempt to find a collection of approximate dominating points, denoted â1 , ..., â p̂ , where
we note that p̂ may not necessarily equal p. In the second stage, we build an IS distribution that consists ofa p̂-mixture of Gaussian distributions, where each mixture component has mean âi and identity covariance.
In the first stage, we can find the approximate dominating points sequentially. We first find â1 =argminxi ,i=1...,n0 {kxi k2 : xi ∈ R}. Then we use a “cutting plane” to find â2 = argminxi ,i=1...,n0 {kxi k2 : xi ∈R, â01 (xi − â1 ) < 0}. We keep adding constraints to these optimization problems in a similar fashion, until
we cover all the rare-event points.

Huang, Lam, and Zhao

Figure 7: Probability estimates under different Figure 8: Confidence interval half-widths undersample budgets.
different sample budgets.
We use similar problem settings as in Section 4.1 to compare the performances of the two-stage approachwith the uniform IS. We modify the rare-event set in Section 4.1 to R = {x : (x1 − 5)2 + (x2 − 5)2 ≤ 4} ∪ {x :(x1 − 5)2 + (x2 − 1)2 ≤ 1}. Note that this set can no longer be contained by any half space in the form of{x : a0 (x − a)}.
We again consider different sample budgets from n = 1, 000 to n = 50, 000. Figure 7 presents theprobability estimates of the two approaches with different sample budgets. In this figure, it is not that
obvious which approach is more stable. Nonetheless, Figure 8 shows that the two-stage approach (redsolid line) is better than the uniform approach (black dash line) in terms of shorter confidence intervalhalf-widths. These results are quite similar to the convex set case, and hints that the two-stage approachshows better performances for such types of rare-event problems as well.
5

CONCLUSION

We have considered the use of IS to estimate rare-event probabilities that lack structural information thatis typically needed to design efficient schemes. We have considered a two-stage approach where thefirst, exploration, stage generates a space-filling set of points to learn the rare event of interest, and thesecond, exploitation, stage estimates the dominating point to guide an exponential tilting. Our derivation
indicates that our two-stage IS estimator is efficient when, roughly speaking, the volume of the space islarge compared to the position of the rare-event set, when compared to an IS based on a naive uniform
sampling. This phenomenon is confirmed in two sets of numerical experiments. In future studies, we will
extend our efficiency analysis to more general settings and compare with other potential methods.
ACKNOWLEDGMENTS
We gratefully acknowledge support from the National Science Foundation under grants CMMI-1542020
and CAREER CMMI-1653339/1834710.
REFERENCES
Adler, R. J., J. H. Blanchet, J. Liu et al. 2012. “Efficient Monte Carlo for High Excursions of Gaussian
Random Fields”. The Annals of Applied Probability 22(3):1167–1214.
Asmussen, S., and P. W. Glynn. 2007. Stochastic Simulation: Algorithms and Analysis, Volume 57. New
York: Springer Science & Business Media.
Au, S., and J. L. Beck. 1999. “A New Adaptive Importance Sampling Scheme for Reliability Calculations”.
Structural Safety 21(2):135–158.

Huang, Lam, and Zhao
Bengio, Y., and J.-S. Senécal. 2008. “Adaptive Importance Sampling to Accelerate Training of A Neural
Probabilistic Language Model”. IEEE Transactions on Neural Networks 19(4):713–722.
Blanchet, J., and C. Li. 2011. “Efficient Simulation for the Maximum of Infinite Horizon Discrete-time
Gaussian Processes”. Journal of Applied Probability 48(2):467–489.
Blanchet, J., and H. Lam. 2012. “State-dependent Importance Sampling for Rare-Event Simulation: An
Overview and Recent Advances”. Surveys in Operations Research and Management Science 17(1):38–
59.
Blanchet, J., P. Glynn, and H. Lam. 2009. “Rare Event Simulation for A Slotted Time M/G/s Model”.
Queueing Systems 63(1-4):33.
Botev, Z. I., P. LEcuyer, and B. Tuffin. 2013. “Markov Chain Importance Sampling with Applications toRare Event Probability Estimation”. Statistics and Computing 23(2):271–285.
Bourinet, J.-M., F. Deheeger, and M. Lemaire. 2011. “Assessing Small Failure Probabilities by Combined
Subset Simulation and Support Vector Machines”. Structural Safety 33(6):343–353.
Bucklew, J. 2013. Introduction to Rare Event Simulation. New York: Springer Science & Business Media.
Chan, J. C., and D. P. Kroese. 2012. “Improved Cross-entropy Method for Estimation”. Statistics andComputing 22(5):1031–1040.
Das, P., and Y. Zheng. 2000. “Cumulative Formation of Response Surface and Its Use in Reliability
Analysis”. Probabilistic Engineering Mechanics 15(4):309–315.
De Boer, P.-T., D. P. Kroese, S. Mannor, and R. Y. Rubinstein. 2005. “A Tutorial on the Cross-entropy
Method”. Annals of Operations Research 134(1):19–67.
Dieker, A. B., and M. Mandjes. 2006. “Fast Simulation of Overflow Probabilities in A Queue with Gaussian
Input”. ACM Transactions on Modeling and Computer Simulation (TOMACS) 16(2):119–151.
Dubourg, V., B. Sudret, and F. Deheeger. 2013. “Metamodel-based Importance Sampling for Structural
Reliability Analysis”. Probabilistic Engineering Mechanics 33:47–57.
Dupuis, P., A. D. Sezer, H. Wang et al. 2007. “Dynamic Importance Sampling for Queueing Networks”.
The Annals of Applied Probability 17(4):1306–1346.
Echard, B., N. Gayton, and M. Lemaire. 2011. “AK-MCS: An Active Learning Reliability Method Combining
Kriging and Monte Carlo Simulation”. Structural Safety 33(2):145–154.
Frank, J., S. Mannor, and D. Precup. 2008. “Reinforcement Learning in the Presence of Rare Events”.
In Proceedings of the 25th International Conference on Machine Learning, edited by McCallum, A.
et al., 336–343. New York: ACM.
Glasserman, P. 2013. Monte Carlo Methods in Financial Engineering, Volume 53. New York: Springer
Science & Business Media.
Glasserman, P., P. Heidelberger, and P. Shahabuddin. 1999. “Asymptotically Optimal Importance Sampling
and Stratification for Pricing Path-Dependent Options”. Mathematical Finance 9(2):117–152.
Glynn, P. W., and D. L. Iglehart. 1989. “Importance Sampling for Stochastic Simulations”. Management
Science 35(11):1367–1392.
Grabaskas, D., M. K. Nakayama, R. Denning, and T. Aldemir. 2016. “Advantages of Variance Reduction
Techniques in Establishing Confidence Intervals for Quantiles”. Reliability Engineering & System
Safety 149:187–203.
Grace, A. W., D. P. Kroese, and W. Sandmann. 2014. “Automated State-dependent Importance Samplingfor Markov Jump Processes via Sampling from the Zero-Variance Distribution”. Journal of Applied
Probability 51(3):741–755.
Heidelberger, P. 1995. “Fast Simulation of Rare Events in Queueing and Reliability Models”. ACM
Transactions on Modeling and Computer Simulation (TOMACS) 5(1):43–85.
Huang, Z., H. Lam, D. J. LeBlanc, and D. Zhao. 2017a. “Accelerated Evaluation of Automated Vehicles
Using Piecewise Mixture Models”. IEEE Transactions on Intelligent Transportation Systems: 1–11.

Huang, Lam, and Zhao
Huang, Z., H. Lam, and D. Zhao. 2017b. “An Accelerated Testing Approach for Automated Vehicles with
Background Traffic Described by Joint Distributions”. In Proceedings of 2017 IEEE 20th International
Conference on Intelligent Transportation Systems (ITSC), 933–938. Piscataway, New Jersey: IEEE.
Huang, Z., Y. Guo, H. Lam, and D. Zhao. 2018. “A Versatile Approach to Evaluating and Testing Automated
Vehicles based on Kernel Methods”.
Hurtado, J. E. 2004. “An Examination of Methods for Approximating Implicit Limit State Functions fromthe Viewpoint of Statistical Learning Theory”. Structural Safety 26(3):271–293.
Juneja, S., and P. Shahabuddin. 2006. “Rare-event Simulation Techniques: An Introduction and Recent
Advances”. Handbooks in Operations Research and Management Science 13:291–350.
Kanj, R., R. Joshi, and S. Nassif. 2006. “Mixture Importance Sampling and Its Application to the Analysis
of SRAM Designs in the Presence of Rare Failure Events”. In Proceedings of Design Automation
Conference, 2006 43rd ACM/IEEE, 69–72. Piscataway, New Jersey: IEEE.
Lam, H. 2012. “Efficient Importance Sampling under Partial Information”. In Proceedings of the 2012
Winter Simulation Conference (WSC), edited by Laroque, C. et al., 1–12. New Jersey: IEEE.
OHagan, A. 2006. “Bayesian Analysis of Computer Code Outputs: A Tutorial”. Reliability Engineering& System Safety 91(10-11):1290–1300.
Papadrakakis, M., and N. D. Lagaros. 2002. “Reliability-based Structural Optimization Using Neural
Networks and Monte Carlo Simulation”. Computer Methods in Applied Mechanics and Engineering 191(32):3491–3507.
Rubinstein, R. Y., and D. P. Kroese. 2004. The Cross-Entropy Method: A Unified Approach to Combinatorial
Optimization, Monte-Carlo Simulation and Machine Learning. New York: Springer.
Rubinstein, R. Y., and D. P. Kroese. 2016. Simulation and the Monte Carlo Method, Volume 10. New
Jersey: John Wiley & Sons.
Sadowsky, J. S., and J. A. Bucklew. 1990. “On Large Deviations Theory and Asymptotically Efficient
Monte Carlo Estimation”. IEEE transactions on Information Theory 36(3):579–588.
Shelton, C. R. 2001. Importance Sampling for Reinforcement Learning with Multiple Objectives. Ph.D.
thesis, Massachusetts Institute of Technology, Cambridge, MA.
Siegmund, D. 1976. “Importance Sampling in the Monte Carlo Study of Sequential Tests”. The Annals of
Statistics: 673–684.
Zhao, D., H. Lam, H. Peng, S. Bao, D. J. LeBlanc, K. Nobukawa, and C. S. Pan. 2017. “Accelerated
Evaluation of Automated Vehicles Safety in Lane-change Scenarios Based on Importance Sampling
Techniques”. IEEE Transactions on Intelligent Transportation Systems 18(3):595–607.
Zhao, Y.-G., and T. Ono. 1999. “A General Procedure for First/second-order Reliability Method(FORM/SORM)”. Structural Safety 21(2):95–112.
AUTHOR BIOGRAPHIESZHIYUAN HUANG is a third-year Ph.D. student in the Department of Industrial and Operations Engineering at the University of Michigan, Ann Arbor. His research interests include simulation and stochastic
optimization. His email address is zhyhuang@umich.edu.
HENRY LAM is an Associate Professor in the Department of Industrial Engineering and Operations Research at Columbia University. His research focuses on Monte Carlo simulation, uncertainty quantification,risk analysis, and stochastic and robust optimization. His email address is khl2114@columbia.edu.
DING ZHAO is an Assistant Professor in the Department of Mechanical Engineering at Carnegie Mellon University. His research focuses on Connected and Automated Vehicles (CAVs) using synthesized
approaches rooted in advanced statistics, modeling, optimization, dynamic control, and big data analysis.
His email address is dingzhao@cmu.edu.