Accelerated Evaluation of Automated Vehicles Using Piecewise Mixture Models

1.4 MB·

by  Zhiyuan Huang

IEEETRANSACTIONSONINTELLIGENTTRANSPORTATIONSYSTEMS,VOL.XX,NO.XX,DECEMBER2016 1
Accelerated Evaluation of Automated Vehicles
Using Piecewise Mixture Models
Zhiyuan Huang1 , Ding Zhao2, Henry Lam1, David J. LeBlanc2
Abstract—The process to certify highly Automated Vehicles both Google and Tesla update their software throughout the
hasnotyetbeendefinedbyanycountryintheworld.Currently, process, which may have improved safety, but the newest
companies test Automated Vehicles on public roads, which is
version of the AV has not accumulated that many miles as
time-consuming and inefficient. We proposed the Accelerated
they have claimed. In summary, today’s best practice adopted
7 Evaluation concept, which uses a modified statistics of the
1 surrounding vehicles and the Importance Sampling theory to by the industry is time-consuming and inefficient. A better
0 reduce the evaluation time by several orders of magnitude, approach is needed.
2 while ensuring the evaluation results are statistically accurate.
In this paper, we further improve the accelerated evaluation
n
concept by using PiecewiseMixture Distribution models, instead
a A. Related Researches
of Single Parametric Distribution models. We developed and
J
applied this idea to forward collision control system reacting BesidestheN-FOT,thetestmatrixapproach[6],[7]andthe
1
to vehicles making cut-in lane changes. The behavior of the
worst-casescenariosapproach[8],[9],[10]aretwoalternative
3
cut-in vehicles was modeled based on more than 403,581 lane
methods for vehicle evaluation. Our approach follows the
changes collected by the University of Michigan Safety Pilot
] ModelDeploymentProgram.Simulationresultsconfirmthatthe AcceleratedEvaluationconceptweproposed[11]toprovidea
Y
accuracy and efficiency of the Piecewise Mixture Distribution brand-new alternative. The basic concept is that as high-level
S method outperformed single parametric distribution methods in AVs just began to penetrate the market, they mainly interact
. accuracy and efficiency, and accelerated the evaluation process
s withhuman-controlledvehicles(HVs).Thereforewefocuson
c by almost four orders of magnitude.
modeling the interaction between the AV and the HV around
[
it. The evaluation procedure involves four steps:
1 Index Terms—automated vehicles, testing, evaluation, safety
v • Modelthebehaviorsoftheprimaryothervehicles(POVs)
5 represented by f(x) as the major disturbance to the AV
I. INTRODUCTION
1 using large-scale naturalistic driving data
9 IT is critical to thoroughly and rigorously test and evaluate • Skew the disturbance statistics from f(x) to modified
8 an Automated Vehicle (AV) before its release. Recent statistics f∗(x) (accelerated distribution) to generate
0 crashes involving a Google self-driving car [1] and a Tesla
more frequent and intense interactions between AVs and
.
1 Autopilotvehicle[2]attractedthepublic’sattentiontoAVtest- POVs
0 ing and evaluation. While these AVs are generally considered • Conduct accelerated tests with f∗(x)
7 as industrial leaders, because they use public road for testing,
1 • Use the Importance Sampling (IS) theory to skew back
statistically they have not yet accumulated enough miles. The
: the results to understand real-world behavior and safety
v TeslaAutopilot,inparticular,wascriticizedforbeingreleased
benefits
i
X too early in the hands of the general public [3].
This approach has been successfully applied to evaluate
Currently, there are no standards or protocols to test AVs
r AVs in the frontal crash with a cut-in vehicle [11] and also
a at automation level 2 or higher. Many companies adopt the
frontalcrashwithaleadvehicle[12],[13].Thisapproachwas
Naturalistic Field Operational Tests (N-FOT) approach [4].
confirmed to significantly reduce the evaluation time while
However, this method is inefficient because safety critical
accurately preserving the statistical behavior of the AV-HV
scenarios rarely happen in daily driving. The Google Self-
interaction. In the previous studies, the evaluation time was
driving cars accumulated 1.9 million driving. This distance,
reduced by two to five orders of magnitudes - the accelerated
although sounds a lot, provides limited exposure to critical
rate depends on the test scenarios, where rarer events achieve
events, given that U.S. drivers encounter a police reported
higher accelerated rate. The non-accelerated models and the
crash every five hundred thousand miles on average and fatal
accelerated models were built based on signal component
crash every one hundred million miles [5]. In the meantime,
distributions. While this method does benefit from its simple
*This work was funded by the Mobility Transformation Center at the mathematical form, it has a few drawbacks as illustrated in
University of Michigan with grant no. N021552. Z. Huang and D. Zhao Fig. 1. i) The fitting of the rare events (usually the tail part
contributedequallytotheresearch.
1Zhiyuan Huang ([email protected]) and Henry Lam of the statistical distributions) would be dominated by the
([email protected])areintheDepartmentofIndustrialandOperations fitting of the normal driving behaviors (the majority part of
EngineeringattheUniversityofMichigan. the distributions), which may induce large errors. ii) The full
2Ding Zhao (corresponding author: [email protected]) and
potential in higher accelerated rate is not achieved due to the
DavidLeBlanc([email protected])areintheUniversityofMichi-
ganTransportationResearchInstitute. lack of flexibility of the modified accelerated models.
IEEETRANSACTIONSONINTELLIGENTTRANSPORTATIONSYSTEMS,VOL.XX,NO.XX,DECEMBER2016 2
Fig. 1. Acceleration evaluation based on single parametric distribution and
PiecewiseMixtureDistribution.
Fig.2. LanechangedatacollectedbySPMDvehicle.
B. Contribution
In this paper, we proposed a more general framework for independently with regard to the value of v . By comparing
L
the Accelerated Evaluation method to overcome the afore-
among 17 types of commonly used distribution templates, we
mentioned limitations based on Piecewise Mixture Distribu- selected the Pareto distribution to model R−1 and used the
L
tion Models as illustrated in Fig. 1 b). In this paper, we exponential distribution for TTC−1 segments.
L
implemented the Accelerated Evaluation method under the
Using the empirical distribution of v and parametric dis-
L
new framework. Comparing to our previous work [14], we
tributions of R and TTC , we drew values from these
L L
thoroughly discuss the Cross Entropy method with proposed
distributions as inputs to simulate the AV-HV interaction.
framework in this paper. We present practical tips to over-
The outcome from the simulation can be considered as an
come numerical issues and reduce computational efforts. We
event indicator function I (x) that returns {1,0} depending
ε
demonstratethismethodbyevaluatingthelongitudinalcontrol
on the event of interest. Given the stochastic distribution of
system reacting to vehicles making cut-in lane changes.
the variables and the event indicator function, we obtained
the optimal exponential distribution for Importance Sampling
C. Paper Structure by implementing the Cross Entropy method [16]. As we have
shown in Fig. 1 a), we used only single parametric distribu-
Section II will introduce the lane change model based on
tions.Inthenextsection,weintroduceournewapproachusing
single parametric distributions. In Section III, we present the
Piecewise Mixture Distributions.
new lane change model with Piecewise Mixture Distributions.
We establish the Accelerated Evaluation in Section IV and
discusstheCrossEntropymethodwithPiecewiseMixtureDis- III. LANECHANGEMODELWITHPIECEWISEMIXTURE
tributionmodelsinSectionV.Simulationresultsarediscussed DISTRIBUTIONS
in Section VI. Section VII concludes this paper. Although many commonly used parametric distributions
have concise and elegant forms, they do not always describe
II. ACCELERATEDEVALUATIONWITHSINGLE the data distribution well. Instead, a better fitting can be
PARAMETRICDISTRIBUTIONS achieved by dividing the dataset into several subsets. We
estimatethemodelparametersusingtheMaximumLikelihood
The lane change events were extracted from the Safety
Estimation (MLE) [17] in each subset. The general process of
Pilot Model Deployment (SPMD) database [15]. With over
MLE is as follow.
2 million miles of vehicle driving data collected from 98
Assume we have a family of distribution with Cumula-
cars over 3 years, we identify 403,581 lane change events.
tive Distribution Function (CDF) F(x|θ), where θ is the
Previously [11], we used 173,692 events with a negative
parametervectorofF.ThecorrespondingProbabilityDensity
range rate to build a statistical model focusing on three key
Function (PDF) of F is f(x|θ). Assuming that data D =
variables that captured the effects of gap acceptance of the
{X ,X ,...,X } is independently and identically distributed
lanechangingvehicle:velocityoftheleadvehicle(v ),range 1 2 N
L and the distribution is in the family of F(x|θ), we want to
totheleadvehicle(RL)andtimetocollision(TTCL).TTCL find the most “likely” parameter θˆ.
was defined as:
R We define the likelihood function [18] as
TTC =− L, (1)
L R˙L L(θ|D)=P(D|θ)=ΠNn=1f(Xn|θ). (2)
where R˙L is the relative speed. We call the estimation of θˆ that maximizes the likelihood
The modeling of these three variables was hard to handle
function the mostly likely estimation MLE.
because of dependency, so we simplified it based on a crucial
For computation convenience, we introduce the log-
observation. Although TTC is dependent on v generally,
L L likelihood function
we split the data into 3 segments: v at 5 to 15 m/s, 15 to 25
L
N
m/sand25to35m/s.Withineachsegment,R isindependent (cid:88)
L L(θ|D)=lnL(θ|D)= lnf(X |θ). (3)
withv andTTC .ThisallowedustomodelTTC andR n
L L L L
n=1
IEEETRANSACTIONSONINTELLIGENTTRANSPORTATIONSYSTEMS,VOL.XX,NO.XX,DECEMBER2016 3
Since the logarithm is monotone, the log-likelihood function 1) MLEforboundedexponentialdistribution: Thebounded
preservestheoptimizeroftheoriginalfunction.[19]Theopti- exponential distribution with rate θ has the form
mizeroflog-likelihoodfunction,θˆ,istheMLEofdistribution
θe−θx
family F. We have the MLE as f(x|γ ≤x<γ )= (10)
1 2 e−θγ1 −e−θγ2
θˆ=argmax L(θ|D). (4)
for γ ≤x<γ .
θ 1 2
In the following, we describe the Piecewise Mixture Dis- For dataset D ={X1,...,XN}, the log-likelihood function
tribution fitting concept based on MLE and we present the is
bounded distribution fitting results. All optimization problems N
(cid:88)
presented in this section are tractable and can be solved by L(D|θ)= lnθ−θX −ln(e−θγ1 −e−θγ2), (11)
n
fminunc in MATLAB. n=1
where L is concave over θ. Although we cannot solve the
A. General Framework of the Piecewise Mixture Distribution maximization analytically, it is solvable through numerical
Lane Change Model methods.
We define Piecewise Mixture Distribution to be distribution Therefore, the MLE of θ is given by the optimization
with CDF in the form of
N
(cid:88)
(cid:88)k max Nlnθ−Nln(e−θγ1 −e−θγ2)− θXn. (12)
F(x)= π F (x|γ ≤x<γ ), (5) θ
i θi i−1 i n=1
i=1
2) MLE for bounded normal distribution: Consider a
where k is the number of truncation, (cid:80)ki=1πi = 1, and bounded normal distribution with mean 0 and variance θ2
Fi(x|γi−1 ≤x<γi)istheconditionalcumulativedistribution conditional on 0≤γ1 ≤x<γ2. The PDF is
function, meaning that F (γ |γ ≤ x < γ ) = 0 and
i i−1 i−1 i 1φ(x)
Fi(γi|γi−1 ≤ x < γi) = 1 for i = 1,...,k. θi denotes the f(x|γ ≤x<γ )= θ θ . (13)
parameter(s) for F . We can consider that π = P(γ ≤ 1 2 Φ(γ2)−Φ(γ1)
i i i−1 θ θ
x < γ ) and when x ≥ 0, we have γ = 0 and γ = ∞. By
i 0 k The MLE of the bounded normal distribution is given by
this definition, the PDF of the Piecewise Mixture Distribution
is
k (cid:80)N X2 γ γ
f(x)=(cid:88)π f (x|γ ≤x<γ ). (6) max− n=1 n −Nlnθ−Nln(Φ( 2)−Φ( 1)). (14)
i θi i−1 i θ 2θ2 θ θ
i=1
3) Fitting mixture model with EM algorithm: Compared to
In our case, θ = {π ,...,π ,θ ,...,θ }. Splitting D into
1 k 1 k single parametric distributions, mixture distribution combines
pieces regarding the truncation points {γ ,...,γ }, gives
1 k−1 several classes of distribution and thus is more flexible. We
data index sets S = {j|γ ≤ X < γ } for i = 1,...,k.
i i−1 j i consider the fitting problem of mixture bounded normal dis-
We can write the log-likelihood function as
tribution.
L(θ|D)=(cid:80)k (cid:80) lnπ The PDF of mixture of m bounded normal distribution can
+(cid:80)ki=1(cid:80) n∈Silnf i(X |γ ≤x<γ ). (7) be written as
i=1 n∈Si θi n i−1 i
m
We obtain the MLE of θ can be obtained by maximizing (cid:88)
f(x|γ ≤x<γ )= p f (x|γ ≤x<γ ) (15)
L(θ|D) over θ. Since L is concave over π , we take 1 2 j j 1 2
i
j=1
∂L
=0 (8) where f is bounded Gaussian distribution with mean
∂π j
i 0 and variance σ2. The parameters here are θ =
j
and get {p ,...,p ,σ2,...,σ2 }. We want to find MLE of p and σ2
πˆ =|S |/N. (9) 1 m 1 m j j
i i for j =1,...,m.
For parameters θi in Fi, it is known (7) to be the same The log-likelihood function for data D ={Xn}Nn=1 is
as computing the MLE of θ with corresponding dataset
i N m
Di = {X|γi−1 ≤ X < γi and X ∈ D}. Since we L(θ|D)= (cid:88)ln(cid:88)pjfj(Xn|γ1 ≤x<γ2). (16)
use bounded distribution for each F , below we explain the
i n=1 j=1
estimationofparametersforthethreedistributionsweapplied
We note that this is hard to solve directly, because there
in later sections.
is a sum within the log function. Therefore, we apply the
TosamplefromaPiecewiseMixtureDistribution,wecould
Expectation-Maximization (EM) [20] algorithm to find the
use the inverse function approach. See Appendix A for the
optimizer, i.e. MLE, for the parameters.
details.
We define Zj to denote whether or not the random number
n
X comesfrommixturedistributionj,j =1,...,m,andZj =
B. Bounded Distribution n n
{0,1}. We also introduce the expectation
Wedevelopthreeboundeddistributionsandusetheminthe
lane change model. E[Zj|X ]:=τj. (17)
n n n
IEEETRANSACTIONSONINTELLIGENTTRANSPORTATIONSYSTEMS,VOL.XX,NO.XX,DECEMBER2016 4
The EM algorithm starts with initial parameters {p ,σ }, with smaller variance. This is known as Importance Sampling
j j
j = 1,...,m. For data D = {X }N , we set com- [21] and F∗ is the IS distribution.
n n=1
plete data as D = {X ,Z }N . The EM algorithm op- For estimating P(X ∈ ε),we note that an optimal IS
c n n n=1
timizes E[L(θ|D )|D] in every step. The E step updates distribution
c
E[L(θ|D )|D], and the M step optimizes this function. The P(X ≤x, ε)
c F∗∗(x)=F(x|ε)= (25)
algorithm iterates E step and M step until reaching the P(x∈ε)
convergence criterion.
couldreducethevarianceofISestimationto0,buttheoptimal
In our case,
requires the knowledge of P(X ∈ ε). However, it guides the
E[L(θ|D )|D]=(cid:80)N (cid:80)m τj(lnp +lnf (X )). selection of the IS distribution.
c n=1 j=1 n j j n
(18)
B. Exponential Change of Measure
SinceobjectiveE[l (θ|D )|D]intheMstepisconcaveover
c c
p and σ , we could maximize the objective function through Exponential change of measure is commonly used to con-
j j
an analytic approach for p : structF∗.Althoughtheexponentialchangeofmeasurecannot
j
guarantee convergence to optimal distribution, it is easy to
(cid:80)N τj
p = n=1 n. (19) implement and the new distribution generally stays within the
j N same class of distribution.
For σ , we can solve the following maximization problem Exponential change of measure distribution takes the form
j
through numerical approach. of
f (x)=exp(θx−κ(θ))f(x), (26)
(cid:18) (cid:19) θ
X
σj =argmσijn−τnjlnσj +τnjlnφ σjn − where θ is the change of measure parameter and κ(θ) is
(cid:18) (cid:19) the log-moment generating function of original distribution f.
γ γ
τjln Φ( 2)−Φ( 1) . (20) When θ =0, we have f (x)=f(x).
n σ σ θ
j j For a bounded exponential distribution, the exponential
See Appendix B for the full EM algorithm. change of measure distribution is
(λ−θ)e−(λ−θ)x
IV. ACCELERATEDEVALUATIONWITHIMPORTANCE fθ(x|γ1 ≤x<γ2)= e−(λ−θ)γ1 −e−(λ−θ)γ2, (27)
SAMPLING
whereλistheparameterforexponentialdistribution.Wenote
Importance Sampling (IS) is thus used to accelerate the that f is still a bounded exponential distribution and λ =
θ θ
evaluationprocess,becausecrudeMonteCarlosimulationsfor λ−θ.
rare events can be time-consuming. Below we describe the IS For a bounded normal distribution, the exponential change
method. of measure distribution is
1φ(x−σ2θ)
A. Important Sampling and Optimal IS distribution fθ(x|γ1 ≤x<γ2)= Φ(γ2−θσσ2)−σΦ(γ1−θσ2), (28)
Let x be a random variable generated from distribution F, σ σ
where the original distribution truncated from a normal distri-
and ε ⊂ Ω where ε is the rare event of interest and Ω is the
bution with parameters µ=0 and σ. We note that the change
sample space. Our objective is to estimate
of measure distribution is still a bounded normal distribution
(cid:90)
P(X ∈ε)=E[I (X)]= I (x)dF (21) with µ=θσ2 and σ.
ε ε
V. CROSSENTROPYMETHODANDIMPLEMENTATION
where
(cid:40)
1 x∈ε, Section IV discussed optimal IS distribution F∗∗ providing
I (x)= (22)
ε 0 variance estimation to the value of interest, whereas this
0 otherwise.
section describes the Cross Entropy method used to estimate
We can write the evaluation of rare events as the sample the “optimal” parameters θ, which minimizes the “distance”
mean of Iε(x) between a parametric distribution Fθ and F∗∗ without know-
ing F∗∗. The description below is based on the Piecewise
N
Pˆ(X ∈ε)= 1 (cid:88)I (X ), (23) Mixture Distribution structure.
N ε n
n=1
A. Introduction
where X ’s are drawn from distribution F.
i
Since we have TheCrossEntropy,whichisalsoknownasKullback-Leibler
(cid:90) (cid:90) dF distance [22], measures the similarity between distributions.
E[Iε(X)]= Iε(x)dF = Iε(x)dF∗dF∗, (24) We define the Cross Entropy between function g and h as
g(X) (cid:90)
dwiestrcibauntiocnomFpu∗t,ewthheichsamhapslethmeeasnamoef sIuεp(Xpo)rtddFFw∗ithoveFr,thtoe D(g,h)=Eg[lnf(X)]= g(x)lng(x)dx−
(cid:90)
obtain an unbiased estimation of P(X ∈ε). By appropriately
g(x)lnh(x)dx. (29)
selecting F∗, the evaluation procedure obtains an estimation
IEEETRANSACTIONSONINTELLIGENTTRANSPORTATIONSYSTEMS,VOL.XX,NO.XX,DECEMBER2016 5
From (25), we know that the PDF of the optimal IS We note that if we have two independent variables where
distribution F∗∗ is f(x,y) = f(x)f(y), we can take a parametric distribution
f∗∗(x)= Iε(x)f(x). (30) fΘor=ea{cθh,vθar}i.abTlheeaonbdjehcativveeffuΘn(cxt,ioyn)c=orrfeθs1p(oxn)dfiθn2g(yt)o,(w33h)eries
P(x∈ε) 1 2
Since P(x ∈ ε) is generally unavailable, we use a para- N
mWeetrwicadnitsttroibfiuntidonthFeθptaoraamppertoearcθh∗ththeaotpmtiminaiml IiSzedsistthreibuCtrioosns. mθax N1 n(cid:88)=1Iε(Xn,Yn)ffΘ(sX(Xnn,Y,Ynn))(lnfθ1(Xn)+
Entropy [23] between f∗∗ and fθ. We denote θ∗ as the lnfθ2(Yn)), (34)
optimal parameter for the parametric distribution. Then the
which can be decoupled into two optimization problem over
minimization problem
θ and θ respectively and I (X ,Y ) f(Xn,Yn) is a known
minD(fθ,f∗∗) (31) co1nstant2given {Xn,Yn}. ε n n fΘs(Xn,Yn)
θ
We implement the Cross Entropy on the Piecewise Mixture
is equivalent to
Distribution with one variable. We note that we can apply the
f(X) results to the lane change model, since the Cross Entropy ob-
max E [I (X) lnf (X)], (32)
θ θs ε fθs(X) θ jective function of independent variables can be implemented
in (34).
where f denotes the sampling distribution with parameters
θs
θ . We note that this is a generalized setting, since we can
s
use any sampling distribution f as long as it has the same B. OptimizationFunctionforPiecewiseMixtureDistributions
θs
support with f. This is the baseline for iterations in the Cross We propose a parametric family of IS distribution for
Entropy method. We use the same form as fθ because in the Piecewise Mixture Distribution
following sections, we use a sampling distribution which is in
k
the same family as the parametric distribution. (cid:88)
f (x)= π˜ exp(θ x−κ(θ ))f (x|γ ≤x<γ ), (35)
Weestimateθ∗ bysolvingthestochasticcounterpartof(32) θ i i i i i−1 i
i=1
max 1 (cid:88)N I (X ) f(Xn) lnf (X ), (33) where we use exponential change of measure for each piece
θ N n=1 ε n fθs(Xi) θ n poafrdaimstertiebrutiisonθa=nd{θad,j.u.s.,tθth,eπ˜pr,o..p.o,πr˜tio}n. parameter to π˜i. The
1 k 1 k
wdihsterriebustiaomnpfleθss.{X1,...,XN} are drawn from the sampling datIan,(s3o3)w,ecnsi=mpIlεif(yXtnh)effθfsu(X(nXncnt)i)onisaasknownconstantgiventhe
We note that if I (X ) = 0 for all n = 1,..,N in (33),
ε n
N
the objective equals to 0 constantly. To avoid this situation, 1 (cid:88)
max c lnf (X ). (36)
we select a sampling distribution which emphasizes the rarer θ N n θ n
n=1
events.
Fig. 3 shows the iteration procedure of the Cross Entropy We split the samples into index sets Si = {j|γi−1 ≤
method. The core part of the Cross Entropy method is to use Xj < γi} for i = 1,...,k for each bounded segment. Since
theoptimizeroftheobjectivefunction(33)intheithiteration, fi(Xn|γi−1 ≤ x < γi) (cid:54)= 0 only if n ∈ Si, for each θi and
θi∗, as the parameters for the sampling distribution in the next π˜i, the optimization function is equivalent to
iteration. The underlying idea is that the IS distribution in
1 (cid:88)
distribution family f should better approach the optimal IS max c ln(π˜ exp(θ X −κ(θ ))
distribution.Thereforθe,asweiterate,weobtainmore“critical” θi,π˜i N n∈Si n i i n i
rareeventsandhaveabetterestimationoftheoptimizerwhich f (X |x<γ ≤x<γ )). (37)
i n i−1 i
leads to even more “critical” rare events in the next iteration.
We define the stopping criterion regarding the parameter or We can further rewrite the optimization function regarding
the objective value. In practice, we want to start with an θi and π˜i respectively. For π˜i, we have
appropriate sampling distribution to get a good solution with 1 (cid:88)
max c lnπ˜, (38)
laesssamiteprlaintigond.isSteriebuseticotnio.n V-C1 for a discussion of initializing π˜i N n∈Si n i
which obtains an analytical form for the optimizer
(cid:80)
c 1{n∈S }
π˜ = n∈Si n i . (39)
i (cid:80) c
n∈Si n
For θ , we have
i
1 (cid:88)
max c lnexp(θ X −κ(θ ))
θi N n∈Si n i n i
f (X |γ ≤x<γ ), (40)
Fig.3. IterationsofCrossEntropy. i n i−1 i
IEEETRANSACTIONSONINTELLIGENTTRANSPORTATIONSYSTEMS,VOL.XX,NO.XX,DECEMBER2016 6
which is an exponential change of measure with D only. 2) Adjusting sample size N: The choice of sample size N
i
We note that we can simplify this optimization function by should not only depend on the total number of rare events
rewriting the log term as obtained in each iteration. For each parameter of interest, we
need sufficient non-zero c ’s to guarantee the qualification of
1 (cid:88) n
max c (lnexp(θ X −κ(θ ))+ theestimation.Wenotethattheparametersestimationdepend
θi N n∈Si n i n i onlyontherareeventinthecorrespondingpiece,soweadjust
lnfi(Xn|γi−1 ≤x<γi)), (41) sample size N to ensure that each piece with large portion π˜i
contains enough rare event samples.
which is equivalent to
3) Setting a lower bound for π˜ : When we update π˜ in
i i
1 (cid:88)
mθaixN n∈Sicn(θiXn−κ(θi)), (42) e(3v9e)n,tisfacmnp=le 0infothrealplience∈, wSie,hmaevaeniπ˜nig=tha0t.tWhehreeniswneohraavree
π˜ = 0, the support of the IS distribution will differ from
since the latter term does not depend on θ . i
i
the original distribution. We note that it might cause bias in
For a bounded exponential distribution with parameter λ,
our simulation analysis. On the other hand, once π˜ hits 0, it
the Cross Entropy iteration solves i
will be 0 in the following iterations. Therefore, we need to
1 (cid:88) keep π˜ > 0. Setting a low bound for π˜ , for example, 0.01,
max c (θ X − i i
θi N n∈Si n i n when there is no rare event for piece i, gives an efficient IS
distribution while avoiding the problems.
e−(λ−θi)γi−1 −e−(λ−θi)γi
ln ). (43) 4) Updating parameter θ : The absence of rare event
λ−θ i
i samples also leads to failures in updating θ . In this case,
i
For a bounded normal distribution with parameters µ = 0
we use either the value of θ in the last iteration, or we set
i
andσ,theoptimizationfunctionfortheCrossEntropyiteration
it to 0, i.e. reset the distribution as the real distribution. We
is
note that we can tolerant some inaccurate estimation if π˜ is
i
max (cid:88) c X θ −((cid:88) c )(σ2θi2+ small, since a small π˜i indicates that this piece might not be
θi n∈Si n n i n∈Si n 2 im5p)orCtahnatntgointhgetrraurnecaevtieonntsγ. : The truncations of the Piece-
i
Φ(γi−θiσ2)−Φ(γi−1−θiσ2) wise Mixture Distribution are fixed throughout the Cross En-
ln σ σ ). (44)
Φ(γi)−Φ(γi−1) tropymethod.Thus,ifthereisabadselectionoftruncationin
σ σ
ouroriginaldistributionmodel,theCrossEntropycannotgive
C. Discussion on Numerical Implementation an efficient IS distribution. The changing of truncation points
WehavepresentedtheoptimizationfunctionsforCrossEn- is hard to implement by optimization, so we use a heuristic
tropy iterations, but we cannot reliably apply these equations approach for adjusting the truncation points to emphasize the
in practice without considering some of the problematical tail part of the Piecewise IS distribution.
numerical details. In this section, we discuss methods to In any iteration, if the number of rare events is not enough
overcome these numerical issues. to properly update the parameters, we check π˜ of the current
i
1) Initializing Cross Entropy Iterations for Rare Events: sampling distribution. If the π˜ of the tail piece is the largest
k
Since rare events occur with small probability, using the possible value, we increase the value of the all truncation
originaldistributionassamplingdistributiontostarttheCross points except γ with a certain value. Shifting the truncation
0
Entropy iterations it becomes computationally burdensome gives more weight to the tail part. Then by sampling from
to sample a single rare event. One possible approach is to the adjusted distribution, we check if the number of events of
initialize with guess of sampling distribution. When we have interest is sufficient. We repeat these actions until we obtain
some rough knowledge about the optimal IS distribution, enough rare events in the iteration.
we can use the knowledge to construct a proper sampling We propose this heuristic approach, since the flexibility of
distribution. the Piecewise Mixture Distribution is not fully exploited if
Forcaseswherewehavelittleknowledgeabouttheoptimal we cannot change the truncation points. We note that finding
IS distribution, we construct adaptive events that gradually a more systematic procedure to locate the knots remains an
reduce the rarity. For rare events denoted by ε, we define open question.
the sequence of events to be ε ⊃ ε ⊃ ... ⊃ ε ⊃ ε,
1 2 n
where ε is not rare for our initializing sampling density. For
1 VI. SIMULATIONANALYSIS
each iteration t, we gradually reduce the rare event set ε and
t
A. Automated Vehicle Model
use ε to replace ε in the objective function. Since ε is a
t t
subsetofε ,theISdistributionforε alsoprovidesmore First,wepresentourPiecewiseMixtureModelsforR−1and
t−1 t−1
chances for samples from ε . We use the optimal solution TTC−1 andthencomparetheresultswiththesingleparamet-
t
in (t−1)th iteration θ∗ as the sampling parameter θ for ric distribution model used in [11]. For both approaches, we
t−1 t
the next iteration and choose ε to have a relatively larger divide the data of TTC−1 into three segments regarding the
t
probability to occur under f . Since ε gradually approaches rangeofv.Sincethethreesegmentsaresimilarindistribution,
θt t
ε as we iterate, eventually we obtain the optimal parameters we only show the results of the segment for v in the range of
for ε. 15 to 25 m/s.
IEEETRANSACTIONSONINTELLIGENTTRANSPORTATIONSYSTEMS,VOL.XX,NO.XX,DECEMBER2016 7
Fig.4. PiecewiseMixtureDistributionfittingforR−1. Fig.6. ComparisonoffittingforR−1.
1) Piecewise mixture models for R−1 and TTC−1: Fig. 4
shows the fitting of R−1 using two bounded exponential dis-
tributionsandthreeboundedexponentialdistributions.Adding
one more truncation point provides a better fitting to the body
part of distribution while having the same fitting of tail.
In Fig. 5, we truncated the data into two parts. For the
tail part, we use the exponential distribution. For the body
part, the mixture of two normal distributions gives a better
fit. The Piecewise Mixture Models enable us to use different
distributions for the body part and the tail part.
2) Comparison with single parametric distribution models:
Fig.6 and Fig.7 compare the new model and the previous Fig.7. ComparisonoffittingforTTC−1 givenvL between15and25m/s.
model. We note that Piecewise Mixture Models provide more
flexibility in data fitting.
lane change model determines whether a crash happens by
checking to see if the value of R, the range between two
B. Cross Entropy Results
vehicles,reaches0.Meanwhile,theTTC alsogoesto0when
Here,weusethelanechangemodeltoexemplifytheCross a crash happens. To construct events less rare than a crash,
Entropy method. For the three variables R,TTC,v, the dis- we relax the criterion for crash to be either R hits t > 0
R
tribution is f(R,TTC,v)=f(v)f(R)f(TTC|v) where f(v) or TTC hits t > 0. By changing these two thresholds,
TTC
is the empirical distribution. Since we have three conditional t and t as shown in Fig. 8, we construct the adaptive
R TTC
distributions of TTC regarding the value of v, we find the rare events sequence for the Cross Entropy iterations. We use
IS distributions independently for each case. We present the sample size N =1000 for each iteration.
results for v from 5 to 15 m/s.
Fig. 9 and 10 show the parameters present in each of the
Weassumethatwehavelessinformationabouttherelation
iterations. We observe that the parameters stabilize gradually.
between the distribution of variables and the rare events.
Fig. 11 shows how the distribution changes gradually from
Our objective is to construct adaptive rare events to help
the original distribution to the IS distribution. We note that
us approach the IS distribution. We recall that our original
the density moves toward the tail part as we iterate.
C. Simulation Results
In our simulation experiments, we set the convergence
criterion as the relative half-width of 100(1−α)% confidence
interval drops below β. In this case, we use α=0.2 and β =
0.2 to study the number of samples needed for convergence.
OurgoalistocomparetheefficiencyofthePiecewiseMixture
Distribution and single exponential distribution models.
Fig. 12 shows that both models give a similar estima-
tion as the number of experiments grows large, and that
the Piecewise Mixture Distribution model converges slightly
faster than the single parametric model. The circles show that
the relative half-width of the Piecewise Mixture Distribution
Fig.5. PiecewiseMixtureDistributionfittingforTTC−1givenvLbetween
15and25m/s. model reaches the target confidence value after 7800 samples,
IEEETRANSACTIONSONINTELLIGENTTRANSPORTATIONSYSTEMS,VOL.XX,NO.XX,DECEMBER2016 8
Fig.8. CrossEntropyiterationswithsequenceofeventswiththresholdsfor
crash.Weleaveiteration1blanktokeepthex-axisconsistentwithFig.9and
10.
whereas the single parametric model needs about 13800 sam- Fig.10. CrossEntropyiterationswithsequenceofeventsofTTC−1 forv
from5to15m/s.
ples.UsingthePiecewiseMixtureDistributionmodelreduced
the sample size by 44%.
To reduce stochastic uncertainty, we repeat the tests 10
times and calculate the average. It takes 7840 samples on
average to obtain a converged estimation using the Piecewise
Mixture Distribution model, whereas it takes 12320 samples
on average using the single accelerated distribution model to
converge. Table I compares the two models with the crude
MonteCarlomethod[24].Weestimatethenumberneededfor
convergence of crude Monte Carlo by using the fact that the
numberofeventsofinterestoccurringisBinomialdistributed.
We compute the standard deviation of the crude Monte Carlo
estimation Pˆ(x∈ε) by
(cid:115) Fig.11. DistributionchangethroughCrossEntropyiterationswithsequence
Pˆ(x∈ε)(1−Pˆ(x∈ε)) ofeventsofTTC−1 forv from5to15m/s.
std(Pˆ(x∈ε))= , (45)
n
which allows us to estimate where zα/2 is the (1−α/2) quantile of normal distribution.
WecalculatetherequiredsamplesizeN ofcrudeMonteCarlo
z2 (1−Pˆ(x∈ε)) in Table I from an estimation Pˆ(x ∈ ε) = 7.4×10−7 with
Nˆ = α/2 , (46)
β2Pˆ(x∈ε) 80% confidence interval (7.0×10−7,7.8×10−7).
Finally, we apply the heuristic approach in Section V-C5 to
thedatasegmentwithv from5to15m/s.Werunsimulations
with this segment and compare the results with the standard
Fig.9. CrossEntropyiterationswithsequenceofeventsofR−1 forvfrom Fig.12. Estimationofcrashprobabilityforonelanechangeusingpiecewise
5to15m/s. andsingleaccelerateddistributions.
IEEETRANSACTIONSONINTELLIGENTTRANSPORTATIONSYSTEMS,VOL.XX,NO.XX,DECEMBER2016 9
TABLEI TABLEII
NUMBEROFSAMPLES(N)NEEDEDTOCONVERGE. COMPARISONOFTHECOMPUTATIONTIMEBETWEENSINGLE
PARAMETRICMODELANDPIECEWISEMODEL.
Piecewise Single Crude
N 7840 12320 5.5×107 Stages Crude Single Piecewise
RatiotoPiecewise 1 1.57 7×103 Fitting - 4parametersto 18parametersto
estimate estimate
30,000simulations 24,000simulations
CrossEntropy -
4parameters 18parameters
5.5×107 12,320 7840
Simulation
simulations simulations simulations
the same confidence level compared to the single parametric
Model.
APPENDIXA
INVERSECDFOFPIECEWISEMIXTUREDISTRIBUTIONS
We can sample from Piecewise Mixture Distribution by the
inverse CDF approach. Here, we derive the inverse CDF for
Piecewise Mixture Distribution.
Fig. 13. Relative half-width of crash probability estimation for one lane
changewithleadingvehicle’sspeedinrangeof5to15m/s,comparingsingle, The CDF of Piecewise Mixture Distribution (5) can split
piecewiseandheuristicaccelerateddistributions. into

...

approach for the ewise Mixture Distribution and single para- F(x)= (cid:80)i−1π +π F (x|γ ≤x<γ ) γ ≤x<γ .
j=1 j i i i−1 i i−1 i
metric distribution models. Fig. 13 shows the convergence of ...
confidence half-width. We note that the relative half-width of (47)
the heuristic, which is smaller than the standard approach for Therefore the inverse function can be written as
the Piecewise Mixture Distribution model, indicates that the 
...
latter model’s performance can be further improved. F−1(y)=F−1(y−(cid:80)ij−=11πj|γ ≤x<γ) (cid:80)i−1π ≤y<(cid:80)i π .
...i πi i−1 i j=1 j j=1 j
VII. CONCLUSIONS (48)
Thispaperproposedanewmodelforacceleratedevaluation where Fi−1 is the inverse conditional CDF of Fi. Below, we
of AVs. The Piecewise Mixture Distribution Models provide give two example of inverse conditional CDF.
more accurate fitting to the surrounding human-controlled FortheinverseCDFofconditionalexponentialdistribution,
vehicle behaviors than the single parametric model used we have
in the literature. The proposed model was more efficient
F−1(y|F (γ )≤y <F (γ ))=
and reduced the evaluation time by almost half than single θ θ 1 θ 2
parametric model. The Cross Entropy procedure described in Fθ−1((Fθ(γ2)−Fθ(γ1))y+Fθ(γ1)), (49)
this paper effectively worked in this scenario analysis. We
whereF andF−1aretheCDFandinverseCDFofexponential
provided practical solutions to deal with the numerical issues
distribution.
which occurred while calculating the optimal parameters. The
For conditional normal distribution, the inverse CDF is
heuristic approach exploited the flexibility of the Piecewise
Mixture Distribution structure. Testing the proposed model F−1(y|F (γ )≤y <F (γ ))=
θ θ 1 θ 2
on a large dataset of cut-in crashes caused by improper lane
γ −θσ2 γ −θσ2
changes, the Piecewise Mixture Distribution model reduced σΦ−1((Φ( 2 )−Φ( 1 ))y+
σ σ
the simulation cases by about 33% compared with the single
γ −θσ2
parametric model under the same convergence requirement. Φ( 1 ))+θσ2. (50)
σ
Moreover, the proposed model was 7000 times faster than the
Crude Monte Carlo method.
Table II summarizes the comparison of the computation
efforts between the models. We note that using the Piecewise
APPENDIXB
Mixture Distribution model increases the number of parame-
EMALGORITHMFORMIXTUREBOUNDEDNORMAL
ters estimated, where the estimation of parameters is almost
DISTRIBUTION
instant.Inthe CrossEntropystage,the numberofsimulations
required for the Piecewise model is not significantly less than Here, we present a numerical MLE algorithm with mixture
thesingleparametricmodel,becauseweassumenoknowledge bounded normal distribution. The steps are as follows.
about the optimal IS distribution for the Piecewise model. ALGORITHM:
Overall,thePiecewisemodelneedsfewersimulationstoreach 1) Initialize {p ,σ }, j =1,...,m.
j j
IEEETRANSACTIONSONINTELLIGENTTRANSPORTATIONSYSTEMS,VOL.XX,NO.XX,DECEMBER2016 10
2) E step: update [16] R. Y. Rubinstein, “Rare Event Simulation via Cross-entropy and Im-
portance Sampling,” in Second International Workshop on Rare Event
p f (X |σ )
τj = j j n j . (51) Simulation,1999,pp.1–17.
n (cid:80)m p f (X |σ ) [17] J. Aldrich, “RA Fisher and the Making of Maximum
j=1 j j n j Likelihood 1912-1922,” Statistical Science, 1997. [Online]. Available:
3) M step: update http://projecteuclid.org/euclid.ss/1030037906
(cid:80)N τj [18] D.CoxandD.Hinkley,TheoreticalStatistics,1979.
p = n=1 n (52) [19] S.BoydandL.Vandenberghe,ConvexOptimization,2004.
j N [20] A. Dempster, N. Laird, and D. Rubin, “Maximum Likelihood from
IncompleteDataviatheEMAlgorithm,”Journaloftheroyalstatistical
and
society,1977.
σ (Φ(γ2)−Φ(γ1)) [21] J. Bucklew, Introduction to Rare Event Simulation. Springer Science
σ =argmin−τjln j σj σj . (53) &BusinessMedia,2004.
j σj n φ(Xσjn) [[2223]] VD..VPa.pKnirko,esSet,atRis.ticYa.lRLuebairnnsitneginT,haenodryP,.1W99.8G.lynn, The Cross-Entropy
4) Repeat 2 and 3 until L(θ|D) converges. MethodforEstimation. ElsevierB.V.,2013,vol.31.
[24] S. Asmussen and P. Glynn, Stochastic Simulation: Algorithms and
Analysis. Springer,2007.
APPENDIXC [25] P.-T. d. Boer, “A Tutorial on the Cross-Entropy Method,” pp. 19–67,
2005.
VANILLACROSSENTROPYMETHOD
ALGORITHM:[25]
1) Initialize θ .
s
2) Sample {X ,...,X } from f and update
1 N θs
θ =argmax 1 (cid:88)N I (X ) f(Xi) lnf (X ). (54) Zprhei-ycuanadnidHatueanPgh.DZ.hisytuuadnenHtuinangIndisusatrisaelcoannddyOepa-r
θ N i=1 ε i fθs(Xi) θ i erationsEngineeringattheUniversityofMichigan,
AnnArbor.Hisresearchinterestsincludesimulation
3) Update andstochasticoptimization.
θ =θ. (55)
s
4) Repeat 2 and 3 until θ “converges”.
REFERENCES
[1] GoogleAutoLLC,“Monthlyreports?GoogleSelf-DrivingCarProject.”
[Online].Available:https://www.google.com/selfdrivingcar/reports/
[2] “PreliminaryReport,?HighwayHWY16FH018.”
[3] A.Evan,“FatalTeslaSelf-DrivingCarCrashRemindsUsThatRobots
Aren’tPerfect,”IEEESpectrum,2016. DingZhaoDingZhaoreceivedthePh.D.degreein
[4] FESTA-Consortium,“FESTAHandbookVersion2DeliverableT6.4of 2016 from theUniversity of Michigan, AnnArbor.
theFieldopErationalteStsupporTAction,”FESTA,Tech.Rep.,2008. HeiscurrentlyaResearchFellowintheUniversity
[5] NHTSA,“TrafficSafetyFacts2014,”DOT,Tech.Rep.,2014.[Online]. of Michigan Transportation Research Institute. His
Available:http://www-nrd.nhtsa.dot.gov/Pubs/811620.pdf\nhttp://www- research interest includes evaluation of connected
nrd.nhtsa.dot.gov/Pubs/809778.pdf and automated vehicles, vehicle dynamic control,
[6] H.PengandD.Leblanc,“EvaluationofthePerformanceandSafetyof driverbehaviorsmodeling,andbigdataanalysis.
AutomatedVehicles,”WhitePap.NSFTransp.CPSWork,2012.
[7] M.Aust,“EvaluationProcessforActiveSafetyFunctions:Addressing
KeyChallengesinFunctional,FormativeEvaluationofAdvancedDriver
AssistanceSystems,”2012.
[8] Y. Kou, “Development and Evaluation of Integrated Chassis Control
Systems,”2010.
[9] A.UngorenandH.Peng,“AnAdaptiveLateralPreviewDriverModel,”
VehicleSystemDynamics,vol.43,no.4,pp.245–259,42005.
[10] W. Ma and H. Peng, “A Worst-case Evaluation Method for Dynamic
Systems,”Journalofdynamicsystems,,1999. HenryLamHenryLamreceivedtheB.S.degreein
[11] D.Zhao,H.Lam,H.Peng,S.Bao,D.J.LeBlanc,K.Nobukawa,and actuarialsciencefromtheUniversityofHongKong
C. S. Pan, “Accelerated Evaluation of Automated Vehicles Safety in in2005,andtheA.M.andPh.D.degreesinstatistics
Lane-Change Scenarios Based on Importance Sampling Techniques,” from Harvard University, Cambridge, in 2006 and
IEEETransactionsonIntelligentTransportationSystems,2016. 2011.
[12] D.Zhao,H.Peng,S.Bao,K.Nobukawa,D.J.LeBlanc,andC.S.Pan, From2011to2014,hewasanAssistantProfessor
“Accelerated evaluation of automated vehicles using extracted natural- in the Department of Mathematics and Statistics
istic driving data,” in Proceeding for 24th International Symposium of at Boston University. Since 2015, he has been an
VehiclesonRoadandTracks,2015. Assistant Professor in the Department of Indus-
[13] D.Zhao,X.Huang,H.Peng,H.Lam,andD.J.Leblanc,“Accelerated trial and Operations Engineering at the University
EvaluationofAutomatedVehiclesinCar-FollowingManeuvers,”ArXiv, of Michigan, Ann Arbor. His research focuses on
p.12,2016. stochastic simulation, risk analysis, and simulation optimization. Dr. Lam’s
[14] Z.Huang,D.Zhao,H.Lam,D.J.LeBlanc,andH.Peng,“Evaluation works have been funded by National Science Foundation and National
of Automated Vehicles in the Frontal Cut-in Scenario - an Enhanced Security Agency. He has also received an Honorable Mention Prize in the
Approach using Piecewise Mixture Models,” 10 2016. [Online]. Institute for Operations Research and Management Sciences (INFORMS)
Available:http://arxiv.org/abs/1610.09450 GeorgeNicholsonBestStudentPaperAward,andFinalistinINFORMSJunior
[15] D. Bezzina and J. R. Sayer, “Safety Pilot: Model Deployment Test FacultyInterestGroupBestPaperCompetition.
Conductor Team Report,” NHTSA, Tech. Rep. June, 2014. [Online].
Available:http://safetypilot.umtri.umich.edu/