This work is motivated by a challenging traffic data problem provided by the 2021 Algorithms for Threat Detection (ATD) program sponsored by the National Science Foundation (NSF) and the National Geospatial-Intelligence Agency (NGA) [1]. The complete dataset (i.e., without missing values) consists of hourly traffic flow (i.e., the number of vehicles) for 400 sensors (locations) in California from 2011 to 2020. A sample of the complete data is provided in Table 1.

The goal of the 2021 ATD challenge is to detect traffic flow anomalies at a sensor s at hour h on weekday w for an incomplete dataset where the majority of the data is not observed. A sample of the incomplete data is provided in Table 2. In fact, the ATD organizers downsampled each sensor by a particular sampling rate (i.e., one minus the missing rate): 1%, 2%, 5%, or 10%. The incomplete dataset has the same form as the complete dataset except that a large number of rows are missing as depicted in Table 2.

The definition of an anomaly is determined by the ATD program organizers. All analyses are at the granular level of sensor s, hour h, and weekday w, which is called a slice (see definition 2.1). The traffic flow data are detrended for each slice as the anomaly definition assumes stationarity. Specifically, a linear trend model is fit on the complete hourly data over time, and its predicted linear trend is then subtracted from the original data. The average traffic flow is added back to the resulting data. We refer to the dataset that results from the aforementioned steps as the complete detrended dataset D.

Let the detrended traffic flow be denoted by f. We refer to the complete detrended data (i.e., no missing values) concerning sensor s at hour h on weekday w as the “slice” which has the following definition.

Consequently, the true mean traffic flow for the slice is: (2.1) μs,h,w=1ns,h,w∑iDs,h,wi where Ds,h,wi is the ith traffic flow for sensor s at hour h on weekday w (i.e. multiple values collected from 2011 to 2020). Similarly, the true standard deviation for the slice is (2.2) σs,h,w=∑i(Ds,h,wi−μs,h,w)2ns,h,w The true anomaly label is determined by the ATD program organizers and is defined as the following: Is,h,wi=1(anomaly),|Ds,h,wi−μs,h,w|≥3σs,h,w0(normal),|Ds,h,wi−μs,h,w|<3σs,h,w.

In this project, the complete dataset D is not available. The sampling rate is provided for each sensor. As the definition of anomaly infers, the anomaly status is defined based on the complete dataset, but only the incomplete data are available to predict the anomaly. The purpose of this study is to use the incomplete data to accurately predict anomalies whose status is determined using the complete data.

Gaussian Process Regression (GPR) models are non-parametric kernel-based models based on the assumption that each observed response value yi is sampled from a multivariate normal distribution. Different from linear regressions, GP regression models are built on random processes with a Gaussian prior instead of a parametric formulation for the latent function. It also provides uncertainty measures over predictions In GPR models, the relationship between a latent function f(xi) and the target (response) yi can be written as follows: yi=f(xi)+ϵi,ϵi∼N(0,σ2), where the random noise ϵi is normally distributed with mean 0 and variance σ2 . The covariance function of the response in a GPR is cov(y)=K+σ2I where the entry at the i-th row and j-th column of K, Ki,j=k(xi,xj) . The covariance matrix characterizes correlations between different responses in the process. If xi and xj have similar kernel function values, then their GPR outputs, f(xi) and f(xj) , shall also be similar. For some new points X∗ , the prediction f(X∗) can be obtained from the joint distribution of f on the input X and f∗ on the new points X∗ , a simplified notation of f(X∗) as follows: ff∗∼Nm(X)m(X∗),KK∗K∗TK∗∗, where K=K(X,X) , K∗=K(X,X∗) , K∗∗=K(X∗,X∗) is a covariance matrix, whose entries are given by the covariance function (kernel), Ki,j=k(xi,xj) , and m(X)=m(X∗)=0 for convenience. The posterior predictive distribution is obtained by marginalizing out the training set latent variables [21]: p(f∗|y)=∫p(f,f∗|y)df=1p(y)∫p(y|f)p(f,f∗)df, and it can be written as p(f∗|y)=N(Kf∗,f(Kf,f+σ2I)−1y,Kf∗,f∗−Kf∗,f(Kf,f+σ2I)−1Kf,f∗). Notice that the covariance does not depend on the observed output y but only on the inputs X and X∗ .

Recently, GPs have become a popular method for non-parametric modeling, and there are several reasons for using GPs. First, the locations of the observations used to train the model do not necessarily correspond to the locations of the points for which we wish to investigate the function at [22]. Inference for the GP at a set of locations X∗ that differ from those corresponding to X can be achieved by evaluating the posterior distribution of f(X∗) . Second, there are various kernel functions that can be proposed to meet the modeling expectations. For example, if we believe the informativeness of past observations in explaining current data is a function of their time difference, the squared exponential function can be used in the covariance function, i.e., f(xi,xj)=h2exp[−(xi−xjλ)2] , where hyperparameters h and λ control the output scale of the function and the input, or time, scale. Kernels can also be combined together (e.g., addition and multiplication) to develop flexible nonlinear correlations.

Other non-parametric models, such as splines are also commonly used for spatiotemporal data analysis as they rely on fewer assumptions about the underlying data generating mechanisms. Although spline models can implement complex nonlinear functions, they are less efficient in modeling the effects of covariate interactions when compared to GPs.

Approximate Bayesian computation (ABC) has been used widely for parameter inference in ecology and biology as it does not require the specification of a likelihood function [8, 4]. The ABC approach estimates the posterior distributions of parameters by simulation data with parameters sampled from the prior distribution. A value for the tolerance (ϵ) is pre-specified to determine if the distance between a simulation dataset and the observed dataset is less than or equal to ϵ, it is retained. Otherwise, it is discarded.

GPs have been used as a replacement for supervised neural networks in nonlinear regression and extended for classification [26, 20]. Previous research has shown that these two approaches were equivalent for a neural network model with a single-layer fully-connected neural network in the limit of infinite width [17, 27, 28]. For such a neural network with independent and identically distributed (i.i.d.) random parameters, as a result of an affine transformation of the final hidden layer, each scalar output contains a sum of i.i.d. terms. According to the Central Limit Theorem, the neural network’s computed function is drawn from a generalized projection (GP) when the width is infinite [27]. Multilayer Perceptron (MLP) neural networks are popular in capturing nonlinear relationships. It has been shown that a MLP neural network and GP have similar behavior when the number of hidden neurons tends to infinity and weight decay is employed. Based on such correspondence, it is possible to provide exact Bayesian inference for neural networks with infinite width by evaluating the corresponding generalized partials. There are kernel functions designed to mimic multi-layer random neural networks, but they are only available outside of a Bayesian network and have not yet been identified as covariance functions for GP in previous work [27].

In this subsection, we describe how to pre-process the incomplete data which includes detrending the data in the same way as the ATD organizer, normalizing the data, and augmenting the data from similar groups: a) nearby hours of the same day and b) same hour of other weekdays. We first detrend the traffic flow following the same strategy described in Section 2.1. The resulting detrended data is denoted by d, and ds,h,w is the detrended total flow for slice (s,h,w) .

In the normalization step, instead of normalizing the data based on the slice of s, h, and w which has high sparsity, the day of week is aggregated into weekdays or weekends to increase the sample size for normalization. A weekend indicator wc is created to recode the day of week variable to a dichotomous variable with value 1 indicating weekends and 0 indicating weekdays, wc=1Dwi∈{Sat., Sun.}∩Doi=1,0Dwi∈{Mon., Tue., Wed., Thu., Fri.}∩Doi=1. Then we normalize the data for each slice of s (sensor), h (hour), and wc (weekdays or weekends). That is, data at the same locations and hours of weekdays are pooled together, and data at the same locations and hours of weekends are pooled together for normalization as follows, ds,h,wc′i=ds,h,wci−d¯s,h,wcsds,h,wc where d¯s,h,wc and sds,h,wc are the mean and standard deviation of slice (s,h,wc) . This approach is advantageous as compared to normalizing the data based on the slice of s, h, and w since there are only a few points available for each slice. It also works better than normalizing based on the slice of s and h because the traffic flows of weekdays and weekends are significantly different and thus not suitable for pooling.

In the data augmentation step, it is assumed that the traffic flow data is smooth so that observed data can be duplicated to other slices. We augment data from similar groups:

nearby hours of the same day.

Similarly, the observed data from the same hour on other weekdays are also transformed for data augmentation. For slice (s,h,w) , its traffic flow data are also used for other days of the week that have the same weekday category (weekday vs. weekend), i.e., slice (s,h,w′)s.t.wc(w′)=wc(w) . The transformation formula is described as follows: (4.2) ys,h,w′′(n+1)=g2(ds,h,w′i) (4.3) =ds,h,w′i−d′¯s,h,wsds,h,w′×sds,h,w′′+d′¯s,h,w′ where n is the current sample size of ys,h,w′′ and updated after each augmentation step, d′¯s,h,w and sds,h,w′ are the sample mean and standard deviation of the data for slice (s,h,w) , d′¯s,h,w′ and sds,h,w′′ are the sample mean and standard deviation of the data for slice (s,h,w′) .

When the augmentation process completes, ys,h,w′ contains the detrended data from the original incomplete dataset, transformed data from the neighboring hours of the same day, and transformed data from other similar days of week at the same hour and same week. ys,h,w′=ds,h,w′∪g1(ds,h±j,w′)∪g2(ds,h,w′′) where j=±1,…,±k and wc(w′)=wc(w) . The augmentation strategy can be extended to incorporate spatial patterns by using data from the closest stations on the same day and same hour. It is not used in our project as it does not improve the model’s performance.

Although the data augmentation step has increased the sample size to a certain degree, it is based on observed data and does not infer what the unobserved looks like. Temporal patterns could be further explored to get a full picture of the data. The purpose of GPR models is to recover the unobserved traffic flow from the augmented data ds,h,w′ . To reconstruct the complete traffic flow data, we apply the fitted GPR models for each slice to predict the unobserved traffic flow from the posterior predictive distribution. With the augmented data ys,h,w′ , we apply GPR to recover the complete traffic flow data from 2011 to 2020 based on its date information, denoted by ts,h,w . The dependent variable in the GPR is the traffic flow, and the dependent variable is the time index x, a new feature introduced to recover the time series pattern of the data. To find x, we first list all the possible dates that correspond to the hour and weekday value in slice (h,w) : Th,w={dates between 2011 and 2020:Th=h∪Tw=w}. Let Nh,w be the size of Th,w . Then the time index for i-th observation in slice (s,w,h) xs,h,wi for ys,h,w′ can be identified by matching its dates with T as follows, xs,h,wi={j:Th,wj=ts,h,wi}, where ts,h,wi is the date of i-th observation in slice (s,h,w) . We have a collection Xs,h,w⊂{1,…,Nh,w} by definition. Let ns,h,w be the size of Xs,h,w . For the set of observed data G={(xs,h,wi,ys,h,w′i),i=1,…,ns,h,w} , we use the GPR to model the relationship between the traffic flow and the time index on the augmented data. Assume that ys,h,w′=f(xs,h,w)+ϵ,ϵ∼N(0,σ2I). A GP prior is placed over f(xs,h,w) and simplifying the notation by removing the subscript (s,h,w) as follows: p(f(x)|x)=N(0,K) where K is the covariance matrix determined by the kernel function. In this project, the rational quadratic kernel is implemented, (4.4) kα,l(xi,xj)=1+d(xi,xj)22αℓ2−α, where α is the scale mixture parameter which determines the weighting of large or small-scale variations, l is the length scale of the kernel, and d(·,·) is the Euclidean distance. The details of the kernel selection are discussed in the following section.

There are several advantages of recovering the complete sequence. First, the true mean and standard deviation could be approximated from the recovered data and fed into the classification process. Second, with the reconstructed complete time series data, other features such as the 5-th, 25-th, 50-th, 75-th, and 95-th percentiles can be computed and used when building the final prediction model. These summary statistics might not be represented well by the original incomplete data which have less than five observations ( n<5 ) on average for each slice.

The choice of kernel profoundly affects the performance of a GP on a given task [23] as it controls properties of latent functions including smoothness and periodicity. In practice, sophisticated kernels are often achieved by composing a few standard kernel functions, which may also lead to overfitting. Commonly used kernels include the Gaussian radial basis function (RBF) kernel, rational quadratic kernel, expsine kernel, and Matérn kernel.

The RBF kernel is also known as the squared-exponential kernel and takes the form as follows: kl(xi,xj)=exp−d(xi,xj)22l2 where l determines the length scale of the associated hypothesis space of functions. A large value of γ shrinks the covariance between nearby sets of observations, which results in a more smooth output.

We use the Python Scikit-learn package [18] to optimize the hyperparameters in the kernel function, which are tuned in order to maximize the log-marginal-likelihood. Specifically, we use 0.1 and 1 as initial values for α and l, respectively. And the bounds are α,l∈(0.1,10) .

In Figure 2, we show the reconstructed time series for an arbitrary sensor (sensor 1) using GPR with a rational quadratic kernel. Sensor 1 has a 2% sampling rate. The blue curve represents the complete traffic flow and the orange is the reconstructed total flow. We only show the first two week’s estimation (2011-01-01 to 2011-01-14) for better resolution. Overall, the estimated total flow resembles the true values, but it does not fully recover the true distribution due to the sparsity of the observed data.

The incomplete traffic flow data has been detrended and normalized in the previous steps. There are additional features from the complete reconstructed data using GPR (i.e., estimated mean, standard deviation and percentiles). In this section, we will use d′ (i.e., the detrended and normalized downsampled data), as well as features from GPR to build a logistic regression model for each sampling rate. All variables are normalized using the estimated mean and standard deviation of its corresponding slice. For example, data for slice (s,h,w) , ds,h,w′i is normalized as follows: (4.5) zs,h,wi=ds,h,w′i−μˆs,h,wσˆs,h,w where i=1,…,|ds,h,w′| , and μˆs,h,w and σˆs,h,w are the mean and standard deviation of the reconstructed complete traffic flow for slice ({s,h,w) from GPR. Denote the normalized data of d′ by z. Features from the reconstructed data are also normalized in the same fashion. For example, the 5-th percentile w.r.t. the reconstructed data for slice (s,h,w) is normalized as follows: ps,h,w;5′=ps,h,w;5−μˆs,h,wσˆs,h,w

With the normalized data {z,p5′,p50′,p75′,p95′} , we train the logistic regression for each sampling rate. Define Nj be the size of sampling rate j where j∈{0.01,0.02,005,0.1} . Omitting the slice notation {s,h,w} , the logistic regression model for sampling rate j can be written as (4.6) logpi1−pi=β0j+β1jzi+β2jp5′i+β3jp25′i (4.7) +β4jp50′i+β5jp75′i+β6jp95′i where i=1,2,…,Nj .

The logistic regression model returns the probability that a traffic flow is anomalous, and a cutoff threshold cj needs to be specified to determine the anomaly status for each sampling rate j where j∈{0.01,0.02,005,0.1} . A grid search is used to choose the optimal value of cj . Note that the three standard deviation method to find anomalies is a special case for logistic regression when β0=β2=β3=β4=β5=β6=0 and cutoff becomes e3β11+e3β1 as zi=3 . So, the logistic regression model encompasses the three-standard deviation criterion. But it is more flexible than the three-standard deviation criterion because it allows the tuning process for cj to obtain better prediction results.

Once the optimal values are selected, an observation is classified as an anomaly if its predicted probability of being anomalous is greater than c. For example, the optimal value of c was found to be 0.31 for a sampling rate of 0.1 so that an observation with predicted anomaly probability ≥0.31 is classified as an anomaly. This is defined in the formula below. (4.8) I(zs,h,wi)=1(anomaly),pˆs,h,wi≥cj0(not anomaly),pˆs,h,wi<cj where cj is the optimal value for sampling rate j that sensor s corresponds to. Then we can compute the corresponding F1 score using Equation (5.1).

The classification results were evaluated with F1 score which is the metric adopted by the challenge organizer. It is defined as follows: (5.1) F1=true positivetrue positive+0.5(false positive + false negative). F1 score is widely used in many applications such as image classification and documentation classification. Unlike classification accuracy which is the proportion of the number of correct predictions from all predictions made, F1 score conveys the balance between precision and recall. For example, in the imbalanced setting when negative cases dominate, the classification accuracy would be high by simply labeling all observations negative. However, the F1 score would be 0 as the predictions for the negative observations are incorrect.

To test the performance of the proposed method, we randomly sample 200 for training and the remaining 200 sensors for testing. The proposed algorithm is run in parallel at the sampling rate level. Within each sampling rate run, it can be further paralleled at the slice level for GPR hyperparameter tuning, which is comparably the most computationally expensive part among all components. Overall, the proposed algorithm is efficient and takes about 1 hour for training and 7 minutes for testing. The test F1 score for each sampling rate is given in Table 4. We also compute the F1 score on the test dataset using the baseline algorithm provided by NSF, which labels an observation as an anomaly if the detrended traffic flow is above or below three standard deviations from the observed mean.

Although this paper focuses on the use of the Gaussian Process in the anomaly detection framework (see Figure 1), other methods may be used to reconstruct the complete total flow. One method applied is stochastic gradient boosting [11, 12]. Specifically, we use the XGboost library [7] to implement the XGBoost method. The GPR has an average F1 score of 0.6911 whereas the XGBoost has an average F1 score of 0.5776. The GPR has a descent F1 score for the sampling rate of 1% (0.5505) as compared to XGBoost (0.4074).

We also test the proposed framework by replacing logistic regression with XGBoost (column XGBoost-2 in Table 4). We observe that the XGBoost tends to overfit the training dataset and logistic regression gives better results for all four sampling rates in the test dataset.

We also applied the ABC algorithm to estimate the mean and standard deviation of each slice through the Python abcpy package [10]. Specifically, the parameter samples are generated from a uniform distribution ( μ∼(−0.5,0.5) , σ∼(0.5,1.5) ) and a Gaussian likelihood is assumed for the distribution of the total flow. The posterior mean and standard deviation are then used to determine if a traffic flow is anomalous based on the three-standard deviation rule. We examine the hyperparameter ϵ (the distance tolerance) at two different values, 8 and 10, after comparing the simulation datasets and the observed data. We find the results for ϵ=8 and ϵ=10 are very close; therefore only the first one is shown in Table 4. It is shown that the proposed algorithm has higher F1 scores than the ABC algorithm for all four sampling rates.

Furthermore, the F1 score on the test dataset using the baseline algorithm is also computed for comparison purposes. The baseline algorithm labels an observation as an anomaly if the detrended traffic flow is above or below three standard deviations from the observed mean. The proposed method has a higher F1 score for all four sampling rates than the baseline approach.

In this section, we show the anomaly probability on a map using the standardized latitude and longitude information to reveal the spatial pattern of anomalies (i.e., the original latitude and longitude values are not provided). Due to the nature of the data, we aggregate the anomalies by year and compute both the observed and predicted anomaly probability for each sensor. The left and right plots in Figure 3 are the true and predicted anomaly probability in 2016, respectively. In general, the predicted anomaly probability is close to the observed anomaly probability. In both plots, it is shown that sensors with large anomaly probabilities are mainly located in the northwestern area which is likely close to downtown Los Angeles. We are not able to show it on the exact map since the latitude and longitude provided by NSF are normalized. Few sensors located in the southeastern area also have relatively high anomaly probabilities. The scatter plot also showed that sensors that are close to each other had similar anomaly probability.

In this section we examine the hourly traffic patterns of holidays and non-holidays and their daily profiles.

Holidays vs. non-holidays. We analyze the hourly anomaly probability to investigate holiday and weekday patterns (see Figure 4). In general, holidays have a higher probability of experiencing unusual traffic than non-holidays. Moreover, holidays generally had anomalous traffic in the morning (4 AM to 10 AM), whereas non-holidays had higher anomaly probabilities from 10 AM to 8 PM. For holidays, people tend to leave home early, so this may explain the peak of the anomaly probability in the early morning [6]. For non-holidays, the anomaly probability patterns coincide with daily work commute patterns. Figure 4 indicated that the daily profile had two peaks for non-holidays, corresponding to two rush hour periods: one in the morning (11 AM to 12 PM) corresponding to lunch breaks and another one in the afternoon (7 PM to 8 PM) which corresponds to the end of a working day.

Daily profiles in holidays and non-holidays. Stratifying the traffic flow data by holidays and non-holidays, the hourly anomaly probabilities are further analyzed at the weekday level. Figure 5 showed the daily profiles by days of the week for holidays and non-holidays. We observe that there are large traffic flows on holidays that occurred on Mondays through Saturdays from 7 AM to 8 PM whereas Sunday has large traffic flows at 1 AM and a lower anomaly probability overall. We also found that Tuesday, Wednesday and Thursday have higher a probability of experiencing anomalous traffic than other days in the holiday group. Also, we found that 1 AM has the second highest probability of experiencing anomalies for most holidays which is consistent with a previous study [15].

For non-holidays, we find that Monday, Tuesday, Wednesday, and Thursday have similar patterns: a) two peaks at 11 AM to 12 PM and 7 PM to 8 PM and b) fewer anomalies at 8 AM. Friday’s pattern differes in that its anomaly probability gradually increases after 1 AM instead of dropping in the morning which is similar to Saturday and Sunday. It is interesting to point out that morning traffic volumes from Monday to Friday may be high but has a relatively smaller chance of reporting anomalies as compared to other time frames.