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

DESIGNING IMPORTANCE SAMPLERS TO SIMULATE MACHINE LEARNING
PREDICTORS VIA OPTIMIZATION
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
We study the problem of designing good importance sampling (IS) schemes to simulate the probability
that a sophisticated predictor, built for instance from an off-the-shelf machine learning toolbox, gives aprediction that exceeds a large threshold. This problem is motivated as a step towards building good learningalgorithms that takes into account the extremal risks of the prediction. We provide a framework to design
IS for two common machine learning models, namely random forest and a basic neural network. Our
approach utilizes some available mathematical programming formulations to optimize over these models
and a simple “cutting plane” idea to look for dominating points under Gaussian input distributions.
1

INTRODUCTION

We consider the problem of estimating P(g(X) > γ) via simulation, where X ∈ Rd is a random input andγ ∈ R is a high threshold, so that the probability is small and leads to a rare-event simulation problem.
We are interested particularly in g(·) that is the response of a sophisticated predictor, built for instancefrom an off-the-shelf machine learning toolbox. This problem is motivated as a step towards building goodlearning methods that take into account the extremal risks of machine learning prediction.
To be more concrete, suppose we are interested in training the parameter θ of a predictor, say Yθ (X), by
minimizing a risk criterion E[L(Yθ (X),Y (X))], where Y (·) is the true response function. Suppose that thisloss function L exerts a big value in some “hidden region” where the occurrence of X is rare. The data inthis case may not reveal satisfactorily this hidden risk. An approach to learn a good risk-conscious predictor
in this situation is to fit the data using some probability distribution, and use it to train the predictor. A
judicious choice of the distribution allows one to generate Monte Carlo samples, and to use importancesampling (IS) to populate more samples in the rare-event set. Of course, this approach requires various
aspects of considerations, both statistically and practically. However, one step towards handling such typeof problems involves how to speed up simulation for a machine learning model, which is precisely what
we address in this paper.

Huang, Lam, and Zhao
As a motivating example, autonomous vehicle (AV) design often uses predictions of the movementsof surrounding cars to decide its own behavior (e.g. Jain et al. 2015). Predictions that seriously deviate
from the reality can post substantial danger and catastrophic failure. However, unusual surrounding driving
movements that lead to serious deviations of the predictions are rare. In this case, one can represent thedriving environment as a stochastic model (the X), which provides inputs for the prediction model (theYθ (X)) (e.g. Wang et al. 2017; Huang et al. 2017), and a risk-sensitive predictive training could involverare-event simulation on functionals associated with Yθ (X).
2

PROBLEM SETTING

We state more precisely our problem setting. Suppose we are given a prediction model g(·), with the input
X ∈ Rd and the output g(X) ∈ R. Suppose that the input follows a standard Gaussian distribution, i.e,
X ∼ N(0, Id ), where Id is the d × d identity matrix. We want to estimate the probability p = P(g(X) ≥ γ),
where γ ∈ R is a threshold that triggers a rare event. We note that the Gaussian assumption can be relaxedwithout much difficulty in our framework to, for instance, mixtures of Gaussians, which can expand ourscope of applicability. For simplicity, however, we will concentrate on the standard Gaussian distribution
in this paper.
When p is small, estimation using crude Monte Carlo is challenging due to large variance. A common
approach to speed up simulation in such contexts is to use IS (see, e.g. the surveys Bucklew 2013; Asmussen
and Glynn 2007; Rubinstein and Kroese 2016; Glasserman 2013; Juneja and Shahabuddin 2006; Blanchetand Lam 2012, among others). Suppose X has a density f . The basic idea of IS is to change the sampling
distribution to say f˜, and output
f (X)
,
(1)
Z = I(g(X) ≥ γ)f˜(X)where X is now sampled from f˜. This output is unbiased if f is absolutely continuous with respect to f˜over the rare-event set {x : g(x) ≥ γ}. Moreover, by choosing f˜ appropriately, one can substantially reduce
the simulation variance.
In the case of Gaussian input distributions, finding a good f˜ is particularly handy and one approachto devise good IS distributions uses the notion of so-called dominating point. Starting from a standard
Gaussian distribution, a well-known IS scheme is as follows (see, e.g., Sadowsky and Bucklew 1990;
Dieker and Mandjes 2006; Blanchet and Li 2011; Adler et al. 2012):
1. If there exists a dominating point a, i.e., a point a such that a = arg minx {kxk2 : g(x) ≥ γ} and{x : g(x) ≥ γ} ⊆ {x : a0 (x − a) ≥ 0}, then we use a Gaussian distribution with mean at a as the IS
distribution f˜.
2. If we can split {x : g(x) ≥ γ} into R1 , ..., Rr , and for each Ri , i = 1, ..., r there exists a dominating
point ai such that ai = arg minx {kxk2 : x ∈ Ri } and Ri ⊆ {x : a0i (x − ai ) ≥ 0}, then we use a Gaussian
mixture distribution with r components as the IS distribution f˜, where the ith component has meanai .
Note that the quantity arg minx {kxk2 : g(x) ≥ γ} is equivalent to arg maxx {φ (x) : g(x) ≥ γ}, where φ isthe standard Gaussian density. The condition 2a0 (x − a) ≥ 0 is the first order condition of optimality forthe optimization minx kxk2 over a convex set for x, and so the first situation above occurs if {x : g(x) ≥ γ}
is a convex set (though it can be more general). The proposals above guarantee the so-called asymptotic
efficiency or logarithmic efficiency (e.g., Asmussen and Glynn 2007) of the IS, meaning that the relative
error, measured by the ratio of standard deviation of a single IS output over the probability of interest,
grows at most polynomially in the location of the rare-event set.
Given the proposal above, one can follow Algorithm 1 to obtain the dominating points a1 , ..., ar to build
an efficient IS distribution in our prediction model context. The procedure uses a sequential “cutting plane”
approach to exhaustively look for all dominating points, by reducing the search space at each iteration via
taking away the regions covered by the existing dominating points. The set A in the procedure serves tostore the dominating points we have found throughout the procedure. At the end of the procedure, one
obtain a set A that contains all the dominating points a1 , ..., ar .

Huang, Lam, and Zhao
Algorithm 1: Procedure to find all dominating points for the set {x : g(x) ≥ γ}.
Input: Prediction model g(x), threshold γ.
Output: Dominating-points set A.
1 Start with A = 0;/
2 While {x : g(x) ≥ γ, a0i (x − ai ) < 0, ∀ai ∈ A} 6= 0
/ do
3Find a dominating point a by solving the optimization problem
a = arg min kxk2
x(2)

s.t. g(x) ≥ γ
a0i (x − ai ) < 0, for ∀ai ∈ A
and update A ← A ∪ {a};4 End

To apply Algorithm 1 to find all dominating points, the key is to be able to solve the optimization
problems (2). We will investigate this for two popular prediction models, random forest and neural network(NN). Section 3 studies the tractable formulations using some recent work on optimization over thesemodels, while Section 4 applies them to design IS schemes and demonstrate some numerical results.
3

TRACTABLE FORMULATION FOR PREDICTION MODELSWe discuss how to formulate the optimization problems in Algorithm 1 as a mixed integer program (MIP)
with quadratic objective function and linear constraints. Sections 3.1 and 3.2 focus on a basic NN andrandom forest respectively. Section 3.3 then further discusses solving larger-scale problems using a simple
bisection method that leverages existing techniques in Benders decompositions in the random forest case.
3.1 Tractable Formulation for Neural Network
We consider a NN with L layers and for each layer, the number of neurons is n1 , ..., nL . (As we focus on
output g(x) ∈ R, we have nL = 1.) We use rectified linear unit (ReLU) for each neuron, i.e., the activation
function for each neuron is max{0, x}. The input of the jth neuron in layer i is weighted from the outputof the previous layer by a vector wij ∈ Rni−1 and is added by a bias bij ∈ R. We use si ∈ Rni , i = 1, ..., L torepresent the output of the ith layer. At the ith layer, given the output from the (i − 1)th layer si−1 , weThave si ∈ Rni , where the jth element of si is given by sij = max{wij si−1 + bij , 0}. For further details onsuch type of NNs, see, e.g., Goodfellow et al. (2016).
Here we reformulate (2) by replacing g(x) ≥ γ with constraints that represent the structure of the NN:sL ≥ γTsij = max{wij si−1 + bij , 0}, i = 1, ..., L, j = 1, ..., nis0 = x.
Tjeng and Tedrake (2017) discussed one approach to express ReLU as a mixed integer model, by
introducing auxiliary binary variables z. For y = max(x, 0), if we have l ≤ x ≤ u, where l ≤ 0 and u ≥ 0 are
the smallest and largest possible values of x, then we set z = I(x ≥ 0). We can reformulate y = max(x, 0)
as a set of linear constraints: y ≤ x − l(1 − z); y ≥ x; y ≤ uz; y ≥ 0; z ∈ {0, 1}.
Using this approach, we rewrite problem (2) with A = 0/ as
minx,s0 ,...,sL ,z1 ...,zL

kxk2

(3)

Huang, Lam, and Zhao
s.t. sL ≥ γsi ≤ Wi T si−1 + bi − l(1 − zi ), i = 1, ..., Lsi ≥ Wi T si−1 + bi , i = 1, ..., Lsi ≤ uzi , i = 1, ..., Lsi ≥ 0, i = 1, ..., Lzi ∈ {0, 1}ni , i = 1, ..., L
s0 = x,where Wi = [w1i , ..., wni i ] ∈ Rni−1 ×ni , bi = [b1i , ..., bni i ]T ∈ Rni , and si = [s1i , ..., sni i ]T ∈ Rni .
Note that this optimization has a quadratic objective and linear constraints. Similarly, we can formulate(2) by adding linear constraints a0i (x − ai ) < 0, ∀ai ∈ A to (3), which arrives at the same optimization class.
Medium-size instances of these problems can be handled by standard solvers.
3.2 Tractable Formulation for Random Forest
To look for dominating points in a random forest or tree ensemble, we follow the route in Mišic (2017)that studies optimization over these models. We consider a random forest as follows. The input x has ddimensions. Suppose the model consists of T trees f1 , ..., fT . In each tree ft , we use ai, j to denote the jthunique split point for the ith dimension of the input x, such that ai,1 < ai,2 < ... < ai,Ki , where Ki is thenumber of unique split points for the ith dimension of x.
Following the notations of Mišic (2017), let leaves(t) be the set of leaves (terminal nodes) of tree t
and splits(t) be the set of splits (non-terminal nodes) of tree t. In each split s, we let left(t) be the set ofleaves that are accessible from the left branch (the query at s is true), and right(t) be the set of leaves thatare accessible from the right branch (the query at s is false). For each node s, we use V(s) ∈ {1, ..., d} todenote the dimension that participate in the node and C(s) ∈ {1, ..., KV(s) } to denote the set of values ofdimension i that participate in the split query of s (C(s) = { j} and V(s) = {i} indicates the query xi ≤ ai, j ).
T
We use λt to denote the weight of tree t (∑t=1
λt = 1). For each l ∈ leaves(t), pt,l denotes the output
for the lth leaf in tree t.
To formulate the random forest optimization as an MIP, we introduce binary decision variables zi, j andyt,l . Firstly, we havezi, j = I(xi ≤ ai, j ), i = 1, ..., d, j = 1, ..., Ki .
(4)
We then use yt,l = 1 to denote that tree t outputs the prediction value pt,l on leaf l, and yt,l = 0 otherwise.
We use z, y to represent the vectors of zi, j and yt,l respectively. For the input x, we assume that x ∈ [−B, B]dand |ai, j | ≤ B. Then (4) is represented by the following constraints
xi ≤ ai, j + 2(1 − zi, j )Bxi > ai, j − 2zi, j B.
Now we formulate (2) with A = 0/ as the following MIP
min kxk2

(5)

x,y,z

T

s.t.

∑ ∑

λt pt,l yt,l ≥ γ

t=1 l∈leaves(t)∑l∈leaves(t)yt,l = 1

Huang, Lam, and Zhao

∑ yt,l ≤ ∑ zV(s), j , ∀t ∈ {1, ..., T }, s ∈ splits(t)l∈left(t)j∈C(s)∑

yt,l ≤ 1 −l∈right(t)∑ zV(s), j , ∀t ∈ {1, ..., T }, s ∈ splits(t)j∈C(s)zi, j ≤ zi, j+1 , ∀i ∈ {1, ..., d}, j ∈ {1, ..., Ki − 1}zi, j ∈ {0, 1}, ∀i ∈ {1, ..., d}, j ∈ {1, ..., Ki }yt,l ≥ 0, ∀t ∈ {1, ..., T }, l ∈ leaves(t)xi ≤ ai, j + 2(1 − zi, j )Bxi > ai, j − 2zi, j B.
This formulation again has a quadratic objective function and linear constraints. Similarly, we can formulate(2) with A 6= 0/ by adding linear constraints a0i (x − ai ) < 0, ∀ai ∈ A to (5).
3.3 Bisection Algorithm and the Benders Decomposition for Solving Larger-scale Problems
The MIPs (5) are already tractable for small- and medium-size problems. Nonetheless, because of thespecial structure of (5), we can obtain, for larger-scale problems, an efficient algorithm based on a bisection
on a “dual” form of the problem. The latter is a linear MIP and can be solved by the Benders decomposition
considered in Mišic (2017).
We consider the “dual” form problem of (2) with A = 0:/h(η) = max g(x)x(6)

s.t. kxk2 ≤ η.
Note that (6) is non-decreasing in η. Moreover, we know that if η is above the optimal objective valuefor (2) with A = 0,
/ then (6) will output an objective value at least γ. On the other hand, if η is below that,then (6) will output a value less than γ. This motivates us to propose a bisection method that does a line
search on the smallest value of η that gives g(x) ≥ γ. This is summarized in Algorithm 2.
Note that Algorithm 2 converges to the solution of (5) as ε goes to 0. Also, the number of iterations
we need is given by Niter < log Mε / log 2. Therefore, as long as we can solve (6) efficiently, the same is true
for the dominating point problem. But, because the quadratic constraint in (6) only results in eliminatingsome of the dummy binary variables in the MIP formulation, we recover the formulation suggested inMišic (2017), but with less variables, to represent (6). This means that (6) is equivalent to
T

max
x,y,z

s.t.

λt pt,l yt,l

∑ ∑

(10)

t=1 l∈leaves(t)yt,l = 1∑l∈leaves(t)∑ yt,l ≤ ∑ zV(s), j , ∀t ∈ {1, ..., T }, s ∈ splits(t)l∈left(t)j∈C(s)∑

yt,l ≤ 1 −l∈right(t)∑ zV(s), j , ∀t ∈ {1, ..., T }, s ∈ splits(t)j∈C(s)

kxk2 ≤ η(11)zi, j ≤ zi, j+1 , ∀i ∈ {1, ..., d}, j ∈ {1, ..., Ki − 1}(12)zi, j ∈ {0, 1}, ∀i ∈ {1, ..., d}, j ∈ {1, ..., Ki }(13)

xi ≤ ai, j + 2(1 − zi, j )B

(14)

Huang, Lam, and Zhao
Algorithm 2: Bisection algorithm for solving (2) with A = 0.
/
Input: Prediction model g(x), threshold γ, tolerance parameters ε.
Output: Optimal solution a.
1 Set η1 = 0 and η2 = M, where M is large enough (e.g. M = dB2 );2 Solve
γ1 = h(η1 )andη2 − η1 ≥ ε do2Update ηm = η1 +η
2 ;
Solve

(7)

γ2 = h(η2 );(8)

γm = h(ηm );(9)3 While4
5

If γm ≥ γ do
Update η2 ← ηm ;8
Else do9
Update η1 ← ηm ;10End;
11 End;12 Update a using the optimal solution x from (8) ;6

7xi > ai, j − 2zi, j B(15)yt,l ≥ 0, ∀t ∈ {1, ..., T }, l ∈ leaves(t).
In (10), let us consider x, z as one set of variables and y as the other set of variables, where y can be further
partitioned as y = (y1 , ..., yT ) and yT consists of yt,l ’s. We observe that for any two trees t 6= t 0 , yt and yt 0
does not appear together in any constraints and are only linked through z. This observation allows us to
use the Benders decomposition as Mišic (2017) suggests.
We rewrite problem (10) as follows:
T

max
x,z

∑ λt Gt (z)

(16)

t=1

s.t. constraints (11)-(15),where Gt (z) is the optimal value of the following subproblem:
Gt (z) = max∑ λt pt,l yt,l
y,z

(17)

l∈leaves(t)s.t.

yt,l = 1∑l∈leaves(t)∑ yt,l ≤ ∑ zV(s), j , ∀t ∈ {1, ..., T }, s ∈ splits(t)l∈left(t)j∈C(s)∑

yt,l ≤ 1 −l∈right(t)∑ zV(s), j , ∀t ∈ {1, ..., T }, s ∈ splits(t)j∈C(s)

yt,l ≥ 0, ∀t ∈ {1, ..., T }, l ∈ leaves(t).
The dual of the linear problem(17) is #"min

αt ,βt ,ζt"#∑ αt,s ∑ zV(s), j + ∑ βt,s 1 − ∑ zV(s), j + ζts∈splitsj∈C(s)

s∈splitsj∈C(s)(18)

Huang, Lam, and Zhaos.t.

αt,s +∑
s:l∈left(t)∑βt,s + ζt ≥ pt,l , ∀l ∈ leaves(t)s:l∈right(t)αt,s , βt,s ≥ 0, s ∈ splits.
We use Dt to denote the set of feasible (αt , βt , ζt ) for subproblem (18) of tree t. We write (16) in termsof the dual variables (αt , βt , ζt ):
T

max
x,z,θ

(19)

∑ λt θt

t=1

"
s.t.

#

"

#

∑ αt,s ∑ zV(s), j + ∑ βt,s 1 − ∑ zV(s), j + ζts∈splitsj∈C(s)

s∈splitsj∈C(s)≥ θt , ∀(αt , βt , ζt ) ∈ Dt ,t ∈ {1, ..., T }(20)

constraints (11)-(15).
With problem (19), we can use a constraint generating scheme. We start with solving problem (19) using
a subset of D̄t ⊆ Dt in the constraint (20). We then solve (18) for each tree t to check if there exists asolution (αt , βt , ζt ) for which the constraint (20) is violated. If so, we add the constraint to (19) and solve
it again; otherwise, we obtain the optimal solution. For a quick method to solve problem (18), please refer
to Proposition 5 in Mišic (2017).
The decomposition of (2) with non-empty A follows the same route as above and the resulting formulationis very similar. The dual of the subproblem is exactly the same as (18), while the main problem is (19)with a set of additional constraints a0i (x − ai ) < 0, ∀ai ∈ A.
4

APPLICATIONS TO IMPORTANCE SAMPLING AND NUMERICAL EXPERIMENTS

In this section, we use several sets of simple experiments to illustrate the IS scheme using dominating
points generated from Algorithm 1, for the described NN and random forest. In the first set of problems,
we consider an example where there is only one dominating point. In the second set of problems, we solvea problem with multiple dominating points.
To illustrate the efficiency of the proposed IS scheme, we compare with a naive use of a uniform IS
estimator as follows. Consider a problem where X follows a distribution f (x), and the set {x : g(x) ≥ γ}
is known to lie inside [l, u]d where d is the dimension of the input variable X. The uniform IS estimatoris given by:
Zuni f orm = I(g(X) ≥ γ) f (X)(u − l)d ,where X is generated from a uniform distribution on [l, u]d . This estimator has a polynomially growing relative
efficiency as the magnitude of the dominating points grow, but the efficiency also depends significantly on
the size of the bounded set, i.e., l, u, d.
4.1 Neural Network Example: IS with A Single Dominating Point
We recall the problem setting that the input X follows a standard Gaussian distribution, and we are interested
in estimating the probability P(g(X) ≥ γ). Here g(x) represents the output of a NN prediction at x and γ
is a real-valued threshold. The NN has 3 layers with 100 neurons in each of the 2 hidden layers, and allneurons are ReLU. We consider only X in the region [0, 5]2 , so that g(x) can be thought of as being set to
0 outside this box.
The NN is trained as follows. We generate 2,601 samples using a uniform grid over the space [0, 5]2
with a mesh of 0.1 on each coordinate. For the input x = [x1 , x2 ], we use the function
y(x) = (x1 − 5)3 + (x2 − 4.5)3 + (x1 − 1)2 + x22 + 500(21)Huang, Lam, and Zhao

Figure 1: Response surface of the neural Figure 2: Rare-event set {x : g(x) ≥ γ} of
network.
the neural network.

Figure 3: Probability estimation with dif- Figure 4: 95% confidence interval halfferent numbers of samples.
width with different numbers of samples.
to generate output value as the “truth” . We obtain the dataset D = {(Xn ,Yn )} and use it to train the NNwith gradient descent. We present the response surface of g(x) in Figure 1. We use γ = 500 in this example
and the shape of the rare-event set {x : g(x) ≥ γ} in this case is presented in Figure 2. We observe that theset is roughly convex and should have a single dominating point. By solving (3), we obtain the dominating
point for the set at (3.3676, 2.6051).
We use the obtained dominating point to construct the IS estimator described in (1). In Figures 3 and4, we compare the performance of the proposed IS scheme with that of the uniform IS estimator. Figure
3 shows the estimated probabilities of the two estimators as the number of samples increases in a single
sample path. We observe that the uniform IS (black dash line) has more fluctuations than the proposed
IS (red solid line), which indicates that the proposed IS gives more stable estimates. This observation is
confirmed in Figure 4 that shows the half-width of the 95% confidence intervals of the two estimators as
the number of samples varies. Our IS appears to have shorter confidence intervals and it only takes about
12,000 samples for our IS to reach the accuracy that the uniform IS reaches with 50,000 samples.
4.2 Neural Network Example: IS with Multiple Dominating Points
We now consider true output values generated according to the function
y(x) = 10 × e 
x −5 2x −5 2
− 13
− 24

2

2

+ 10 × e−x1 −(x2 −4.5) .

(22)

We use a uniform grid over [0, 5]2 with a mesh of 0.1 on each coordinate to train a neural network with 2
hidden layers, 100 neurons in the first hidden layer and 50 neurons in the second hidden layer. All neurons
in the neural network are ReLU. The response surface of the trained model g(x) is shown in Figure 5. We

Huang, Lam, and Zhao

Figure 5: Response surface of the neural Figure 6: Rare-event set {x : g(x) ≥ γ} of
network.
the neural network.

Figure 7: Probability estimation with dif- Figure 8: 95% confidence interval halfferent numbers of samples.
width with different numbers of samples.
set γ = 8. The shape of the rare-event set is shown in Figure 6. We observe that the set now consists of
two separate regions and therefore we expect to obtain multiple dominating points. Using Algorithm 1
with the formulation in Section 3.1, we obtain two dominating points, (0.113, 4.162) and (4.187, 3.587).
We use these dominating points to construct a mixture distribution, as discussed in Section 2, as the IS
distribution.
The comparison between the proposed IS scheme and the uniform IS is shown in Figures 7 and 8.
In Figure 7, the probability estimates of the proposed IS (red solid line) appear more stable than that of
the uniform IS (black dash line). Figure 8 further shows that the confidence intervals of the proposed IS
are shorter. The efficiency of the proposed IS is roughly 4 times better than the uniform IS, consideringthat the confidence interval half-width of the proposed IS at 50,000 samples is similar to the uniform IS
at 200,000 samples.
Thus, comparing with uniform IS, the estimation accuracy of the proposed IS is better in both experiments,
demonstrating the efficiency of the IS scheme obtained in the formulation in Section 3.1.
4.3 Random Forest Example: IS with A Single Dominating Point
Consider now that g(x) represents a random forest. The random forest is trained from samples generated
uniformly over the space [0, 5]2 with a mesh 0.25 on each coordinate, with output values from (21). The
trained random forest model consists of 3 trees with 563, 535, 565 nodes respectively. The response
surface of the model g(x) is shown in Figure 9. Here we use γ = 500 and the shape of the rare-event
set {x : g(x) ≥ γ} in presented in Figure 10. We use Algorithm 2 to obtain the dominating point of theHuang, Lam, and Zhao

Figure 9: Response surface of the random Figure 10: Rare-event set {x : g(x) ≥ γ} offorest.
the random forest.

Figure 11: Probability estimation with dif- Figure 12: 95% confidence interval halfferent numbers of samples.
width with different numbers of samples.
rare-event set. In the algorithm, we choose the tolerance parameters to be ε = 0.01. This gives us thedominating point for the set at (3.00001, 2.62501).
We use the proposed IS and the uniform IS to estimate the probability P(g(x) ≥ γ). In Figure 11, we
observe that the estimates from the uniform IS (black dash line) are relatively unstable, compared to thosefrom the proposed IS (red solid line). In Figure 12, the confidence intervals of the uniform IS (black dashline) appear wider than the proposed IS. These comparisons are similar to the neural network case.
4.4 Random Forest Example: IS with Multiple Dominating Points
We now generate input samples uniformly over [0, 5]2 with a mesh 0.2 on each coordinate, with the output
values drawn from (22). The random forest has 2 trees in this example. The trained model using thegenerated dataset has 865 nodes in the first tree and 835 nodes in the second tree. Figure 13 shows the
response surface of the random forest. We consider γ = 8 and the shape of the rare-event set is shown in
Figure 14, which shows that the rare event set consists of two separate regions.
We find dominating points by implementing Algorithm 1 combined with Algorithm 2. Again, we use
ε = 0.01 and δ = 0.01 for the tolerance level. Note that the tolerance level we choose could generate anerror on the dominating point we obtain. For ε = 0.01, the error is given by kak2 − kâk2 ≤ 0.02, wherea is the true dominating point and â is the dominating point we obtain from the bisection. This error
might bring an issue to Algorithm 1. That is, the constraint â(x − â) < 0 we add may not cut off thecorresponding half-space. This can possibly lead to an infinite loop for Algorithm 1. Here we resolve thisissue
√ by modifying the additional constraint as â(x − â) < −δ , where δ > 0 is approximately in the scaleof ε. In this case, we use δ = 0.1. We obtain two dominating points, (0.1, 3.9) and (4.1, 3.3).

Huang, Lam, and Zhao

Figure 13: Response surface of the random Figure 14: Rare-event set {x : g(x) ≥ γ} offorest.
the random forest.

Figure 15: Probability estimation with dif- Figure 16: 95% confidence interval halfferent numbers of samples.
width with different numbers of samples.
The obtained dominating points allow us to implement the proposed IS. Compared to the uniform IS,
Figures 15 and 16 show that the proposed approach is more efficient as in the previous experiments. The
estimates of the proposed approach are more stable (red solid line in Figure 15) and the confidence intervals
are shorter (red solid line in Figure 16). The uniform IS (black dash line in Figure 16) requires roughly 3
times more samples to achieve the same level of accuracy as the proposed IS.
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.
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.

Huang, Lam, and Zhao
Bucklew, J. 2013. Introduction to Rare Event Simulation. New York: Springer Science & Business Media.
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.
Glasserman, P. 2013. Monte Carlo Methods in Financial Engineering, Volume 53. New York: Springer
Science & Business Media.
Goodfellow, I., Y. Bengio, A. Courville, and Y. Bengio. 2016. Deep Learning, Volume 1. Massachusetts:
MIT press Cambridge.
Huang, Z., H. Lam, D. J. LeBlanc, and D. Zhao. 2017. “Accelerated Evaluation of Automated Vehicles
Using Piecewise Mixture Models”. IEEE Transactions on Intelligent Transportation Systems:1–11.
Jain, A., H. S. Koppula, B. Raghavan, S. Soh, and A. Saxena. 2015. “Car That Knows Before You
Do: Anticipating Maneuvers via Learning Temporal Driving Models”. In Proceedings of the IEEE
International Conference on Computer Vision, 3182–3190. Piscataway, New Jersey: IEEE.
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.
Mišic, V. V. 2017. “Optimization of Tree Ensembles”. Working Paper: arXiv preprint arXiv:1705.10883.
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.
Tjeng, V., and R. Tedrake. 2017. “Verifying Neural Networks with Mixed Integer Programming”. Working
Paper: arXiv preprint arXiv:1711.07356.
Wang, W., D. Zhao, J. Xi, D. J. LeBlanc, and J. K. Hedrick. 2017. “Development and Evaluation of Two
Learning-Based Personalized Driver Models for Car-following Behaviors”. In Proceedings of American
Control Conference (ACC), 2017, 1133–1138. Piscataway, New Jersey: IEEE.
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.