11institutetext: Instituto de Astrofísica de Canarias (IAC), Calle Vía Láctea
s/n, E-38200 La Laguna, Tenerife, Spain
11email<EMAIL_ADDRESS>22institutetext: Departamento de Astrofísica, Universidad
de La Laguna (ULL), E-38206 La Laguna, Tenerife, Spain 33institutetext:
Consejo Superior de Investigaciones Científicas (CSIC), E-28006 Madrid, Spain
44institutetext: Departamento de Química, Universidad de La Laguna (ULL),
E-38206 La Laguna, Tenerife, Spain 55institutetext: Center for Astrophysics
and Space Science, University of California San Diego, 9500 Gilman Drive, La
Jolla, CA 92092, USA

# Ammonia-methane ratios from H-band near-infrared spectra of late-T and Y
dwarfs

E. L. Martín 112233 J.-Y. Zhang 11 P. Esparza 44 F. Gracia 11 J. L. Rasilla 11
T. Masseron 1122 A. J. Burgasser 55

(,)

###### Abstract

Context.

Aims. Our goals are to investigate the relative absorption strengths of
ammonia and methane using low-resolution H-band (1.5-1.7 microns) spectra
obtained in the laboratory and compared with observational spectra of late-T
and Y dwarfs, and to estimate what can be expected from the wide-angle low-
resolution near-infrared spectroscopic survey that will be provided by the
upcoming Euclid space mission.

Methods. Gas cells containing ammonia and methane at atmospheric pressure were
custom-made in our chemical laboratory. Low-resolution near-infrared
spectroscopy of these gas cells were collected in our optical laboratory. They
are compared with simulated spectra using the high-resolution transmission
molecular absorption database (HITRAN) for temperatures of 300 K and 500 K,
and with near-infrared spectra of late-T and Y dwarfs, Jupiter and Saturn. We
selected for this investigation the spectral region between 1.5 and 1.7
microns (H-band) because it is covered by the Euclid red grism, it is
particularly sensitive to the relative proportions of ammonia and methane
opacity, and it is free from strong contributions of other abundant molecules
such as water vapour.

Results. The laboratory spectra showed that the ammonia and methane features
that are present in the simulations using the HITRAN database are incomplete.
Using our laboratory spectra, we propose a modified version of the NH3-H
spectral ratio with expanded integration limits that increases the amplitude
of variation of the index with respect to spectral type. Combinations of our
laboratory spectra were used to find the best fits to the observed spectra
with the relative absorption ratio of ammonia to methane as a free parameter.
A relationship was found between the Teff and the ratio of ammonia to methane
from spectral classes T5 to Y2 (1100 K – 350 K), in qualitative fairly good
agreement with theoretical predictions for high gravity objects and
temperatures from 1100 K to 500 K. The ammonia to methane ratios in late-T and
Y dwarfs are similar to that of Jupiter, suggesting a similar chemical
composition. Simulations of the spectroscopic performance of Euclid suggest
that it will yield Teff values and ratios of ammonia to methane for over 103
late T dwarfs in the whole wide survey.

###### Key Words.:

brown dwarfs – chemistry – exoplanets – optical – near-infrared – spectroscopy

##  1 Introduction

Brown dwarfs (BDs) do not stabilize on the Hydrogen (H) burning main-sequence,
and hence they keep cooling down to very low temperatures. The coolest BDs
have been spectroscopically classified as T dwarfs and Y dwarfs, based on the
presence of methane and ammonia in near-infrared spectra, respectively
(Burgasser et al., 2002; Delorme et al., 2008; Cushing et al., 2011), although
the presence of ammonia has been established not only in Y dwarfs but also in
late-T dwarfs (Bochanski et al., 2011; Canty et al., 2015).

The Wide Field Infrared Survey Explorer (WISE; Wright et al., 2010) has been
very effective to discover Y dwarfs and low-resolution spectroscopic follow-up
has been carried out from the ground and from space. In particular,
homogeneous sets of low-resolution near-infrared spectra have been collected
with Keck Near Infrared Spectrometer (NIRSPEC) and with the Hubble Space
Telescope (HST) for late-T and Y dwarfs down to spectral class Y2 (McLean et
al., 2003; Schneider et al., 2015; Cushing et al., 2021).

Estimates of the effective temperature (Teff) of Y dwarfs are difficult
because theoretical fits to their spectral energy distribution are of modest
quality, which indicates that the chemistry calculations and/or the sources of
opacity may be incomplete, and values estimated from their luminosity and
assumed radii using the Stefan-Boltzmann law are hampered by uncertainties due
to possible effects of unresolved binarity (Leggett et al., 2013; Luhman &
Esplin, 2016).

The appearance of ammonia in the near-infrared spectra of ultracool dwarfs has
been discussed by several authors (Saumon et al., 2000; Leggett et al., 2007;
Warren et al., 2007; Bochanski et al., 2011; Burgasser et al., 2012; Canty et
al., 2015). A spectral index named NH3-H was defined by Delorme et al. (2008)
to quantify the effects of ammonia absorption from 1.53 to 1.56 microns. Their
NH3-H index has been adopted for spectral type quantification from T5 to Y2
dwarfs (Cushing et al., 2021).

The upcoming Euclid111http://sci.esa.int/euclid/ space mission is expected to
provide the next order of magnitude leap in the numbers of detected ultracool
dwarfs. At the end of six years of nominal mission, the Euclid wide survey is
expected to reach 15,000 square degrees of the extragalactic sky observed at a
single epoch, and the Euclid deep surveys are expected to cover between 40 and
50 square degrees at multiple epochs. The Near Infrared Spectrometer and
Photometer (NISP) is a two-channel instrument for the Euclid space mission.
The photometric channel is equipped with three broad-band filters, and the
spectroscopic channel is equipped with three identical low-resolution grisms
covering the spectral range from 1250 nm to 1850 nm with a spectral resolving
power of R = 380, and one grism covering the spectral range from 920 nm to
1250 nm. The latter grism is expected to be used only for the Euclid Deep
Fields, and not for the wide survey. The spectroscopic observations will be
made in survey mode without any slit. The NISP field of view is 0.53 square
degrees covered with a matrix of 4x4 Teledyne detectors, each of them with
2040x2040 18 micron pixels with a scale of 0.3 arcsec per pixel.

Updated simulations of the number counts of ultracool dwarfs in the Euclid
wide survey using the most recent information on the photometric passbands
were provided in Solano et al. (2021). In this paper we will provide updated
simulations for the number of T and Y dwarfs expected to be detected with the
NISP spectroscopic mode in the Euclid wide survey. We will show that Euclid
will provide useful low-resolution H-band spectra for an unprecedented large
number of T dwarfs ($\geq 10^{3}$ for T5-T8 subclasses). Consequently, it is
important to investigate efficient methods to derive Teff values for T dwarfs
using NISP spectra. Our first goal is to confirm that the NH3-H spectral index
is really an indicator of the absorption of ammonia in late-T and Y dwarfs,
and to adjust its integration limits using our laboratory spectra. Our second
goal is the determination of the ammonia to methane ratio as a probe to the
atmospheres of ultracool dwarfs using low-resolution H-band spectra.

This paper presents laboratory spectra of gas cells filled with ammonia and
methane, which were made with the aim of helping with the identification of
these molecules in the atmospheres of ultracool dwarfs. The rest of the paper
is organized as follows: Section 2 describes the experimental setup for the
laboratory measurements. Section 3 compares the laboratory spectra with
simulated spectra using the HITRAN database and redefines the NH3-H index.
Section 4 compares the laboratory spectra with observed near-infrared spectra
from IRTF, Hubble, and Keck for dwarfs with spectral types from T5 to Y2, and
for Jupiter and Saturn. Ammonia to methane ratios are derived using a best
matching algorithm between the laboratory and the observed spectra. Section 5
summarizes the results and the implications.

##  2 Low-resolution laboratory spectra

Custom made gas cells were made to fit into the spectrophotometer Cary 5000
located at the optical laboratory of the Instituto de Astrofísica de Canarias
(IAC). They are G2 infrasil cells from International Crystal Lab, 13 cm long
and with window size 38 mm in diameter and 6 mm thick.

The cells were filled with different molecular gas using a vacuum pump system
in our laboratory located at the chemistry department of the Universidad de La
Laguna. The procedure for filling the cells with gas is described in more
detail in Valdivielso et al. (2010). The gas cells with ammonia (NH3), argon
(Ar), and methane (CH4) were completely filled with each of those gases at the
atmospheric pressure.

The gas cells were taken to the IAC and they were left to settle to the
ambient conditions inside the Cary 5000 for at least five hours before taking
the measurements. The Cary 5000 has two channels that allow for two gas cells
to be measured simultaneously. In one channel we placed a gas cell filled with
ammonia or methane, and in the other channel we placed an identical gas cell
filled with Ar so that the instrumental response would be corrected by the
software of the instrument.

The laboratory spectra of all the gas cells were taken on 8/8/2021 with an
integration time of 0.2 s for each step of 0.5 nm. The total time to measure a
spectrum from 800 nm to 3000 nm was 6.3 hours. The ambient conditions in the
laboratory were temperature 298 K and humidity 57 % .

Figure 1:  Our laboratory spectrum for ammonia and methane compared to
simulated spectra of the same molecules using HITRAN for temperatures of 500 K
and 300 K. The wavelength regions proposed for the new NH3-H index are marked
in color bands: green for the numerator, red for the denominator.

##  3 Comparison between laboratory and simulated spectra

Our laboratory spectrum for ammonia and methane were compared with simulated
spectra for the same molecules using the HITRAN database (Gordon et al., 2017)
for temperatures of 300 K and 500 K. The comparisons in the full spectral
range from 1.0 to 2.5 microns are shown in Fig. 1. The agreement is quite good
between 2.0 and 2.5 microns, but the differences increase towards shorter
wavelengths. This is not a surprise because it has been reported before that
HITRAN is very incomplete for hot bands and high-J molecular transitions
(Bailey & Kedziora-Chudczer, 2012; Wong et al., 2019).

We use our laboratory spectrum of ammonia to propose a redefinition of the
integration limits for the spectral index NH3-H, which was originally defined
using HITRAN opacity data (Delorme et al., 2008). These authors used as
numerator the spectral region from 1.53 to 1.56 microns, and as denominator
the region from 1.57 to 1.60 microns. Based on the appearance of the ammonia
band in our laboratory spectra, we propose to redefine the NH3-H index using
as numerator the region from 1.50 to 1.57 microns, and as denominator the
region from 1.57 to 1.61 microns. These integration limits are shown in the
bottom part of Fig. 1. They sample a region dominated by the ammonia
absorption and not by methane.

The increased wavelength range chosen for our proposed modification of the
NH3-H spectral index is motivated by a better understanding of the ammonia
absorption in this region and it has the advantage of increasing the flux
collected in each passband for the classification of ultra faint objects of T
and Y type detectable with NISP. The comparison between previous NH3-H index
and ours of dwarfs from spectral type T5 to Y2 are shown in Fig. 2.

Figure 2: The redefined NH3-H index using our laboratory spectra has a higher
sensitivity to spectral class (larger amplitude of variation) than the
original definition by Delorme et al. (2008) based on HITRAN data. Third order
polynomial fits to the relationship between the indeces and the spectral types
are denoted with dashed lines of different colors as labelled.

##  4 Comparison between laboratory and HST spectra

Different combinations of the relative strengths of ammonia and methane were
computed using our laboratory spectra. An example of the typical combination
that provides the best fit is shown in Fig. 3 together with HST spectra of a
T8 dwarf and two Y dwarfs (Schneider et al., 2016; Cushing et al., 2014). This
combination, denoted with a blue line, corresponds to the case where we scaled
our laboratory spectrum of ammonia by a factor of 0.33 and we added it to our
laboratory spectrum of methane.

The spectral region from 1.5 to 1.7 microns is particularly sensitive to the
effects of changing the relative contribution of ammonia with respect to
methane. Furthermore this spectral region is thought to be free from
significant contribution of other molecules such as water vapour (Delorme et
al., 2008; Canty et al., 2015).

Least square fit and a confidence interval estimate were computed using the
`lsqcurvefit` and `nlparci` routines available in the Optimization Toolbox and
Statistics and Machine Learning Toolbox in the MATrix LABoratory (MATLAB). We
found that the best fit to the HST and NIRSPEC spectra of late-T and Y dwarfs
in the spectral window from 1.5 to 1.7 microns were achieved with the ratios
of ammonia to methane shown in Fig. 4.

Using the same method as described above, we also derived the ammonia to
methane ratio in the spectra of Jupiter and Saturn that we downloaded from the
NASA Infrared Telescope Facility (IRTF) archive (Vacca et al., 2003). Our
ratio for Jupiter is plotted in Fig. 4 and it is consistent (within the 1
sigma uncertainty) with that found in the literature using models and
observations of the Galileo spacecraft (0.26, Taylor et al. (2007)). Our ratio
for Saturn agrees with Owen et al. (1977) who reported on non-detection of
ammonia in high-resolution spectra of Saturn obtained at Palomar observatory
in the wavelength region relevant for this work.

In Fig. 4 (lower panel) we show the comparison of our ammonia to methane
ratios with the theoretical predictions for high-gravity objects
($g=10^{5}\mathrm{cm\cdot s^{-2}}$) from (Zahnle & Marley, 2014). To compare
the behaviour of the models with other data, we used the ratios of theoretical
column densities scaled by a constant factor to make two values reaching an
agreement at 600K. There is a good qualitative agreement in the sense that the
models predict a ratio that has a minimum around 800 K and increases towards
cooler and hotter temperatures.

Figure 3:  Our laboratory spectrum for the typical best fitting combination of
ammonia and methane, and pure methane, compared to HST spectra of a T8, two Y
dwarfs and two giant gas planets.

Figure 4:  Ammonia to methane ratios derived for late T dwarfs and Y dwarfs
with 95% confidence level, from NIRSPEC data (blue) and HST data (red)
compared with Jupiter and Saturn from IRTF data (red) as functions of spectral
type (upper panel) and Teff (lower panel). The Teff of T, Y dwarfs, Jupiter
and Saturn adopted in this work are listed in Table 1. The pink dashed curves
are the 3rd-degree polynomial fittings for all object (except Saturn) against
spectral type and Teff and the black dashed line is the linear fitting later
than T8. The results are compared with the scaled theoretical ratios estimated
from Zahnle & Marley (2014), which are shown in light blue color.

##  5 Discussion and final remarks

Our laboratory spectra of ammonia and methane contain absorption features that
are too weak or not present in the HITRAN database. We used our laboratory
spectrum of ammonia to redefine the spectral index NH3-H in late-T and Y
dwarfs. This index is a sensitive indicator of the spectral type of ultracool
brown dwarfs and it can be used to calibrate the data provided by the slitless
spectroscopic mode of the NISP instrument on Euclid. Our new defined NH3-H
index has broader integration limits and an amplitude of variation with
spectral type almost twice larger than the one based on HITRAN and originally
proposed by Delorme et al. (2008) from spectral type T5 to Y2.

There is a smooth relation between the ammonia to methane ratio derived from
our best fits to the observed spectral region from 1.5 to 1.7 microns and Teff
of the objects considered, which is qualitatively in agreement with
theoretical predictions taken from the ratio of the theoretical column
densities of ammonia and methane under a high surface gravity environment
($g=10^{5}\mathrm{cm\cdot s^{-2}}$) (Zahnle & Marley, 2014). All the late T
and Y dwarfs considered in this work have ammonia to methane ratios consistent
with that of Jupiter within the uncertainties, suggesting that they have
similar chemical abundances.

Using the same assumptions as in Solano et al. (2021) and the densities of
objects of different spectral types in the solar neighborhood Kirkpatrick et
al. (2021), we provide here estimates of the numbers of dwarfs with spectral
types from T0 to Y2 that could be detected by NISP spectroscopy with SNR$>$10
in the Euclid wide survey (Fig. 5). Note that the numbers of T5 to T8 dwarfs
spectroscopically characterized by NISP are expected to be $\geq$103, and they
rapidly decrease into the Y dwarf domain due to the intrinsic faintness at
near-infrared wavelengths of those sources. The results of this work indicate
that ammonia to methane ratios can be measured in late-T and Y dwarfs using
low-resolution near-infrared spectra such as those that will be provided by
Euclid/NISP. These ratios measured in a large sample with NISP will be useful
to explore the diversity of chemical and physical properties that are likely
to exist among the coolest objects in the solar vicinity.

Figure 5:  Simulated number counts of T and Y dwarfs detectable (SNR$\geq$10) in the full Euclid wide survey with the spectroscopic mode of NISP for different values of the galactic scale height of the brown dwarf population in the think disk of the Galaxy.  Table 1:  Adopted Teff values and uncertainties for the different spectral types and giant gas planets considered in this work.  Type | Teff | References
---|---|---
| (K) |
T5 | 1100$\pm$100 | Dupuy & Liu (2017)
T6 | 1000$\pm$100 | Dupuy & Liu (2017)
T8 | 600$\pm$50 | Schneider et al. (2015)
T9 | 550$\pm$50 | Schneider et al. (2015)
Y0 | 450$\pm$50 | Schneider et al. (2015)
Y1 | 400$\pm$50 | Schneider et al. (2015)
Y2 | 300$\pm$50 | Cushing et al. (2021)
Jupiter | 134$\pm$4 | Aumann et al. (1969)
Saturn | 97$\pm$4 | Aumann et al. (1969)

###### Acknowledgements.

ELM acknowledges support from the Agencia Estatal de Investigación del
Ministerio de Ciencia e Innovación (AEI-MCINN) under grant
PID2019-109522GB-C53. J.-Y. Zhang acknowledges a summer grant from the
Instituto de Astrofisica de Canarias. ELM thanks Mike Cushing for sending the
HST infrared spectrum published in Cushing et al. (2021). This research has
made use of the Simbad database, operated at the centre de Données
Astronomiques de Strasbourg (CDS), and of NASA’s Astrophysics Data System
Bibliographic Services (ADS).

## References

  * Aumann et al. (1969) Aumann, H., Gillespie Jr, C., & Low, F. 1969, The Astrophysical Journal, 157, L69
  * Bailey & Kedziora-Chudczer (2012) Bailey, J. & Kedziora-Chudczer, L. 2012, Monthly Notices of the Royal Astronomical Society, 419, 1913
  * Bochanski et al. (2011) Bochanski, J. J., Burgasser, A. J., Simcoe, R. A., & West, A. A. 2011, AJ, 142, 169
  * Burgasser et al. (2012) Burgasser, A. J., Gelino, C. R., Cushing, M. C., & Kirkpatrick, J. D. 2012, ApJ, 745, 26
  * Burgasser et al. (2002) Burgasser, A. J., Kirkpatrick, J. D., Brown, M. E., et al. 2002, ApJ, 564, 421
  * Canty et al. (2015) Canty, J. I., Lucas, P. W., Yurchenko, S. N., et al. 2015, MNRAS, 450, 454
  * Cushing et al. (2011) Cushing, M. C., Kirkpatrick, J. D., Gelino, C. R., et al. 2011, ApJ, 743, 50
  * Cushing et al. (2014) Cushing, M. C., Kirkpatrick, J. D., Gelino, C. R., et al. 2014, AJ, 147, 113
  * Cushing et al. (2021) Cushing, M. C., Schneider, A. C., Kirkpatrick, J. D., et al. 2021, arXiv e-prints, arXiv:2107.00506
  * Delorme et al. (2008) Delorme, P., Delfosse, X., Albert, L., et al. 2008, A&A, 482, 961
  * Dupuy & Liu (2017) Dupuy, T. J. & Liu, M. C. 2017, ApJS, 231, 15
  * Gordon et al. (2017) Gordon, I. E., Rothman, L. S., Hill, C., et al. 2017, Journal of Quantitative Spectroscopy and Radiative Transfer, 203, 3
  * Kirkpatrick et al. (2021) Kirkpatrick, J. D., Marocco, F., Caselden, D., et al. 2021, ApJL, 915, L6
  * Leggett et al. (2013) Leggett, S. K., Morley, C. V., Marley, M. S., et al. 2013, ApJ, 763, 130
  * Leggett et al. (2007) Leggett, S. K., Saumon, D., Marley, M. S., et al. 2007, ApJ, 655, 1079
  * Luhman & Esplin (2016) Luhman, K. L. & Esplin, T. L. 2016, AJ [arXiv:1605.06655]
  * McLean et al. (2003) McLean, I. S., McGovern, M. R., Burgasser, A. J., et al. 2003, ApJ, 596, 561
  * Owen et al. (1977) Owen, T., McKellar, A., Encrenaz, T., et al. 1977, Astronomy and Astrophysics, 54, 291
  * Saumon et al. (2000) Saumon, D., Geballe, T. R., Leggett, S. K., et al. 2000, ApJ, 541, 374
  * Schneider et al. (2015) Schneider, A. C., Cushing, M. C., Kirkpatrick, J. D., et al. 2015, ApJ, 804, 92
  * Schneider et al. (2016) Schneider, A. C., Greco, J., Cushing, M. C., et al. 2016, ApJ, 817, 112
  * Solano et al. (2021) Solano, E., Gálvez-Ortiz, M. C., Martín, E. L., et al. 2021, MNRAS, 501, 281
  * Taylor et al. (2007) Taylor, F. W., Atreya, S. K., Encrenaz, T., et al. 2007, The Composition of the Atmosphere of Jupiter, 59
  * Vacca et al. (2003) Vacca, W. D., Cushing, M. C., & Rayner, J. T. 2003, PASP, 115, 389
  * Valdivielso et al. (2010) Valdivielso, L., Esparza, P., Martín, E. L., Maukonen, D., & Peale, R. E. 2010, ApJ, 715, 1366
  * Warren et al. (2007) Warren, S. J., Mortlock, D. J., Leggett, S. K., et al. 2007, MNRAS, 375, 213
  * Wong et al. (2019) Wong, A., Bernath, P. F., Rey, M., Nikitin, A. V., & Tyuterev, V. G. 2019, The Astrophysical Journal Supplement Series, 240, 4
  * Wright et al. (2010) Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868
  * Zahnle & Marley (2014) Zahnle, K. J. & Marley, M. S. 2014, ApJ, 797, 41


reaction coordinates of dynamical systems, journal = Applied and Computational
Harmonic Analysis, volume = 21, number = 1, pages = 113-127, year = 2006, note
= Special Issue: Diffusion Maps and Wavelets, issn = 1063-5203, doi =
10.1016/j.acha.2005.07.004, url = https://doi.org/10.1016/j.acha.2005.07.004,
author = B. Nadler and S. Lafon and R. R. Coifman and I. G. Kevrekidis,
abstract = A central problem in data analysis is the low dimensional
representation of high dimensional data and the concise description of its
underlying geometry and density. In the analysis of large scale simulations of
complex dynamical systems, where the notion of time evolution comes into play,
important problems are the identification of slow variables and dynamically
meaningful reaction coordinates that capture the long time evolution of the
system. In this paper we provide a unifying view of these apparently different
tasks, by considering a family of diffusion maps, defined as the embedding of
complex (high dimensional) data onto a low dimensional Euclidean space, via
the eigenvectors of suitably defined random walks defined on the given
datasets. Assuming that the data is randomly sampled from an underlying
general probability distribution p(x)=e−U(x), we show that as the number of
samples goes to infinity, the eigenvectors of each diffusion map converge to
the eigenfunctions of a corresponding differential operator defined on the
support of the probability distribution. Different normalizations of the
Markov chain on the graph lead to different limiting differential operators.
Specifically, the normalized graph Laplacian leads to a backward Fokker–Planck
operator with an underlying potential of 2U(x), best suited for spectral
clustering. A different anisotropic normalization of the random walk leads to
the backward Fokker–Planck operator with the potential U(x), best suited for
the analysis of the long time asymptotics of high dimensional stochastic
systems governed by a stochastic differential equation with the same potential
U(x). Finally, yet another normalization leads to the eigenfunctions of the
Laplace–Beltrami (heat) operator on the manifold in which the data resides,
best suited for the analysis of the geometry of the dataset regardless of its
possibly non-uniform density.

@ArticleKobak2021tsneglobal, author=Kobak, Dmitry and Linderman, George C.,
title=Initialization is critical for preserving global data structure in both
t-SNE and UMAP, journal=Nature Biotechnology, year=2021, month=Feb, day=01,
volume=39, number=2, pages=156-157, issn=1546-1696,
doi=10.1038/s41587-020-00809-z, url=https://doi.org/10.1038/s41587-020-00809-z

@ARTICLEsenanayake2020song, author=Senanayake, Damith A. and Wang, Wei and
Naik, Shalin H. and Halgamuge, Saman, journal=IEEE Transactions on Neural
Networks and Learning Systems, title=Self-Organizing Nebulous Growths for
Robust and Incremental Data Visualization, year=2021, volume=32, number=10,
pages=4588-4602, doi=10.1109/TNNLS.2020.3023941

@articlemcinnes2018umap-software, title=UMAP: Uniform Manifold Approximation
and Projection, author=McInnes, Leland and Healy, John and Saul, Nathaniel and
Grossberger, Lukas, journal=The Journal of Open Source Software, volume=3,
number=29, pages=861, year=2018

@miscnaim2019life-expt, author = Mohamed Amine Mairech, Naim1997, title =
DataC’EPT : Life expectancy prediction , publisher = Kaggle, year = 2019, url
= https://kaggle.com/competitions/datacept-life-expectancy-prediction

@inproceedings rhodes2023supervised, title=Supervised Manifold Learning via
Random Forest Geometry-Preserving Proximities, author=Jake Slater Rhodes,
booktitle=Fourteenth International Conference on Sampling Theory and
Applications, year=2023, url=https://openreview.net/forum?id=t6E4dZjp-e

@miscrosaler2023enhanced, title=Towards Enhanced Local Explainability of
Random Forests: a Proximity-Based Approach, author=Joshua Rosaler and Dhruv
Desai and Bhaskarjit Sarmah and Dimitrios Vamvourellis and Deran Onay and
Dhagash Mehta and Stefano Pasquali, year=2023, eprint=2310.12428,
archivePrefix=arXiv, primaryClass=stat.ML

@articleanand2011vne, title = Shannon and von Neumann entropy of random
networks with heterogeneous expected degree, author = Anand, Kartik and
Bianconi, Ginestra and Severini, Simone, journal = Phys. Rev. E, volume = 83,
issue = 3, pages = 036109, numpages = 8, year = 2011, month = Mar, publisher =
American Physical Society, doi = 10.1103/PhysRevE.83.036109, url =
https://link.aps.org/doi/10.1103/PhysRevE.83.036109

@InProceedingspierre2012autoencoders, title = Autoencoders, Unsupervised
Learning, and Deep Architectures, author = Baldi, Pierre, booktitle =
Proceedings of ICML Workshop on Unsupervised and Transfer Learning, pages =
37–49, year = 2012, editor = Guyon, Isabelle and Dror, Gideon and Lemaire,
Vincent and Taylor, Graham and Silver, Daniel, volume = 27, series =
Proceedings of Machine Learning Research, address = Bellevue, Washington, USA,
month = 02 Jul, publisher = PMLR, pdf =
http://proceedings.mlr.press/v27/baldi12a/baldi12a.pdf, url =
https://proceedings.mlr.press/v27/baldi12a.html, abstract = Autoencoders play
a fundamental role in unsupervised learning and in deep architectures for
transfer learning and other tasks. In spite of their fundamental role, only
linear autoencoders over the real numbers have been solved analytically. Here
we present a general mathematical framework for the study of both linear and
non-linear autoencoders. The framework allows one to derive an analytical
treatment for the most non-linear autoencoder, the Boolean autoencoder.
Learning in the Boolean autoencoder is equivalent to a clustering problem that
can be solved in polynomial time when the number of clusters is small and
becomes NP complete when the number of clusters is large. The framework sheds
light on the different kinds of autoencoders, their learning complexity, their
horizontal and vertical composability in deep architectures, their critical
points, and their fundamental connections to clustering, Hebbian learning, and
information theory.

@articletheis2017lossy, title=Lossy image compression with compressive
autoencoders, author=Theis, Lucas and Shi, Wenzhe and Cunningham, Andrew and
Huszár, Ferenc, journal=arXiv preprint arXiv:1703.00395, year=2017

@articleromero2017quantum, title=Quantum autoencoders for efficient
compression of quantum data, author=Romero, Jonathan and Olson, Jonathan P and
Aspuru-Guzik, Alan, journal=Quantum Science and Technology, volume=2,
number=4, pages=045001, year=2017, publisher=IOP Publishing

@misctitanic2016kaggle, author = Creator: Kaggle, title = Titanic: Machine
Learning from Disaster, howpublished = https://www.kaggle.com/c/titanic/data,
year = 2016, note = Accessed: 11/13/2023,

@articlelecun2010mnist, title=MNIST handwritten digit database, author=LeCun,
Yann and Cortes, Corinna and Burges, CJ, journal=AT&T Labs [Online].
Available: http://yann.lecun.com/exdb/mnist, year=2010

@article Rhodes2023RFPHATE, author = Jake S Rhodes and Adrien Aumon and Sacha
Morin and Marc Girard and Catherine Larochelle and Boaz Lahav and Elsa Brunet-
Ratnasingham and Wei Zhang and Adele Cutler and Anhong Zhou and Daniel E
Kaufmann and Stephanie Zandee and Alexandre Prat and Guy Wolf and Kevin R.
Moon, title = Gaining Biological Insights through Supervised Data
Visualization, elocation-id = 2023.11.22.568384, year = 2023, doi =
10.1101/2023.11.22.568384, publisher = Cold Spring Harbor Laboratory, abstract
= Dimensionality reduction-based data visualization is pivotal in
comprehending complex biological data. The most common methods, such as PHATE,
t-SNE, and UMAP, are unsupervised and therefore reflect the dominant structure
in the data, which may be independent of expert-provided labels. Here we
introduce a supervised data visualization method called RF-PHATE, which
integrates expert knowledge for further exploration of the data. RF-PHATE
leverages random forests to capture intricate feature-label relationships.
Extracting information from the forest, RF-PHATE generates low-dimensional
visualizations that highlight relevant data relationships while disregarding
extraneous features. This approach scales to large datasets and applies to
classification and regression. We illustrate RF-PHATEs prowess through three
case studies. In a multiple sclerosis study using longitudinal clinical and
imaging data, RF-PHATE unveils a sub-group of patients with non-benign
relapsing-remitting Multiple Sclerosis, demonstrating its aptitude for time-
series data. In the context of Raman spectral data, RF-PHATE effectively
showcases the impact of antioxidants on diesel exhaust-exposed lung cells,
highlighting its proficiency in noisy environments. Furthermore, RF-PHATE
aligns established geometric structures with COVID-19 patient outcomes,
enriching interpretability in a hierarchical manner.Competing Interest
StatementThe authors have declared no competing interest., URL =
https://www.biorxiv.org/content/early/2023/11/23/2023.11.22.568384, eprint =
https://www.biorxiv.org/content/early/2023/11/23/2023.11.22.568384.full.pdf,
journal = bioRxiv

@articlexiao2017fashion, title=Fashion-mnist: a novel image dataset for
benchmarking machine learning algorithms, author=Xiao, Han and Rasul, Kashif
and Vollgraf, Roland, journal=arXiv preprint arXiv:1708.07747, year=2017

@articlecoifman2006geometric, title=Geometric harmonics: a novel tool for
multiscale out-of-sample extension of empirical functions, author=Coifman,
Ronald R and Lafon, Stéphane, journal=Applied and Computational Harmonic
Analysis, volume=21, number=1, pages=31–52, year=2006, publisher=Elsevier

@articlebengio2003out, title=Out-of-sample extensions for lle, isomap, mds,
eigenmaps, and spectral clustering, author=Bengio, Yoshua and Paiement, Jean-
françcois and Vincent, Pascal and Delalleau, Olivier and Roux, Nicolas and
Ouimet, Marie, journal=Advances in neural information processing systems,
volume=16, year=2003


# Adaptive Splitting of Reusable Temporal Monitors for Rare Traffic Violations

Craig Innes and Subramanian Ramamoorthy1 Authors are with the School of
Informatics, University of Edinburgh, 10 Crichton St, EH8 9AB, United Kingdom.
Corresponding author<EMAIL_ADDRESS>
###### Abstract

Autonomous Vehicles (avs) are often tested in simulation to estimate the
probability they will violate safety specifications. Two common issues arise
when using existing techniques to produce this estimation: If violations occur
rarely, simple Monte-Carlo sampling techniques can fail to produce efficient
estimates; if simulation horizons are too long, importance sampling techniques
(which learn proposal distributions from past simulations) can fail to
converge. This paper addresses both issues by interleaving rare-event sampling
techniques with online specification monitoring algorithms. We use adaptive
multi-level splitting to decompose simulations into partial trajectories, then
calculate the distance of those partial trajectories to failure by leveraging
robustness metrics from Signal Temporal Logic (stl). By caching those partial
robustness metric values, we can efficiently re-use computations across
multiple sampling stages. Our experiments on an interstate lane-change
scenario show our method is viable for testing simulated av-pipelines,
efficiently estimating failure probabilities for stl specifications based on
real traffic rules. We produce better estimates than Monte-Carlo and
importance sampling in fewer simulations.

##  I Statistical Simulation for AV Testing

Autonomous Vehicles (avs) typically undergo rigorous simulated testing before
deployment [39]. A standard set of steps for testing is as follows: First we
define a scenario (e.g., a highway lane-change expressed in OpenScenario
[42]). Next, we define a safety specification (e.g., _“avoid impeding traffic
flow”_) in a formal language like Signal Temporal Logic (stl). Then, we run
stochastic simulations to estimate the probability our av-system violates our
specification [11]. This _statistical simulation_ approach is used because
modern avs contain “black box” components like Neural-Network perception
modules and non-linear solvers [45]. Such components provide few analytical
guarantees over their behaviour.

A core problem plaguing statistical simulation is estimating _rare events_.
Consider a stochastic simulation scenario where there exists a $10^{-4}$
probability that random noise in the sensors will cause our av to “fail”
(i.e., to violate our safety specification). If we ran $100$ simulations, it
is likely none would produce a failure. Even if sampling did produce a
failure, estimation variance would be unacceptable [24].

Many works address rare-event problems for avs via _Importance Sampling_ [25,
6]: Importance samplers draw simulations from a proposal distribution where
the factors leading to a failure occur more frequently. The final estimate of
failure probability is then re-weighted to reflect the original distribution.
Since we do not know in advance all combinations of states which result in
failure, such techniques must _learn_ a good proposal. This learning step has
no convergence guarantees, and probability estimates from such adaptive
techniques can have unbounded error [3]. Importance sampling also tends to
fare better when failures are caused by instantaneous single-state errors, but
in the av domain, failures often occur as a result of accumulated errors over
dependent states [9].

_“Preserve Traffic Flow”_

$\centering\color[rgb]{.25,.25,.25}\scriptsize\square_{[0,\infty]}(\neg\text{slow\\_leading\\_vehicle}(x^{\textit{ego}},x^{o_{1}\dots
o_{3}})\implies\text{preserves\\_flow}(x^{\textit{ego}}))\@add@centering$

Figure 1: Lane-change. Moving vehicles (blue) shown with trajectory. ‘Ego’
vehicle must avoid static obstacles (red). We monitor the safety constraint
shown in English and stl.

This paper instead proposes an approach to av rare-event simulation based on
merging _Adaptive Multi-level Splitting_ (ams) [9] with stl monitoring. ams
relies on estimating probabilities for a sequence of decreasing failure
thresholds $\gamma_{1}>\gamma_{2}\dots>\gamma_{M}$, where the final
$\gamma_{M}$ is equivalent to the rare failure event of interest. The key idea
is that estimating any intermediate $\gamma_{i}$ (given $\gamma_{i-1}$) is
easier than estimating $\gamma_{M}$ outright. To adapt ams from estimating
isolated phenomena (e.g., particle transport [30]) to estimating complex av-
system failures, we face two issues:

The main issue is how to consistently produce simulations which fall below
those intermediate failure thresholds $\gamma_{0\dots M}$, and how to
efficiently measure the distance to failure in the first place. Our approach
measures failure using metrics for evaluating stl specification robustness. By
leveraging online monitoring [13], we can cache the robustness values of
partial trajectories, stop simulations at the point where they fall below the
current threshold, and re-sample from this point onwards to produce
trajectories which fall below subsequent failure thresholds.

A secondary issue is generating stochastic av perceptual errors. Approaches
which assume noise follows a well-known (e.g., Gaussian) state-independent
distribution [7] are insufficient to capture the perceptual variety of a
typical av-system—a LiDAR detector may be great for close range traffic, but
terrible for long range or occluded traffic [36]. We therefore use a
_Perception Error Model_ (pem) [37]—a surrogate trained on real sensor data
which mimics perceptual errors encountered in regular operation (Sec II-A).

The main contribution of this paper is a new method for assessing failure
probability in av-scenarios (Sec III), combining ams (Sec III-B), pems (Sec
II-A), and online stl monitoring (Sec III-C). Our experiments focus on the
case study of a highway lane-change scenario, and show our method can be used
to test a full av-pipeline—perception down to control (Sec IV). Our approach
outperforms Monte-Carlo sampling, as well as fixed and adaptive importance
sampling, across various stl specifications (Sec IV-A).

To limit the scope of experiments, this paper exclusively considers
_probabilistic noise in the perception system_ as the primary source of
simulation stochasticity (As is standard in other works [10]). However, our
proposed sampling method can easily be applied to simulators which consider
other sources of stochasticity such as those arising from traffic behaviour or
physical uncertainties.

##  II Probability of Failure in Black Box Simulation

Consider the lane-change maneuver in Fig (1). Our car (the left-most _ego
vehicle_), must change to the left lane to avoid an obstacle, then re-merge.
Formally, let’s assume our scenario takes place over a total of $T$ time
steps. We denote the $d$-dimensional state of our ego vehicle at time $t$ as
$x^{ego}_{t}\in\mathds{R}^{d}$; other vehicles as $x^{o_{i}}_{t}$. For
succinctness, we write
$x_{t}=\\{x^{ego}_{t},x^{o_{0}}_{t},\dots,x^{o_{M}}_{t}\\}$ for the combined
state. The state $x_{t}$ contains the position, velocity and rotation of each
vehicle. At each time step $t$, the ego vehicle’s control system takes an
action $a_{t}\in\mathds{R}^{2}$ (desired acceleration and turn-velocity) with
the aim of minimizing costs associated with competing driving goals (e.g.,
maintaining a reference velocity and minimizing abrupt steering), and subject
to constraints (e.g., limits on acceleration, avoiding collisions, staying
within road boundaries). We can run a _simulation_ of this system to generate
a trajectory $\tau=[(x_{0},a_{0})\dots(x_{T},a_{T})]$, where
$\tau_{[i{:}j]}=[(x_{i},a_{i})\dots,(x_{j},a_{j})]$) denotes a partial slice.
For a given scenario, we wish to test whether our above av-system will violate
an stl safety specification $\varphi$. Due to probabilistic noise in the
perception system, our simulator is inherently stochastic. Therefore our aim
is to calculate the probability that, for a random run of our simulator, our
av-system will violate $\varphi$:

$P_{\textit{fail}}=\mathds{E}\left[\mathds{1}\\{\tau\nvDash\varphi\\}\right]$
(1)

Where $\mathds{1}\\{\tau\nvDash\varphi\\}$ is an indicator function which
returns $1$ if $\tau$ violates $\varphi$ and $0$ otherwise. To explain how our
method efficiently calculates (1), we first cover the pre-requisites for
perception, tracking, and control for simulating our av (Sec II-A-II-B). We
then cover defining safety specifications $\varphi$ in stl, and how to
quantify their satisfaction using a robustness metric (III-A). We can then
describe our main contribution—interleaving online monitoring and Adaptive
Multi-level Splitting to estimate a failure probability for av-systems via
statistical simulation (III-B).

Figure 2: Lane change with pem observations, tracking, and prediction.
Green/Orange crosses show pem obstacle observations. Purple dots/lines mark
estimated/predicted positions.

###  II-A Simulated State Estimation with PEMs

Most av problems assume our system does not have access to the true state
$x_{t}$. Instead our av must estimate this state via observations from its
sensors (see Fig (2)). At time $t$, let us denote a full-snapshot of the world
by $w_{t}$. This snapshot implicitly contains the relevant ground-truth state
$x_{t}$, but also other scene information (e.g., vehicle types, dimensions
etc.). A sensor $\mathcal{S}$ can be thought of as a function which takes
$w_{t}$ and produces high-dimensional raw sensor data $\mathcal{S}(w_{t})$
(e.g., a LiDAR point-cloud). This sensor data is then passed to a perception
function $f$ (e.g., a neural-network obstacle detector [38]), to produce an
_observation_ $y_{t}\in\mathbb{R}^{d^{\prime}}$ (e.g., the bounding boxes of
other vehicles relative to the ego):

$y_{t}=f(\mathcal{S}(w_{t}))$ (2)

By using standard tracking algorithms [48], we can use these observations to
get an estimate $\hat{x}_{t}$ of the current state:

$\hat{x}_{t}=\mathds{E}\left[x\mid y_{0\dots t}\right]$ (3)

If we have a one-step vehicle dynamics model $f_{\textit{dyn}}$, we can use it
to predict the state in future time steps:

$\displaystyle\hat{x}_{t+i\mid t}=\underbrace{f_{\textit{dyn}}\circ\dots\circ
f_{\textit{dyn}}}_{\text{$i$ times}}(\hat{x}_{t})$ (4)

Here, $\hat{x}_{t+i|t}$ denotes the predicted state at time $t$+$i$ given an
estimate at $i$. In a real-world system, this setup allows us to sense,
estimate, track and predict the state; in simulation, we have a problem: $f$
is typically a data-driven perception module trained on real sensors, but most
simulators _cannot generate_ high fidelity sensor inputs (e.g., photo-
realistic images). We can resolve this issue by re-framing our perception
system as a noisy projection from the state space to observation space:

$y_{t}=f(\mathcal{S}(w_{t}))=Hx_{t}+\epsilon(g(w_{t})))$ (5)

In this view, $f$ is composed of a $d^{\prime}\times d$ projection matrix $H$
on state $x_{t}$, plus stochastic error dependent on the current world state
$w_{t}$. The $\epsilon$ function is a surrogate model known as a _Perception
Error Model_ [37]. This is a probabilistic model of the original av’s
perception noise, dependent on _salient features_ $g(w)$ extractable from
simulated $w$. Salient features can include obstacle positions, dimensions,
occlusion, or environment factors. We model $\epsilon$ as a _gaussian process_
[49], where $m$, $\kappa$ are mean/kernel functions111The nuances of gp-
inference and kernel choice are beyond the scope of this paper, but see [17]
for discussion.:

$\displaystyle\epsilon(w)\sim\textsc{GP}(m(g(w)),\kappa(g(w),g(w^{\prime})))$
(6)

Now, instead of using real-world sensor inputs directly, our simulator applies
probabilistic noise from the pem based on the current simulated state. This
makes each run of the simulator inherently stochastic, as different amounts of
perceptual noise may be applied to observations on each run.

###  II-B Model Predictive Control for Highway Maneuvers

For a sufficiently complex control task such as lane changing, choosing the
best actions $a_{t}$ at each time step is a _non-linear constrained control_
task. We phrase the controller of our av system under test as a _Receding-
Horizon Model-Predictive Control_ optimization [41]. Eq (7) provides a formal
definition of the optimization problem: At each time step $t$, the controller
aims to choose actions $a_{t:(t+H)}$ over a finite time horizon $H$ which
minimize a cost function $J(x,a)$. Actions must be chosen subject to a set of
constraints $c_{j}(x,a)$ and obey physical dynamics $f_{dyn}$:

$\begin{array}[]{@{}ll}\min\limits_{\begin{subarray}{c}x_{t:(t+H)}\\\
a_{t:(t+H)}\end{subarray}}&\sum\limits_{k=t}^{t+H}J(x_{k},a_{k})\\\
s.t.&x_{t}=\hat{x}_{t}\\\ &\forall k,\
x_{k+1}=f_{\textit{dyn}}(x_{k},a_{k})\\\ &\forall j,\ c_{j}(x_{k},a_{k})<0\\\
\end{array}$ (7)

Cost function $J(x,a)$ balances multiple factors such as tracking a reference
velocity, minimizing abrupt movement, and staying close to the lane centre.
The constraint functions $c_{j}(x,a)$ ensure states and actions remain
feasible (e.g., that the car stays within road bounds and acceleration
limits).

We give a further breakdown of the cost function and implementation in Section
(IV). However, the purpose of describing Eq (7) here in the context of our
testing problem is to highlight that our av-controller represents yet another
“black-box” component of our system: Despite behaving deterministically,
solvers for nonlinear control problems are not guaranteed to find a global
solution, and can perform arbitrarily poorly [18]. Other typical methods (such
as reinforcement learning), pose the same problem.

##  III Estimating Specification Failure Probability

We have defined our testing goal and outlined the components of our system-
under-test. Now we can show how we formalize our safety properties, how we
draw samples from the simulator, and how those aspects interact.

###  III-A Specifying Safety with Signal Temporal Logic

We can express av traffic rules involving statements about continuous values
over time using Signal Temporal Logic [33]. stl has grammar:

$\displaystyle\varphi:=$
$\displaystyle\top\mid\eta\mid\neg\varphi\mid\varphi_{1}\wedge\varphi_{2}\mid\square_{I}\varphi\mid$
(8)
$\displaystyle\Diamond_{I}\varphi\mid\varphi_{1}\mathcal{U}_{I}\varphi_{2}\mid
H_{I}\varphi\mid O_{I}\varphi$

Here, $\eta$ is any predicate $\rho(s)-b\leq 0$ (with $b\in\mathbb{R}$ and
$\rho:\mathcal{X}\rightarrow\mathbb{R}$). $\square_{I}\varphi$ means $\varphi$
is _always_ true at every future time in interval $I$. $\Diamond_{I}\varphi$
means $\varphi$ must _eventually_ be true in $I$,
$\varphi_{1}\mathcal{U}_{I}\varphi_{2}$ means $\varphi_{1}$ must remain true
within $I$ _until_ $\varphi_{2}$ is true. $H_{I}$, $O_{I}$ are past versions
of $\square_{I}$, $\Diamond_{I}$.

We can convert $\tau$ and $\varphi$ to a _robustness metric_
$\mathcal{L}(\varphi,\tau,i)$ over trajectories [19]. This measures how
strongly $\tau$ satisfied $\varphi$ (starting from time $i$). Large positive
values indicate robust satisfaction; negative values a strong violation; near-
zero values a trajectory on the boundary of satisfaction/violation. Eq (9)
shows a subset of $\mathcal{L}$’s semantics. Other operators follow similar
definitions [13].

$\begin{array}[]{ll}\mathcal{L}(\rho(\tau)>0,\tau,i)&=\rho(x_{i})\\\
\mathcal{L}(\neg\varphi,\tau,i)&=-\mathcal{L}(\varphi,\tau,i)\\\
\mathcal{L}(\varphi_{1}\wedge\varphi_{2},\tau,i)&=\min(\mathcal{L}(\varphi_{1},\tau,i),\mathcal{L}(\varphi_{2},\tau,i))\\\
\mathcal{L}(\square_{I}\varphi,\tau,i)&=\inf\limits_{i^{\prime}\in
i+I}\mathcal{L}(\varphi,\tau,i^{\prime})\end{array}$ (9)

###  III-B Estimating the Rare Event with Splitting

With our pem, av-controller, specification $\varphi$ and metric $\mathcal{L}$,
we now describe our adaptive sampling contribution: Given $N$ simulated
trajectories $\tau^{(1\dots N)}$, we wish to calculate Eq (1)—the probability
a simulation violates $\varphi$. A naive approach might use _Monte-Carlo_
sampling:

$\mathds{E}\left[\mathds{1}\\{\mathcal{L}(\varphi,\tau,0)<\gamma\\}\right]\approx\frac{1}{N}\sum_{i=1}^{N}\mathds{1}\\{\mathcal{L}(\varphi,\tau_{i},0)<\gamma\\}$
(10)

However, when violating $\varphi$ is rare, the number of simulations needed to
achieve low relative error rapidly becomes infeasible [24]. We instead take a
_multi-level splitting_ approach [9]. Rather than immediately estimating
failure ($\gamma=0$) we instead estimate decreasing thresholds
$\gamma_{1}>\gamma_{2}>\dots>\gamma$. At each stage $m$, starting with $N$
trajectories we take two steps. First, we _discard_ trajectories that do not
fall below threshold $\gamma_{m}$. Second, we _replenish_ back up to $N$
trajectories. To replenish discarded trajectories, we clone a random un-
discarded $\tau^{(i)}$ up to $t^{\prime}$—the first time step where
$\mathcal{L}(\tau^{(i)}_{[0{:}t^{\prime}]},\varphi,0)<\gamma_{m}$. Then, we
re-simulate $\tau^{(i)}$ from $t^{\prime}$ to $T$. This ensures all $N$
trajectories at stage $m$ are below $\gamma_{m}$. With staged, partial re-
samplings, Eq (1) becomes:

$\prod\limits_{m=1}^{M}P(\mathcal{L}(\varphi,\tau,0)<\gamma_{m}|\mathcal{L}(\varphi,\tau,0)<\gamma_{m-1})$
(11)

Given enough levels, each conditional probability should be significantly
larger than $P(\mathcal{L}(\varphi,\tau,0)<\gamma)$. The final computation:

$\hat{p}_{ams}=\left\\{\prod^{M}_{m=1}\frac{N-K_{m}}{N}\right\\}\times\frac{1}{N}\sum^{N}_{i=1}\mathds{1}\\{\mathcal{L}(\varphi,\tau^{(i)},0)\\}$
(12)

has $M$ stages and $N$ initial simulations, where $K_{m}$ is the number
discards per stage. Unlike adaptive importance sampling, ams guarantees
convergence as $N\rightarrow\infty$ [8].

###  III-C Adaptive Splits via Online STL Monitors

To achieve our high-level goal of adapting ams to av testing of stl
specifications, we currently face a slight computational dilemma: Our
robustness metric $\mathcal{L}(\tau,\varphi,0)$ defines a single batch
robustness value from the start to the end of a complete trajectory, and takes
computation time proportional to the length of $\tau$. However in Section
(III-B), we saw that the discard and replenishment steps require access to the
robustness values at _arbitrary prefixes_ of trajectories. In other words, ams
re-simulation requires online computation of _all partial trajectories_ :

$\left\\{\mathcal{L}(\tau^{(i)}_{[0{:}t^{\prime}]},\varphi,0)\mid
t^{\prime}\in[0,T],i\in[0,N]\right\\}$ (13)

We resolve this dilemma by taking a key insight from the online monitoring
literature—we can interleave partial computations into the ams process. To
achieve this interleaving, we must first slightly alter our definition of
$\mathcal{L}$ from (9) to define our metric in terms of partial rather than
full trajectories. Eq (14) is a modified metric $\mathcal{L}_{n}$ (where $n$
references the “nominal semantics” of [13], augmented with past operators).
This definition makes explicit that we only partially evaluate trajectory
$\tau$ up to fixed time step $t$:

$\begin{array}[]{@{}l@{}l@{}}\mathcal{L}_{\textit{n}}(\rho(x)>0,\tau_{[0{:}t]},i)&=\
\rho(\tau_{[i]})\\\
\mathcal{L}_{\textit{n}}(\neg\varphi,\tau_{[0{:}t]},i)&=-\mathcal{L}_{\textit{n}}(\varphi,\tau_{[0{:}t]},i)\\\
\mathcal{L}_{\textit{n}}(\square_{I}\varphi,\tau_{[0{:}t]},i)&=\inf\limits_{i^{\prime}\in(i+I\cap[0,t])}(\mathcal{L}_{\textit{n}}(\varphi,\tau_{[0{:}t]},i^{\prime}))\\\
\mathcal{L}_{\textit{n}}(\Diamond_{I}\varphi,\tau_{[0{:}t]},i)&=\sup\limits_{i^{\prime}\in(i+I\cap[0,t])}(\mathcal{L}_{\textit{n}}(\varphi,\tau_{[0{:}t]},i^{\prime}))\\\
\mathcal{L}_{\textit{n}}(O_{I}\varphi,\tau_{[0{:}t]},i)&=\sup\limits_{i^{\prime}\in(i-I\cap[0,t])}(\mathcal{L}_{\textit{n}}(\varphi,\tau_{[0{:}t]},i^{\prime}))\\\
\mathcal{L}_{\textit{n}}(H_{I}\varphi,\tau_{[0{:}t]},i)&=\inf\limits_{i^{\prime}\in(i-I\cap[0,t])}(\mathcal{L}_{\textit{n}}(\varphi,\tau_{[0{:}t]},i^{\prime}))\\\
\mathcal{L}_{\textit{n}}(\varphi_{1}\mathcal{U}_{I}\varphi_{2},\tau_{[0{:}t]},i)&=\sup\limits_{i_{2}\in(i+I\cap[0,t])}\min\left(\mathcal{L}_{\textit{n}}(\varphi_{2},\tau_{[0{:}t]},i_{2}),\vphantom{\inf\limits_{i_{i}\in[i,i_{2}]}}\right.\\\
\end{array}$ (14)

With the above definition, we could now naively compute all members of (13) by
re-evaluating $\varphi$ at every $i$, up to every $t$, at every re-sampling
step. This is computationally wasteful, as the robustness values of many
partial trajectories share many operations with the computations of their
prefixes. To take advantage of this fact, we instead maintain a _work-list_
for each $\tau$. A _work-list_ stores a mapping from specification $\varphi$
and time step $t$ to robustness value
$\mathcal{L}_{n}(\varphi,\tau_{[0{:}t]},0)$. By using a work-list, we can
obtain a robustness value for the partial trajectories at time $t$+$1$ using
just the newly available state $x_{t+1}$ and the previous values of the work-
list, rather than repeating computations over the entire trajectory.

Alg (1) takes as input the current work-list at time step $t$ and sketches how
it is updated online using the newly arrived state $x_{t+1}$: For predicates
$\rho(x)$, incoming state $x_{t+1}$ is added only if it is in $\varphi$’s
_time horizon_. For example, for $\varphi=\square_{[0,2]}(\rho(x)>0)$, state
$x_{3}$ would not be added, as it falls outside the relevant interval. For
formulas like negation and conjunction, pointwise operators are leveraged to
combine the existing results from previous sub-formula computations. Similarly
for temporal operators, we can use the sliding min-max algorithm of [27] to
compute running maxes over the relevant sub-formula intervals. Further
optimizations can be added (e.g., replacing chunks of a work-list with
‘summaries’ as sufficient information arrives [13]) but we omit the details
here. We can access the robustness of a partial trajectory from the updated
work-list by querying $\textit{w-list}[\varphi][0]$.

1 Function  _upd-wl(_w-list, $\varphi$, $x_{t+1}$_)_

2 switch  _$\varphi$_ do

3 case  _$\rho(x) >0$_ do

4 if  _$t$ +$1$ is within time horizon of $\varphi$_ then
$\textit{w-list}[\varphi][t+1]\leftarrow\rho(x_{t+1})$ ;

5

6 case  _$\neg\psi$_ do

7 upd-wl(_w-list, $\psi$, $x_{t+1}$_)

8 $\textit{w-list}[\varphi]\leftarrow$ Pointwise negation of
$\textit{w-list}[\psi]$

9 case  _$\psi_{1}\wedge\psi_{2}$_ do

10 upd-wl(_w-list, $\psi_{1}$, $x_{t+1}$_)

11 upd-wl(_w-list, $\psi_{2}$, $x_{t+1}$_)

12 $\textit{w-list}[\varphi]\leftarrow$ Pointwise mins of
$\textit{w-list}[\psi_{1}]$ and $\textit{w-list}[\psi_{2}]$

13 case  _$\square_{I}\psi$_ do

14 upd-wl(_w-list, $\psi$, $x_{t+1}$_)

15 $\textit{w-list}[\varphi]\leftarrow\ $Sliding min window of width $|I|$
across $\textit{w-list}[\psi]$

16 case  _$\Diamond_{I}\psi$_ do

17 upd-wl(_w-list, $\psi$, $x_{t+1}$_)

18 $\textit{w-list}[\varphi]\leftarrow\ $ Sliding max window of width $|I|$
across $\textit{w-list}[\psi]$

19

20 case  _$\psi_{1}\mathcal{U}_{I}\psi_{2}$_ do

21 upd-wl(_w-list, $\psi_{1}$, $x_{t+1}$_)

22 upd-wl(_w-list, $\psi_{2}$, $x_{t+1}$_)

23 $\textit{lr-mins}\leftarrow$ Pointwise mins of $\textit{w-list}[\psi_{1}]$
and $\textit{w-list}[\psi_{2}]$

// Calculate backwards inductively

24

25 for  _$i$ in descending timesteps_ do

26 $\textit{us}[i]\leftarrow$$\max\left(\textit{lr-
min}[i],\min(\textit{w-list}[\psi_{1}][i],us[i+1])\right)$

27 $\textit{w-list}[\varphi]\leftarrow us$

28 case  _$H_{I}\psi$_ do

29 upd-wl(_w-list, $\psi$, $x_{t+1}$_)

30 $\textit{w-list}[\varphi]\leftarrow$ Sliding min window of width $|I|$
across (reversed) $\textit{w-list}[\psi]$

31 case  _$O_{I}\psi$_ do

32 upd-wl(_w-list, $\psi$, $x_{t+1}$_)

33 $\textit{w-list}[\varphi]\leftarrow$ Sliding max window of width $|I|$
across (reversed) $\textit{w-list}[\psi]$

34

Algorithm 1 Update Work List (Adapted from [13])

1

2 Function  _stl-ams(_$x_{0}$ , $\varphi$, $\gamma$, $T$, $K$, $N$_)_

3 for  _$i\in[1,N],t\in[0,T]$_ do

4 $\epsilon_{t},y_{t},\hat{x}_{t},a_{t}\leftarrow$ Sample observations from
PEM (5-6), track via (3) and choose actions by solving (7)

5 Append $\langle x_{t},a_{t}\rangle$ to trajectory $\tau^{(i)}$

// Maintain work-list per trajectory (Alg 1)

6

7 $\textit{w-list}^{(i)}_{t+1}\leftarrow$ upd-wl(_$\textit{w-list}^{(i)}_{t}$
, $\varphi$, $x_{t+1}$_)

$L^{(i)}_{t+1}\leftarrow\textit{w-list}^{(i)}_{t+1}[\varphi][0]$ // Robustness
of $\tau^{(i)}_{[0:t+1]}$

8

9 $x_{t+1}\leftarrow$ Step forward simulation

10 Sort $\\{L^{(0)}_{T}\dots L^{(N)}_{T}\\}$ then set $\gamma_{0}\leftarrow
L^{(K)}_{T}$

11 $m\leftarrow 0$

12 while  _$\gamma_{m} >\gamma$_ do

13 $m\leftarrow m+1$

14 Discard all trajectories trajectories $\tau^{(i)}$ where
$L^{(i)}_{T}\geq\gamma_{k}$

15 $\mathcal{I}_{k}\leftarrow$ Indices of remaining un-discarded trajectories

16 for  _$i\in[0,N]\setminus\mathcal{I}_{k}$_ do

17 Select a random $j\in{I}_{k}$

18 Find the first time step $t^{\prime}$ where
$L^{(j)}_{t^{\prime}}<\gamma_{k}$

19 $\tau^{(i)}_{[0{:}t^{\prime}]},L^{(i)}_{[0{:}t^{\prime}]}\leftarrow$ Copy
values from $\tau^{(j)}$

20 $\tau^{(i)}_{[t^{\prime}{:}T]},L^{(i)}_{[t^{\prime}{:}T]}\leftarrow$ Re-
simulate $\tau^{(j)}$ from time $t^{\prime}$

21 Sort $\\{L^{(0)}_{T}\dots L^{(N)}_{T}\\}$ then set $\gamma_{m}\leftarrow
L^{(K)}_{T}$

22 return $\hat{p}_{ams}$ via Eqn (12)

Algorithm 2 Online stl-ams

Now that we can compute and cache stl robustness values for partial
trajectories, we can interleave this with our ams sampler. Alg (2) provides an
overview of our sampling technique for computing the probability that our av
violates an stl specification in a stochastic simulation. It takes as input a
starting state $x_{0}$, specification $\varphi$, failure threshold $\gamma$,
initial simulation amount $N$ and discard rate $K$.

Lines (1-6) generate $N$ initial trajectories by simulating perceptual
observations, control actions, and forward dynamics as outlined in sections
(II-A-II-B). Lines (7-8) track the stl robustness value of trajectories at
every intermediate time step by maintaining up-to-date work-lists as described
in Alg (1). Lines (9, 17) adaptively set discard thresholds $\gamma_{m}$ such
that the $K$ safest trajectories are discarded at each stage. To replenish
those discarded trajectories, lines (14-16) copy one of the remaining un-
discarded $\tau^{(i)}$ up until the first time step $j$ where the robustness
value of the partial trajectory falls below $\gamma_{m}$. We then re-simulate
starting from $j$ to produce a new trajectory. We repeat this process until
the discard threshold $\gamma_{m}$ falls below the desired failure threshold
$\gamma$, then calculate a final estimate $\hat{p}_{\textit{ams}}$ using Eq
(12).

##  IV Experiments

The following experiments use our running example of an av lane-change
maneuver to evaluate our method. They demonstrate that Alg (2) can provide
failure estimates for a full black-box av-system across multiple common
traffic rules. When compared to baselines, we find Alg (2) provides more
accurate failure estimates in fewer simulations. Further, we investigate how
sampling performance differs across discard-rates and rule types. Fig (1)
shows our CommonRoad simulation setup [2]: The ego starts in the centre-lane
at $15m$ with velocity $20\ m/s$. Its primary goal is to track a reference
velocity of $v_{g}$=$30m/s$. Three obstacles impede it—a static obstacle at
$40m$ (forcing the ego to change lanes); A centre-lane vehicle at $50m$ with
velocity $5\ m/s$, which cuts into the overtake lane after $0.6$ seconds
(slowing the ego or forcing a lane change); a right-lane vehicle at $50m$,
velocity $10m$, which merges after $1$ second (preventing the ego from slowing
abruptly).222Full scenario specification and code at
github.com/craigiedon/CommonRules Simulations last $T$=$40$ steps ($4s$,
$\Delta t$=$0.1s$). Dynamics evolve according to a kinematic single-track
model [40].

Our av-system under test uses a standard pre-trained lidar-based obstacle
detector—OpenPCDet’s _Multi-Head PointPillar_ [26, 47]. As described in
Section (II-A), we train a surrogate pem to replicate the behaviour of this
detector in simulation. The pem is composed of two separate Gaussian
Processes: The first is a binary classifier, which predicts whether the lidar
perception system would have failed to detect a given obstacle. The second is
a regression model, which predicts how much noise would typically be added to
the true location of an obstacle’s bounding box. The gps were fitted using
_Pyro_ [4], with an rbf kernel and sparse variational regression with 100
inducing points. As training data for fitting the gps, we used $65500$ entries
from the NuScenes Lidar Validation Set.333https://www.nuscenes.org/nuscenes
The gp input features were a 7-d vector—the x/y obstacle position, rotation,
length/width/height dimensions, obstacle “visibility category” (where “$1$”
means $\leq 40\%$ occlusion; “$4$” means $\geq 80\%$). The gp outputs consist
of a 1-d binary variable for successful/unsuccessful obstacle detection, and a
3-d real-valued variable for the offsets between the obstacle’s true
x/y/rotation and the PointPillar estimate.

For tracking and predicting future vehicle locations, we used the interacting
multiple models (imm) filter for lane-changes from [7]. This method operates
similarly to a typical kalman filter, but instead of making estimates based on
a single model, it maintains estimates from multiple models (i.e., for whether
the vehicles will stay in the current lane, switch to the left lane, or switch
to the right lane) and merges those estimates based on each model’s current
likelihood.

For model predictive control (Eq (7)) we use the lane-change controller from
[29]. At a high level, its cost function $J$ is comprised of 8 sub-goals:
Reach a target destination; track a reference velocity, minimize acceleration,
turn velocity, jerk and heading angle; stay close to the centre of the nearest
lane; and avoid entering the “potential field” of other obstacles. To solve
(7), we used the Gurobi [21] optimizer. To ensure a feasible control action is
always available, we first pre-solve a convex simplification of (7) with CVXPy
[14] (following [18]). For full implementation details of the cost-function
sub-goals (and the weights used to balance them) see [29] and the associated
code for our paper.

TABLE I: Interstate traffic rules (Predicate definitions in [32]). Rule |  Description | STL
---|---|---
$\varphi_{1}$ |  Safe Dist from Vehicles | $\begin{aligned} &\square_{[0,\infty]}(\texttt{in\\_same\\_lane}(x^{\textit{ego}},x^{o_{i}})\wedge\texttt{in\\_front\\_of}(x^{\textit{ego}},x^{o})\wedge\\\ &\ \neg O_{[0,t_{\textit{cut}}]}(\texttt{cut\\_in}(x^{o},x^{\textit{ego}})\wedge H_{[1,\infty]}(\neg\texttt{cut\\_in}(x^{o},x^{\textit{ego}})))\\\ &\ \implies\texttt{keeps\\_safe\\_distance\\_prec}(x^{\textit{ego}},x^{o}))\end{aligned}$
$\varphi_{2}$ |  Unnecessary Braking | $\square_{[0,\infty]}\left(\neg\texttt{unnecessary\\_braking}(x^{\textit{ego}},\\{x^{o_{1},\dots,o_{3}}\\})\right)$
$\varphi_{3}$ |  Preserve Traffic Flow | $\begin{aligned} &\square_{[0,\infty]}(\neg\texttt{slow\\_leading\\_vehicle}(x^{\textit{ego}},x^{o_{1}\dots o_{3}})\implies\texttt{preserves\\_flow}(x^{\textit{ego}}))\end{aligned}$
$\varphi_{4}$ |  Don’t Drive Faster than Left Traffic | $\begin{aligned} &\square_{[0,\infty]}(\texttt{left\\_of}(x^{o_{i}},x^{\textit{ego}})\wedge\texttt{drives\\_faster}(x^{\textit{ego}},x^{o_{i}})\\\ &\implies(\texttt{in\\_slow\\_traffic}(x^{o_{i}},x^{\\{o_{1},\dots o_{3}\\}\setminus x^{o_{i}}})\wedge\texttt{slightly\\_higher\\_speed}(x^{\textit{ego}},x^{o_{i}}))\\\ &\quad\vee(\texttt{on\\_access\\_ramp}(x^{\textit{ego}})\wedge\texttt{on\\_main\\_carriageway}(x^{o_{i}})))\end{aligned}$

We test our sampling method with respect to 4 formalizations of rules from the
Vienna Convention on Road Traffic (taken from [32]). Table (I) shows the stl
for each rule. Full definitions of individual predicates (in_same_lane,
drives_faster etc.) are in [32], but high level descriptions of each rule are
as follows: $\bm{\varphi_{1}}$—maintain a minimum distance from vehicles in
front (proportional to vehicle speed). If a vehicle “cuts in” from an adjacent
lane, the ego gets $t_{\textit{cut}}$ seconds to re-establish distance.
$\bm{\varphi_{2}}$—never drop acceleration below “unnecessary” levels
(relative to vehicles in front). $\bm{\varphi_{3}}$—velocity should never fall
below some minimum level (unless stuck in traffic). $\bm{\varphi_{4}}$—do not
exceed the speed of left-lane vehicles unless merging from an access lane, or
left-lane traffic is slow moving.

Our algorithm (listed as STL-AMS below) uses $N$=$250$ starting simulations, a
discard amount of $K$=$25$, and final failure threshold of $\gamma=0$. We
compare against three baselines: First, a Monte-Carlo sampler (Raw-MC), which
runs $N$ simulations, estimating via (10). Second, an importance sampler with
a fixed proposal (Imp-Naive). We choose a proposal distribution which
deliberately fails to detect 50% of obstacles, and applies gaussian noise
($\mu$ = $0$, $\sigma^{2}$= $1$) to the bounding boxes of those it does
detect. Third, a neural network-based importance sampler with an adaptive
proposal learned via the _cross-entropy method_ (Imp-CE) from [22]. Its inputs
and outputs are the same as those of the gp-pem described previously. Similar
to ams, adaptive importance sampling proceeds in stages: At each stage $m$
(for a total of $M=10$ stages), $N_{m}$=$250$ trajectories are sampled and
sorted by robustness. Imp-CE then takes the lowest $K$=$0.1*N_{m}$
trajectories444Standard practice sets $K$ in range $[0.01,0.2]*N_{m}$ [11].
and minimizes their KL-divergence under the target pem versus current
proposal. The intuition is that biasing our proposal towards the least robust
trajectories in each stage should train a proposal which samples failure
events with increasing probability.

###  IV-A Results and Discussion

TABLE II: Estimated Failure Probabilities (5 Repetitions) Method | $\varphi_{1}$ | $\varphi_{2}$ | $\varphi_{3}$ | $\varphi_{4}$
---|---|---|---|---
$\text{Raw-}\textsc{mc}_{250}$ | $1.2\text{e-}02\ (\pm 4.0\text{e-}03)$ | $0.0\ (\pm 0.0)$ | $0.0\ (\pm 0.0)$ | $0.0\ (\pm 0.0)$
Imp-Naive |  $1.4\text{e-}08$ $(\pm 2.5\text{e-}08)$ |  $2.8\text{e-}08$ $(\pm 8.4\text{e-}08)$ |  $2.0\text{e-}08$ $(\pm 3.90\text{e-}08)$ |  $0.0$ $(\pm 0.0)$
Imp-ce |  $7.1\text{e-}06$ $(\pm 2.1\text{e-}05)$ |  $7.6\text{e-}22$ $(\pm 2.2\text{e-}21)$ |  $4.4\text{e-}18$ $(\pm 8.9\text{e-}18)$ |  $0.0$ $(\pm 0.0)$
stl-ams (ours) |  $8.8\text{e-}03$ $(\pm 6.2\text{e-}03)$ |  $1.5\text{e-}03$ $(\pm 2.1\text{e-}03)$ |  $4.7\text{e-}03$ $(\pm 3.4\text{e-}03)$ |  $3.5\text{e-}04$ $(\pm 1.8\text{e-}03)$
Ground Truth | $9.1\text{e-}03$ | $2.0\text{e-}03$ | $3.6\text{e-}03$ | $4.8\text{e-}05$

Table (II) shows estimated probabilities per rule. For “ground-truth”, we
computed Eq (10) using $100000$ simulations. Across rules, our method produces
the most accurate estimates.

Raw-MC gives a reasonable estimate for $\varphi_{1}$ (the least rare). For
other rules with lower probability, raw sampling yields zero simulation
failures, resulting in $0\%$ estimates. Imp-Naive produces non-zero estimates
for all specifications except $\varphi_{4}$, but vastly underestimates failure
probability. This underscores the difficulty of fixed proposals—if a proposal
distribution is too far from the original target, the likelihood weights per
sample vary wildly, with estimates dominated by a tiny number of simulations.
Imp-CE also produces unreliable estimates. Figs (3(a)-3(c)) provide insights
why: As learning progresses for $\varphi_{3}$, the number of failures produced
goes up, yet failure probability goes down. This suggests Imp-CE learns to
bias its distribution towards a small set of unlikely failures, rather than
approaching the true failure distribution. For $\varphi_{4}$, robustness
thresholds initially decrease, but flatline around $m$=$6$—no failures found.
We observe that this “flatlining” continues even past the $M=10$ stages
documented in this paper. This highlights how challenging it can be to learn
relevant features given longer horizons and state-spaces.

(a) Robustness Threshold

(b) Log Probs (-$\infty$ below axis)

(c) Failures

Figure 3: Imp-CE baseline performance over 10 stages of proposal learning.

The results of Table (II) are encouraging, but we found Alg (2) was sensitive
to discard rate $K$. Fig (4) shows how threshold levels evolve at each stage
(for $K$ values from $2$ to $225$). For $\varphi_{4}$ (the rarest), we found
that too low or high $K$s caused unacceptable numbers of
“extinctions”[9]—stages where all trajectories have identical robustness,
rendering replenishment impossible.

(a) $K=2$

(b) $K=25$

(c) $K=125$

(d) $K=225$

Figure 4: STL-AMS robustness thresholds by stage.

Experiments demonstrate Alg (2) is viable for accurately estimating
specification failure for a black-box av-system. However, this case study
looks only at a single interstate traffic scenario, and our experiments
necessarily have limitations: We considered perceptual disturbance as the sole
source of simulation stochasticity; vehicle starting configurations remained
fixed. Such experiments could be extended by placing a prior distribution over
starts [35], without altering the method. To trust our estimates, we also
assume our simulator accurately represents reality. Whilst outside this
paper’s scope, clearly this assumption may not hold: Our pem may be an
inaccurate surrogate of the test domain (i.e., it would be beneficial to
incorporate work on _ml -uncertainty calibration_ [20]). Our scenario also had
fixed traffic behaviour, but real traffic is reactive and stochastic. Other
works explore these issues in detail [28]. Finally, while it can be seen as an
advantage that our method adaptively selects an appropriate number of
simulations, and re-uses results from previous simulations, these advantages
complicate comparisons of our method to others in terms of sample efficiency.
In future work, we aim to compare performance across a wider range of
scenarios in terms of _fixed computational effort_ across the full sampling
pipeline.

##  V Related and Future Work

This paper estimates failure probabilities. Similar tasks include
falsification (find one failure) and adaptive stress testing (find the most-
likely failure) [11]. Such tasks do not directly accomplish our goal, but may
contain insights for rapidly guiding initial simulations towards failure
areas. A related task is _synthesis_ —constructing controllers which
explicitly obey $\varphi$ [1]. While synthesis can enforce adherence to
specifications expressed in tractable stl subsets, the perceptual and control
uncertainty in av scenarios means testing remains necessary.

Combining splitting and logic has been attempted previously [23]. Rather than
use stl robustness, such works restrict themselves to heuristic decompositions
of linear temporal logic formulae. This renders them unsuitable for cyber-
physical domains like av.

One of our experiment baselines was importance sampling. Proposals are often
represented by exponential distributions, since those have analytic solutions
[44, 51]. Neural networks have also been used to represent the proposal (as we
did), [34]. Most work considers rarity in the context of vehicle
configurations or behaviour. Our work instead considers rarity in the context
of perceptual disturbances. A sampler category unexplored in this paper are
markov-chain methods [5]. With appropriate assumptions on metric smoothness
and system linearizability, such techniques have shown promise in domains with
long chains of dependent states [46, 43]. While out of scope, integrating such
techniques with online stl monitoring may prove fruitful.

Despite the asymptotic normality of ams, both splitting and sampling lack
guarantees on estimation error for _fixed_ $N$. _Certifiable sampling_
addresses this with _efficiency certificates_ [3]—customized samplers with a
bound on sampling error relative to $N$. However, certification methods depend
on first being able to find the failure regions, so techniques from this paper
remain relevant.

Online and offline algorithms exist to calculate stl robustness [16, 13].
Typically, efficiency is not considered in the context of a sampling regime.
Yet recent advances in online monitoring could be leveraged within ams to
discard infeasible trajectories early. For example, incorporating system
dynamics, or causality [50, 12].

Our experiments target interstate lane changes. Others encode rules for
intersections, and situational awareness [15, 31]. In future work, we aim to
assess sampling effectiveness across this diversity of specifications.

## References

  * [1] Aasi, E., Vasile, C.I., Belta, C.: A control architecture for provably-correct autonomous driving. In: 2021 American Control Conference (ACC). pp. 2913–2918. IEEE (2021)
  * [2] Althoff, M., Koschi, M., Manzinger, S.: Commonroad: Composable benchmarks for motion planning on roads. In: 2017 IEEE Intelligent Vehicles Symposium (IV). pp. 719–726. IEEE (2017)
  * [3] Arief, M., Bai, Y., Ding, W., He, S., Huang, Z., Lam, H., Zhao, D.: Certifiable deep importance sampling for rare-event simulation of black-box systems. arXiv preprint arXiv:2111.02204 (2021)
  * [4] Bingham, E., Chen, J.P., Jankowiak, M., Obermeyer, F., Pradhan, N., Karaletsos, T., Singh, R., Szerlip, P.A., Horsfall, P., Goodman, N.D.: Pyro: Deep universal probabilistic programming. J. Mach. Learn. Res. 20, 28:1–28:6 (2019), http://jmlr.org/papers/v20/18-403.html
  * [5] Botev, Z.I., L’Ecuyer, P., Tuffin, B.: Markov chain importance sampling with applications to rare event probability estimation. Statistics and Computing 23, 271–285 (2013)
  * [6] Bugallo, M.F., Elvira, V., Martino, L., Luengo, D., Miguez, J., Djuric, P.M.: Adaptive importance sampling: The past, the present, and the future. IEEE Signal Processing Magazine 34(4), 60–79 (2017)
  * [7] Carvalho, A., Gao, Y., Lefevre, S., Borrelli, F.: Stochastic predictive control of autonomous vehicles in uncertain environments. In: 12th international symposium on advanced vehicle control. vol. 9 (2014)
  * [8] Cérou, F., Guyader, A.: Fluctuation analysis of adaptive multilevel splitting. The Annals of Applied Probability 26(6), 3319 – 3380 (2016). https://doi.org/10.1214/16-AAP1177, https://doi.org/10.1214/16-AAP1177
  * [9] Cérou, F., Guyader, A., Rousset, M.: Adaptive multilevel splitting: Historical perspective and recent results. Chaos: An Interdisciplinary Journal of Nonlinear Science 29(4) (2019)
  * [10] Corso, A., Lee, R., Kochenderfer, M.J.: Scalable autonomous vehicle safety validation through dynamic programming and scene decomposition. In: 2020 IEEE 23rd International Conference on Intelligent Transportation Systems (ITSC). pp. 1–6. IEEE (2020)
  * [11] Corso, A., Moss, R., Koren, M., Lee, R., Kochenderfer, M.: A survey of algorithms for black-box safety validation of cyber-physical systems. Journal of Artificial Intelligence Research 72, 377–428 (2021)
  * [12] Deng, Z., Eshima, S.P., Nabity, J., Kong, Z.: Causal signal temporal logic for the environmental control and life support system’s fault analysis and explanation. IEEE Access 11, 26471–26482 (2023)
  * [13] Deshmukh, J.V., Donzé, A., Ghosh, S., Jin, X., Juniwal, G., Seshia, S.A.: Robust online monitoring of signal temporal logic. Formal Methods in System Design 51, 5–30 (2017)
  * [14] Diamond, S., Boyd, S.: CVXPY: A Python-embedded modeling language for convex optimization. Journal of Machine Learning Research 17(83), 1–5 (2016)
  * [15] Dokhanchi, A., Amor, H.B., Deshmukh, J.V., Fainekos, G.: Evaluating perception systems for autonomous vehicles using quality temporal logic. In: Runtime Verification: 18th International Conference, RV 2018, Limassol, Cyprus, November 10–13, 2018, Proceedings 18. pp. 409–416. Springer (2018)
  * [16] Donzé, A., Ferrere, T., Maler, O.: Efficient robust monitoring for stl. In: Computer Aided Verification: 25th International Conference, CAV 2013, Saint Petersburg, Russia, July 13-19, 2013. Proceedings 25. pp. 264–279. Springer (2013)
  * [17] Duvenaud, D.: Automatic model construction with Gaussian processes. Ph.D. thesis, University of Cambridge (2014)
  * [18] Eiras, F., Hawasly, M., Albrecht, S.V., Ramamoorthy, S.: A two-stage optimization-based motion planner for safe urban driving. IEEE Transactions on Robotics 38(2), 822–834 (2021)
  * [19] Fainekos, G.E., Pappas, G.J.: Robustness of temporal logic specifications for continuous-time signals. Theoretical Computer Science 410(42), 4262–4291 (2009)
  * [20] Guo, C., Pleiss, G., Sun, Y., Weinberger, K.Q.: On calibration of modern neural networks. In: International conference on machine learning. pp. 1321–1330. PMLR (2017)
  * [21] Gurobi Optimization, LLC: Gurobi Optimizer Reference Manual (2023), https://www.gurobi.com
  * [22] Innes, C., Ramamoorthy, S.: Testing rare downstream safety violations via upstream adaptive sampling of perception error models. In: 2023 IEEE International Conference on Robotics and Automation (ICRA). pp. 12744–12750. IEEE (2023)
  * [23] Jegourel, C., Legay, A., Sedwards, S.: An effective heuristic for adaptive importance splitting in statistical model checking. In: International Symposium On Leveraging Applications of Formal Methods, Verification and Validation. pp. 143–159. Springer (2014)
  * [24] Juneja, S., Shahabuddin, P.: Rare-event simulation techniques: An introduction and recent advances. Handbooks in operations research and management science 13, 291–350 (2006)
  * [25] Kim, Y., Kochenderfer, M.J.: Improving aircraft collision risk estimation using the cross-entropy method. Journal of Air Transportation 24(2), 55–62 (2016)
  * [26] Lang, A.H., Vora, S., Caesar, H., Zhou, L., Yang, J., Beijbom, O.: Pointpillars: Fast encoders for object detection from point clouds. In: Proceedings of the IEEE/CVF conference on computer vision and pattern recognition. pp. 12697–12705 (2019)
  * [27] Lemire, D.: Streaming maximum-minimum filter using no more than three comparisons per element. arXiv preprint cs/0610046 (2006)
  * [28] Li, J., Sun, L., Zhan, W., Tomizuka, M.: Interaction-aware behavior planning for autonomous vehicles validated with real traffic data. In: Dynamic Systems and Control Conference. vol. 84287, p. V002T31A005. American Society of Mechanical Engineers (2020)
  * [29] Liu, C., Lee, S., Varnhagen, S., Tseng, H.E.: Path planning for autonomous vehicles using model predictive control. In: 2017 IEEE Intelligent Vehicles Symposium (IV). pp. 174–179. IEEE (2017)
  * [30] Louvin, H., Dumonteil, E., Lelièvre, T., Rousset, M., Diop, C.M.: Adaptive multilevel splitting for monte carlo particle transport. In: EPJ Web of Conferences. vol. 153, p. 06006. EDP Sciences (2017)
  * [31] Maierhofer, S., Moosbrugger, P., Althoff, M.: Formalization of intersection traffic rules in temporal logic. In: 2022 IEEE Intelligent Vehicles Symposium (IV). pp. 1135–1144. IEEE (2022)
  * [32] Maierhofer, S., Rettinger, A.K., Mayer, E.C., Althoff, M.: Formalization of interstate traffic rules in temporal logic. In: 2020 IEEE Intelligent Vehicles Symposium (IV). pp. 752–759. IEEE (2020)
  * [33] Maler, O., Nickovic, D.: Monitoring temporal properties of continuous signals. In: International Symposium on Formal Techniques in Real-Time and Fault-Tolerant Systems. pp. 152–166. Springer (2004)
  * [34] Müller, T., McWilliams, B., Rousselle, F., Gross, M., Novák, J.: Neural importance sampling. ACM Transactions on Graphics (ToG) 38(5), 1–19 (2019)
  * [35] O’Kelly, M., Sinha, A., Namkoong, H., Tedrake, R., Duchi, J.C.: Scalable end-to-end autonomous vehicle testing via rare-event simulation. Advances in neural information processing systems 31 (2018)
  * [36] Pandharipande, A., Cheng, C.H., Dauwels, J., Gurbuz, S.Z., Ibanez-Guzman, J., Li, G., Piazzoni, A., Wang, P., Santra, A.: Sensing and machine learning for automotive perception: A review. IEEE Sensors Journal 23(11), 11097–11115 (2023). https://doi.org/10.1109/JSEN.2023.3262134
  * [37] Piazzoni, A., Cherian, J., Dauwels, J., Chau, L.P.: Pem: Perception error model for virtual testing of autonomous vehicles. IEEE Transactions on Intelligent Transportation Systems (2023)
  * [38] Qi, C.R., Yi, L., Su, H., Guibas, L.J.: Pointnet++: Deep hierarchical feature learning on point sets in a metric space. Advances in neural information processing systems 30 (2017)
  * [39] Rajabli, N., Flammini, F., Nardone, R., Vittorini, V.: Software verification and validation of safe autonomous cars: A systematic literature review. IEEE Access 9, 4797–4819 (2020)
  * [40] Rajamani, R.: Vehicle dynamics and control. Springer Science & Business Media (2011)
  * [41] Rawlings, J.B., Mayne, D.Q., Diehl, M.: Model predictive control: theory, computation, and design, vol. 2. Nob Hill Publishing Madison, WI (2017)
  * [42] Riedmaier, S., Ponn, T., Ludwig, D., Schick, B., Diermeyer, F.: Survey on scenario-based safety assessment of automated vehicles. IEEE access 8, 87456–87477 (2020)
  * [43] Scher, G., Sadraddini, S., Tedrake, R., Kress-Gazit, H.: Elliptical slice sampling for probabilistic verification of stochastic systems with signal temporal logic specifications. In: Proceedings of the 25th ACM International Conference on Hybrid Systems: Computation and Control. pp. 1–11 (2022)
  * [44] Schmerling, E., Pavone, M.: Evaluating trajectory collision probability through adaptive importance sampling for safe motion planning. arXiv preprint arXiv:1609.05399 (2016)
  * [45] Schwarting, W., Alonso-Mora, J., Rus, D.: Planning and decision-making for autonomous vehicles. Annual Review of Control, Robotics, and Autonomous Systems 1, 187–210 (2018)
  * [46] Sinha, A., O’Kelly, M., Tedrake, R., Duchi, J.C.: Neural bridge sampling for evaluating safety-critical autonomous systems. Advances in Neural Information Processing Systems 33, 6402–6416 (2020)
  * [47] Team, O.D.: Openpcdet: An open-source toolbox for 3d object detection from point clouds. https://github.com/open-mmlab/OpenPCDet (2020)
  * [48] Thrun, S.: Probabilistic robotics. Communications of the ACM 45(3), 52–57 (2002)
  * [49] Williams, C.K., Rasmussen, C.E.: Gaussian processes for machine learning, vol. 2. MIT press Cambridge, MA (2006)
  * [50] Yu, X., Dong, W., Yin, X., Li, S.: Online monitoring of dynamic systems for signal temporal logic specifications with model information. In: 2022 IEEE 61st Conference on Decision and Control (CDC). pp. 1553–1559. IEEE (2022)
  * [51] Zhao, D., Lam, H., Peng, H., Bao, S., LeBlanc, D.J., Nobukawa, K., Pan, C.S.: 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 (2016)