Detecting Traffic Incidents Using Persistence Diagrams

: We introduce a novel methodology for anomaly detection in time-series data. The method uses persistence diagrams and bottleneck distances to identify anomalies. Speciﬁcally, we generate multiple predictors by randomly bagging the data (reference bags), then for each data point replacing the data point for a randomly chosen point in each bag (modiﬁed bags). The predictors then are the set of bottleneck distances for the reference/modiﬁed bag pairs. We prove the stability of the predictors as the number of bags increases. We apply our methodology to trafﬁc data and measure the performance for identifying known incidents. windows near each incident in the data sets and the thresholds we determined from training data sets. In each phase, we use the thresholds determined by the training data to classify the windows in the corresponding testing data. We measure the performance of these predictions using precision, recall and an F-score.


Introduction
Traffic incidents are severely detrimental to society in terms of financial costs which are estimated in the U.S. by the National Highway Traffic Safety Administration [1]. Consequently, an important focus of data analysis concerns detecting incidents from traffic data for the management of response to accidents which can have significant benefits for society [2]. The type of data we consider is a time series of volumetric traffic counts. We propose a novel methodology for analyzing this data for the purpose of identifying traffic incidents.. Our approach is to view the identification of incidents as a problem of anomaly detection within the time-series data. Our method uses tools from statistical analysis-bootstrap aggregation or bagging-and from topological data analysis (TDA)-persistence diagrams-to form an ensemble of predictors to determine whether data points are normal or anomalous. To each data point and random bag, we associate two persistence diagrams, one reference and one modified. The predictors then consist of a score-the bottleneck distance between the diagrams-and for each data point, the set of scores are aggregated into several summary statistics. We then identify data points as being incidents or anomalies by percentile scores of the summary statistics.
Our algorithm using randomized bagging of the data and resultant persistence diagrams as a feature extractor for anomaly detection can be viewed as a semi-supervised learning algorithm. Indeed, our method trains the summary statistics on a small amount of labeled data. Moreover, our algorithm can be applied to any vectorized data and, thus, can be adapted for many other data analytic problems.

Description of the Data and Challenge Problem
The problem we address is to identify incidents in a data set of volumetric traffic counts obtained from multiple inductive loop road sensors, [3], supplied by the California Department of Transportation. The counts were aggregated over 5 min intervals, and we had access to one full calendar year of traffic counts for each sensor. See Table 1 for a sample from the data. We apply our method to this data set in Section 4. Our task to identify incidents was part of an anomaly detection problem hosted by the Joint Algorithms for Threat Detection (ATD) and Algorithms for Modern Power Systems (AMPS) Annual Workshop; the challenge problem that we specifically addressed was curated for the ATD program [4]. The challenge problem consisted of two phases: in Phase 1, we were asked to detect incidents within data provided by 10 road sensors that were spatially located to be independent from one another (we were not informed of the locations); in Phase 2, we were asked to repeat Phase 1 provided additional information on the location of the sensors-see Figure 1. For each phase, we were supplied with a training data set and a testing data set with labeled incidents that were hand-selected by the ATD program coordinators. The training data set for Phase 1 contained 9 labeled incidents, and the training data set for Phase 2 contained 10 labeled incidents. There were numerous unlabeled incidents in the training data sets. The timeline of the challenge problem was as follows: Supplementary materials in the form of data sets and code are available at the following repository: https://bitbucket.org/esweber/randomized-persistence-landscapes/src/master/.
Data sets for the challenge problem were originally released as files in MATLAB. We have converted the data to comma-separated values files and all experiments discussed in the paper can be reproduced by executing our code in R-Studio. Our solutions were submitted as a binary output for each 5 min interval with 1 indicating an incident and 0 otherwise. Three evaluation metrics were then used to evaluate the performance of our algorithm: precision, recall, and an F-score.

Survey of Literature
Detecting incidents in traffic data is important for the safety and welfare of society [2]. Significant efforts have been made to automate incident detection [5][6][7][8]. Traditional automatic incident detection data sources include radar sensors [9], video cameras [10], probes [11]. Loop detectors are also a very common source of traffic data [3,12], which is the source of the traffic data analyzed in our method.
Detecting incidents in traffic data can be formulated as an anomaly detection problem. Approaches in this formulation include: dynamic clustering [13] or Support Vector Machines [14] to detect anomalous trajectories, autoencoders to learn regularity of video sequences [15], motion reconstruction [16], and convolutional neural networks to predict congestion in traffic networks [17]. In this paper, we also view incident detection as an anomaly detection problem, with our approach using TDA to identify outliers. For an overview of anomaly detection, see [18].
Topological data analysis has been used previously in the context of traffic data: [19] uses TDA as a model for tracking vehicles and [20] uses TDA to understand individual travel behaviors. TDA has also been used previously for anomaly detection in [21,22]. TDA provides a framework for studying the "shape" of data in geometric and topological terms rather than a solely statistical framework [23][24][25].
Our algorithm is motivated by the intuition that if we take a data set of vectors that represent typical behavior and one of the vectors were replaced with an anomalous observation, the topology of the data set would be significantly altered. This intuition is supported by randomized persistence diagrams [46,47]. For example, [48] finds that samples from a manifold corrupted by noise with a Gaussian distribution will not in general significantly deviate from the persistence diagram of the underlying manifold, but other distributions could. Further results in [49] emphasize this by showing that extending the stability results for noise with a Gaussian distribution in [48] to other noise distributions should not "be expected".

Main Contributions
Our main contributions in this paper are two-fold. First, we introduce a novel machine learning algorithm that uses TDA for detecting anomalies in time-series data. This is done by establishing a baseline of typical data distribution to serve as a reference through the method of bagging, and deviations from the reference for various data points are measured using the bottleneck distance from TDA. This procedure is repeated multiple times to reduce the variation that naturally arises from random sampling. The algorithm requires selection of several hyperparameters. Second, we address the problem of identifying incidents in traffic data that consists only of traffic counts. Our data set does not include any functional data, such as velocity, nor does it include video feed which are the most common data sources used in traffic incident detection problems.

Outline
The rest of the paper is organized as follows. In Section 2 we present the necessary background for understanding persistence diagrams, the fundamental tool from TDA that our algorithm implements. In Section 3 we present the algorithm. In Section 4 we present the results of our algorithm when applied to the traffic data set from the ATD challenge problem.

Topological Data Analysis (TDA)
Here, we provide the necessary background for understanding the TDA that our algorithm implements. For more details on the subject, see [50,51]. The appeal of TDA is that it offers techniques from topology to ascertain the "shape" of high-dimensional data sets that is not accessible by other methods. In particular, persistent homology is ubiquitous in TDA for measuring certain topological features of a space (e.g., connectivity, holes, voids). These features are often summarized using Betti numbers or a persistence diagram [23]. Our focus will be on the latter and the bottleneck distance which serves as a metric on a set of those diagrams.

Persistence Diagrams
Persistence diagrams are computed based on the notion of a simplicial complex. It is natural to think of a simplicial complex in a Euclidean space, but a more abstract definition adequately serves our purposes.

Definition 1 (Abstract Simplicial Complex).
A finite collection of sets A is called an abstract simplicial complex if A is closed under subset inclusion, i.e., if β ∈ A whenever α ∈ A and β ⊂ α. Furthermore, we call α ∈ A a combinatorial p-simplex if |α| = p + 1.
In our consideration, each simplex is a subset of data points. There are several useful methods by which to construct simplicial complexes from data. We have chosen the Vietoris-Rips (VR) complex for its computational convenience. Definition 2 (Vietoris-Rips Complex). Let Y be a finite subset of data points contained in some metric space X. For every r ≥ 0, the Vietoris-Rips Complex of Y at scale r is defined as We note that the VR complex is ascending in the sense that VR r 0 (Y) ⊂ VR r (Y) for r 0 ≤ r. Thus, as the radius r increases, we produce a filtration of simplicial complexes on top of the data points in Y. This filtration yields a persistent homology that quantifies the apparent topological features for different dimensions. 0-, 1-, and 2-dimensional features refer to connected components, holes, and voids, respectively. Higher-dimensional features would be the analogues of these features. The persistence diagram quantifies the importance of these features in terms of their persistence over various radii. Persistence gives us a sequence of homology groups for each dimension of the ambient metric space containing the data. By identifying the radii in the sequence when the topological features appear and disappear we obtain a collection of birth and death times for each feature of each dimension.
The birth-death pairs as a multiset in R 2 make up the persistence diagram corresponding to Y. The persistence diagrams in our algorithm only contain information on connected components, the zero-dimensional features. We were unable to find a meaningful use for persistence diagrams computed using higher-dimensional features.

Bottleneck Distance
We use the bottleneck distance to make comparisons between two persistence diagrams. The bottleneck distance can be thought of more generally as a pseudometric on multisets P and Q that quantifies the cost of mapping P to Q. To be precise, suppose P and Q are plain multisets of points in the extended plane R 2 . We say M ⊂ P × Q is a partial matching between P and Q if it satisfies the following: • For each q ∈ Q, there exists at most one p ∈ P such that (p, q) ∈ M. • For each p ∈ P there exists at most one q ∈ Q such that (p, q) ∈ M.
We say that We define the cost function for each matching M to be Definition 3 (Bottleneck Distance). Let M(P, Q) be the set of all partial matchings between P and Q. We define the bottleneck distance between P and Q as

Remark 1.
We refer to the bottleneck distance as a pseudometric because there are multisets P and Q such that P = Q yet W(P, Q) = 0. Take for example P = Q 2 \∆ and Q = (Q + √ 2) 2 \∆. Although c(M) > 0 for every M ∈ M(P, Q), we find W(P, Q) = 0 by taking an infimum over M(P, Q). It has been shown that in the case where P and Q are finite multisets in R 2 \∆ then W(P, Q) = 0 implies P = Q, [52].

Description of the Algorithm
Our algorithm is motivated by the intuition that if we take a data set of vectors that represent typical behavior and one of the vectors were replaced with an anomalous observation, the topology of the data set would be significantly altered. The topology of the set is understood by computing persistence diagrams of randomly chosen subsets β ⊂ {x i } m i=1 . When we modify a bag β by randomly replacing one of its vectors with x j we denote the resulting modified bag by β (j) . Figure 2 depicts the data points in the reference and modified bags. We then let D and D (j) denote the persistence diagrams for vectors in β and β (j) , respectively. Figure 3 depicts the persistence diagrams for the reference and modified bags for a specfic data point.  The bottleneck distance quantifies how much the topology of a set has changed. We say there is evidence that There are some undesirable cases that may occur due to the random sampling used to form β. It may be that the observations contained in β do not adequately represent the entire data set, or if the reference bag contains an anomaly, we might replace an anomaly with a non-anomalous vector when forming the modified bag. To mitigate issues such as these, the algorithm uses multiple bags of the same size. We select hyperparameters S, N ∈ N and repeat N times the process of randomly choosing S vectors to form β. This creates a collection of reference bags {β k } N k=1 . For each of the N bags, we form m distinct modified bags {β (j) k } m j=1 by randomly choosing y ∈ β k and replacing it with x j . We calculate the summary statistics mean, median, and standard deviation for the bottleneck distances, and choose a function f : R 3 → {0, 1} that identifies anomalies from the summary statistics. The entire algorithm is presented in Algorithm 1.

Algorithm 1: Anomaly Detection using Persistence Diagrams of Random Subsamples
Compute and store d j,k = W(D k , D The summary statistics of the bottleneck distances {d j,k } N k=1 can be used in various ways depending on the application. In our application to traffic data, we set thresholds for each summary statistic based on training data. Thus, we set thresholds τ 1 , τ 2 , τ 3 for the three summary statistics, In addition to those previously mentioned, we may also encounter the issue when the modified bag β (k) is formed by replacing the same vector x k . This would make an anomaly look much less anomalous based on our bottleneck distance measurement. Fortunately, if N is chosen to be sufficiently large, then the effect can be very small since the summary statistics converge in probability. Let us formalize the probability distribution here. The sample space consists of ordered pairs We place the uniform distribution on Ω and denote this by P . Theorem 1. Let d j , s j , and h j be denote the mean, standard deviation, and median respectively of bottleneck distance associated with x j as computed in Algorithm 1. Then there exist µ j , σ j , η + j , and η − j such that Proof. Define the random variable d j on Ω by To prove the second limit, we start by writing as N tends to infinity. Adding up the terms in (1), this means s j → σ j in probability. For the third equation, first define Since 1 2 − α > 0, we apply the Weak Law of Large Numbers to b j,1 , ..., b j,N to see that . Much like what we had in (2), we find Since 1 2 − δ > 0, we again have from the Weak Law of Large Numbers applied to b j,1 , ..., b j,N that

Remark 2.
Since we are sampling from a discrete probability distribution, it is not guaranteed that the distribution producing the {d j,k } has a true median. Hence we must settle for convergence of the sample median to some interval. In the case where there exists η such that P[d j,k ≥ η], P[d j,k ≤ η] ≥ 1 2 , then η + j = η − j = η and the sample median converges in probability to η. We define the sample median the usual way where if X 1 ≤ X 2 ≤ ... ≤ X N , the sample median would be X (N+1)/2 if N is odd or 1 2 (X N/2 + X N/2+1 ) if N is even.

Remark 3.
We note that the values guaranteed in Theorem 1 are data dependent, and that for the purposes of Theorem 1 as well as Algorithm 1, the data are fixed.

Results
In this section, we present our contribution to the ATD challenge problem through the application of our method on traffic data collected from major highways in the Sacramento area. The problem was divided into two phases-depending on whether the location information of the sensors was provided-where each phase consisted of a training and testing data set. The objective of each phase of the challenge problem was to predict the time and location of hand-selected traffic incidents in the testing data using the training data which had very few labeled incidents. The data was collected as a count of cars that passed a given sensor during each 5-min interval throughout the 2017 calendar year. An example of this can be seen in Table 1.
The training sets included details on certain incidents reported during the year that include the nearest sensor, the timestamp, and the duration of each incident. In each data set there are a few instances in which a particular sensor was not operating so that there are no counts reported during those five-minute intervals. Table 2 contains additional information on the data sets that were provided. To apply the algorithm to the volumetric data, we considered the embedding of the data in R 12 as sequences of 12 consecutive 5-min counts from a sensor. This means each vector represents 1 hour's worth of counts. We index each vector by a 3-tuple (p, t, w). We use p to denote the sensor ID from which the counts were collected. We let t = (t, d), where t denotes the starting time of the five-minute window corresponding to the first of the 12 counts and d ∈ {1, ..., 7} indicates the day of the week with d = 1 corresponding to Sunday, d = 2 corresponding to Monday, etc. We let w ∈ {1, ..., 52} denote the week in which these counts were collected. (We note that in 2017, there were 53 Sundays, and so for d = 1, w ∈ {1, ..., 53}.) The vectors were then sorted by sensor, starting time, and day of week, meaning the set of all vectors was partitioned into smaller data sets of 52 vectors. We let D p,t denote the collection of count vectors from sensor p collected at the time (hour and day of the week) corresponding to t. For two different starting times from the same weekday, say t = (t, d) and t * = (t * , d), it is possible that the components of vectors in D p,t overlap with those of vectors in D p,t * . For example, if t = 8:00 AM, and t * = 8:05 AM, then the vectors in D p,t * are nearly the same vectors as in D p,t , but the entries are shifted left by one component and the final entry is the count from 9:00-9:05 AM. Since there are 288 five-minute intervals throughout the day and we required 12 consecutive counts to form each vector, there are 277 possible vectors of counts each day.
After sorting the vectors, we applied Algorithm 1 to each collection of counts D p,t with S = 30 and N = 30. Thus, each vector in D p,t , was assigned 30 bottleneck distances. For each time window (p, t, w), we recorded the mean, median, and standard deviation of these bottleneck distances which we denote by d p,t,w , h p,t,w , and s p,t,w , respectively.
As in [53], we expected different periods of the day to exhibit different behavior. Therefore, after obtaining the summary statistics from Algorithm 1, our vectors were classified again according to three factors: day of week, sensor, and time of day. Time of day had 5 levels based on the time of the vector's first count. The levels were Early Morning, Morning, Mid Day, Evening, and Late Evening, and they corresponded to one-hour windows with start times in the ranges: The windows were ranked by each of their summary statistics within their classification. These rankings were given in the form of a percentile and used as anomaly scores.
We expected that traffic patterns throughout the year should be fairly consistent given a particular sensor, day of the week, and period of the day; however, it is possible that traffic behavior could be unusual for an entire day due to a holiday, road construction, or a weather event. Then, according to this window classification, it cannot be readily ascertained if a particular window with a relatively large mean bottleneck distance is due to something acute, such as a traffic incident, or to something predictable, like Labor Day or a major sporting event. We took this into account by reclassifying the windows by sensor, time of day, and day of the year. We then ranked the windows a second time within each of these treatment groups. Each window was assigned 6 percentile scores, based on three different summary statistics each ranked two different ways.
We describe our procedure using an example from the Phase 2-Train data. To measure the likelihood that there was an incident occurring near sensor p 0 := S314402 on Monday, 6 February 2017, during the window from 7:55 AM to 8:55 AM, we set t 0 = (7 : 55, 2). We apply our algorithm to D p 0 ,t 0 . See Figure 3 for an example of the persistence diagrams related to this particular timestamp. If an incident occurred on February 6th, the sixth Monday of 2017 in the time window classification, we would expect that d p 0 ,t 0 ,6 would be large compared to the rest of the collection {d p 0 ,t 0 ,w } 52 w=1 . When this is done, we find that d p 0 ,t 0 ,6 = 41.95, which is the second largest of the 52 average bottleneck distances in {d p 0 ,t 0 ,w } 52 w=1 . Even though this was not the largest average bottleneck distance observed for this sensor, day of week, and time of day, it did rank above the mid-98th percentile among all Mid Day windows according to average bottleneck distance for S314402 on a Monday. Similarly, s p 0 ,t 0 ,6 = 27.56 ranked above the 96th percentile and of standard deviations of bottleneck distance, and h p 0 ,t 0 ,6 = 31.16 ranked just above the 99th percentile of median bottleneck distances observed from Mid Day, Monday windows at sensor S314402. If we consider this window among only the observations from Mid Day on February 6th at sensor S314402, then d p 0 ,t 0 ,6 ranks just above the mid-89th percentile, s p 0 ,t 0 ,6 ranks above the mid-90th percentile, and h p 0 ,t 0 ,6 is at least tied as the highest ranking sample median of its class.
After the six rankings are determined for each window in our data set, we are ready to apply selection criteria to determine which windows overlapped with traffic incidents. The selection consists of six thresholds for the six rankings. The thresholds are determined by the rankings of windows near the starting times of the labeled incidents in the training data using the procedure outlined below:

1.
Identify all windows starting within 30 min of the timestamp of a reported incident at the same sensor where the incident was located.

2.
For each of the 6 types of rankings, identify a minimum ranking needed to include at least one window corresponding to each labeled incident. If all 6 minimum rankings seem too small, the minimum ranking to include one window from all but one incident could be used instead.

3.
Set the first threshold, τ 1 as the largest of the six minimums.

4.
Identify which of the windows found in step 1 satisfy the threshold of τ 1 .

5.
For each of the 5 types of rankings that were not used to determine τ 1 , identify the minimum ranking met by all the windows identified in step 4. 6.
Once the thresholds τ 1 , ..., τ 6 are set, we classify all windows that have all six rankings each above their respective thresholds as incidents. In the rest of this section, we present the results when this procedure is applied to each of the four data sets. We present the rankings of windows near each incident in the data sets and the thresholds we determined from training data sets. In each phase, we use the thresholds determined by the training data to classify the windows in the corresponding testing data. We measure the performance of these predictions using precision, recall and an F-score.

Remark 4.
Our performance evaluation using Algorithm 1 exclude any incidents detected on Sundays since the training data did not include any reported incidents on Sunday.
To provide some comparison, we also apply a standard normal deviates (SND) algorithm to detect incidents in the phase-1 data. The SND algorithm detects an incident in the kth 5-min time window if the count for that window, x k , satisfies |x k − x| > τσ.
We let x andσ respectively denote the mean count and standard deviation of counts from the same treatment group and τ denotes some threshold that applies across all treatment groups. In this case, a treatment group is made up of all counts belonging to the same sensor, on the same day of week, at the same time of day. For example, in Phase-1, there were 52 counts taken at sensor S312425, on Monday mornings at 8:00 AM, {x k } 52 k=1 . We can compute the mean and standard deviation of this treatment group, For a better description of the SND algorithm and its application to incident detection see [11], where incident detection was performed using vehicle velocity rather than traffic flow. The threshold τ was determined using the Phase 1-Train data by computing the deviates for each window occurring during a labeled incident. For each of the 15 incidents, we computed the 85th percentile of the counts that took place during the incident. We chose τ to be the second smallest of the 15 percentiles we computed. Thus, for the phase 1 data, we found τ = 1.06, meaning any 5-min window where the count exceeded 1.06 standard deviations from the mean was detected as an incident.

Data without Sensor Locations
In the Phase 1-Train data set, we were able to find windows near the incidents that had very large average bottleneck distances. When sorted by sensor, day of week, and time of day, all of the 15 labeled incidents overlapped with a one-hour window starting within half an hour of some window that was ranked in the 85th percentile. When sorting these sensors further by day of year, the start time of six of the incidents fell within half an hour of the window with the highest average bottleneck distance in their category. Tables 3 and 4 present the percentile scores of the incidents in the Phase 1-Train dataset. Table 5 contains the 6 thresholds determined using the rankings of windows near the labeled incidents. The quality of the fit is given in Table 6; as a comparison, the quality of fit using the standard normal derivates is given in Table 7. To better describe the quality of the anomaly scores, we provide the receiver operating characteristic (ROC) curve based on the three rankings in Figure 4. We also provide the ROC curve based on the standard normal deviates.
Of course, these performance scores say very little about the method since the thresholds were determined using the incidents in the data set. Rather, we apply the thresholds in Table 5 to classify the Phase 1-Test data. Table 8 has the performance scores from this experiment, whereas Table 9 has the performance for the standard normal deviates on this dataset. Table 3. Percentiles of incidents from the Phase1-Train data when sorted by sensor, day of week, and time of day.  Table 3. Cont.

Data with Sensor Locations
In Phase 2 we were given the locations of the sensors in terms of latitude and longitude. Tables 10 and 11 display the nine incidents reported in Phase 2-Train and the largest percentile recorded for each summary statistic in windows starting within half an hour of the starting time of each incident at the sensor where the incident was reported. Depending on the severity of the incident, it is possible that a traffic incident occurring near one sensor will cause some abnormal behavior in the counts of the adjacent sensors. For example, in the network of sensors used in the Phase 2-Train data, traffic flows directly from S314402 to S318593. A map of all the sensor locations in this data set is provided in Figure 1. If one of the lanes is obstructed near S318593 it is likely to slow down traffic between the two sensors. This might cause the counts from S314402 to look fairly large when compared to S318593. On the other hand, if a motorist was able to anticipate this problem in the traffic a mile ahead, they might deviate their route before even passing S314402. If enough cars did this, there would be a lower count at the first sensor, but the count at S318593 might still be higher because of all the unfortunate cars that got caught between the sensors at the time of the incident. In either scenario, we would expect some unusual behavior when we compare the counts of both sensors together. Not knowing if the difference should be large or small, we apply the random bagging algorithm to the sequence of differences between the counts, with the intuition that if an incident occurred during a particular hour, especially one with lots of traffic, the mean bottleneck distance for that time window would also be large. We refer to the summary statistics obtained from this procedure as adjacency statistics. In Tables 12 and 13 we provide the highest rankings near the reported incidents based on the adjacency statistics.
The most noticeable difference made by the adjacency statistics can be observed by the percentile of the July incident when windows are sorted by sensor, day of year, and time of day. The highest any window near that time ranks by average bottleneck distance of raw counts is in the 66th percentile of Mid Day counts for that day, but if we consider the differences in counts between the two adjacent sensors, we find the highest ranking Mid Day window for July 24th starts within half an hour of the July incident. With the addition of the adjacency statistics, we had 12 percentiles to consider for each window. Table 14 presents the 12 thresholds determined from the incidents in the data for Phase 2-Train. The quality of the fit is given in Table 15. We only report the performance statistics for sensors S314402 and S318593 since those were the only sensors in the data set where incidents were labeled. Table 13. Percentiles of Phase 2-Train incidents when sorted by sensor, day of year, and time of day based on their adjacency statistics.  As in Phase 1, we use the thresholds learned from the training data to classify the windows in the Phase 2-Test data. In this data set, there were 8 sensors. The 8 sensors formed 8 pairs of adjacent sensors, meaning we were able to compute adjacency scores for the windows on every sensor. In Table 16 we present the performance scores for this classification. The data for Phase 2-Test was very different from any of the other data sets. There were 1409 incidents with an average duration of 42 min. Surprisingly, no true incidents were reported at sensor S313386.

Conclusions
Detecting traffic incidents using volumetric data is challenging, as reflected in our performance scores even though our method performed the best of all ATD challenge problem participants [54]. When compared to the SND algorithm, which is a commonly used method for anomaly detection in traffic data, our TDA approach performed slightly better based on the area under the ROC curve. We think there is much room for improvement in using TDA for traffic incident detection. One possibility is to apply unsupervised learning techniques to the collection of bottleneck distances produced for each vector in Algorithm 1. Another possibility is to use persistence landscapes rather than bottleneck distances.