Contextual anomaly detection on time series: a case study of metro ridership analysis

Abstract

The increase in the amount of data collected in the transport domain can greatly benefit mobility studies and create high value-added mobility information for passengers, data analysts, and transport operators. This work concerns the detection of the impact of disturbances on a transport network. It aims, from smart card data analysis, to finely quantify the impacts of known disturbances on the transportation network usage and to reveal unexplained statistical anomalies that may be related to unknown disturbances. The mobility data studied take the form of a multivariate time series evolving in a dynamic environment with additional contextual attributes. The research mainly focuses on contextual anomaly detection using machine learning models. Our main goal is to build a robust anomaly score to highlight statistical anomalies (contextual extremums), considering the variability within the time series induced by the dynamic context. The robust anomaly score is built from normalized forecasting residuals. The normalization of the residuals is carried out using the estimated contextual variance. Indeed, there are complex dynamics on both the mean and the variance in the ridership time series induced by the flexible transportation schedule, the variability in transport demand, and contextual factors such as the station location and the calendar information. Therefore, they should be considered by the anomaly detection approach to obtain a reliable anomaly score. We investigate several prediction models (including an LSTM encoder–decoder of the recurrent neural network deep learning family) and several variance estimators obtained through dedicated models or extracted from prediction models. The proposed approaches are evaluated on synthetic data and real data from the smart card riderships of the Quebec Metro network. It includes a basis of events and disturbances that have impacted the transport network. The experiments show the relevance of variance normalization on prediction residuals to build a robust anomaly score under a dynamic context.

Appendix

1.1 Bias-variance Estimators

First, we propose two ways to learn and estimate the bias-variance based on the prediction residues produced by the forecasting models.

1. EMP: Empirical estimation on a prior sampling.

The estimation model is based on prior knowledge. We segment the contextual attribute space $\mathbf{c }$ into prior subspaces (subsamplings) defined by a set of constraints ($V^{inf}, V^{sup}$) given by expert knowledge. The bias ${\hat{B}}$ and variance ${\hat{\sigma }}$ estimators are summarized in three steps, as follows:

• 1. Extract from each prior $E_k$ the subsampling bias and variance.

• 2. Associate each time step t to its subsampling $E_k$.

• 3. Return the bias ${\hat{B}}_t$ and variance ${\hat{\sigma }}_t$ for each time step t.

$$\begin{aligned}&\{E_k : {t \in E_k\ |\ V_{k}^{inf}> \mathbf{c }(t) >V_{k}^{sup}} \}\\&{\hat{B}}_{E_k} = \sum _{t\in E_k}\frac{r_t}{\#E_k}={\hat{r}}_{E_k} \qquad \ {\hat{\sigma }}_{E_k} = \sqrt{\sum _{t\in E_k}\frac{(r_t-{\hat{B}}_{E_k})^2}{\#E_k}} \end{aligned}$$

2. ML: Machine learning-based estimation.

The estimation model can be learned by a machine learning algorithm. We train two prediction models to learn the bias and variance of the residues of the predictions from the contextual attributes.

The two models are similar in terms of estimating a type of mean (absolute for the bias and quadratic-centered for the variance) on a learned contextual subsample.

$$\begin{aligned}&Bias\ :\ \theta = {{\,\mathrm{argmin}\,}}_{\theta } \sum _{t} |M^{{\hat{B}}}_\theta (X_t)- r_t| {\hat{B}}(t) = M^{{\hat{B}}}_\theta (X_t) = {\hat{r}}_t \\&Variance\ :\ \theta = {{\,\mathrm{argmin}\,}}_{\theta }\sqrt{\sum _{t} |M^{{\hat{\sigma }}}_\theta (X_t) - (r_t-{\hat{B}}(t))^2|} \\&\quad {\hat{\sigma }}(t) = \sqrt{M_\theta ^{{\hat{\sigma }}}(X_t)} \end{aligned}$$

Second, we propose directly extracting an estimation of the bias and variance from a forecasting model. We propose exploring the extraction for a random forest and a deep neural network. Often, extracting the estimated bias from the model itself will lead to a result of zero since the model has been optimized to minimize this bias.

• RF: Random forest extraction

  In [35], the authors show that we often exploit valuable information about the distribution learned from a random forest by considering only the mean of the subsamples. From this assumption, we propose extracting the variance based on a learned subsampling of our random forest forecasting model.

  Let M be a random forest composed of $(T^1,..,T^n)$ binary trees. Each tree $T^k$ is composed of a set of leaves $L^k$. Values $j_i$ are assigned to each leaf during the learning phase according to their attribute modalities $X_i$. We define a tree walk operator $F^k(X_t)$ that takes attributes $X_t$ and returns for the associated leaf $L^k_i$, the set of assigned values.

  $$\begin{aligned} M(X_t) = \frac{1}{n} * \sum _{k \in [1,n]} \left( \sum _{j \in F^k(X_t)} \frac{j}{\#F^k(X_t)}\right) = {\hat{y}}_t \end{aligned}$$

  The prediction of an element by an RF model is similar to the weighted mean of a subsample formed by elements sharing a leaf. The weighting depends on the shared leaf number and shared element tree number. Shared leaf elements can be considered contextual neighbors on the basis of their attributes. Then, we can extract the bias (equal to 0) and variance from this contextual subsampling.

  $$\begin{aligned} {\hat{B}}(t) = 0\quad {\hat{\sigma }}(t) = \sqrt{ \frac{1}{n} * \sum _{k \in [1,n]}\left( \sum _{j\in F^k(X_t)} \frac{(j- {\hat{y}}_t)^2}{\#F^k(X_t)}\right) } \end{aligned}$$

• DEEP: Neural network extraction

  A second form of extraction is based on variational dropout [38], which aims to approximate Bayesian behavior in a deterministic network. A study in [12] applies this technique to an LSTM neural network to extract the confidence in the prediction model. Following the same line of research, we use the variational dropout to estimate the variance from our LSTM encoder-predictor model.

  Let $M_{{\hat{\theta }}}$ be a neural network that infers $y_t$ from $X_t$.

  $$\begin{aligned} \theta ={{\,\mathrm{argmin}\,}}_{\theta } \sum |M^{\theta }(X)-y|^2 \quad M_{\theta }(x_t) = \hat{y_t} \end{aligned}$$

  The neural network aims to capture the link between the attributes and prediction targets through an embedding of the attribute space into the prediction space. Successive nonlinear projections in the abstract space Z are used to this end. These abstract spaces give us abstract representations $z_t$ of our elements that capture the topological structure of our data. We can exploit such spaces to perform contextual subsampling by defining a neighborhood in Z space. The contextual subsampling will be based on the contextual information captured by M. The main issue comes from the definition of a neighborhood ${\mathcal {B}}(z_t)$ in Z space.

  $$\begin{aligned}&\{ {\mathcal {B}}(z_t) : k\ tq\ z_k \in [z_t \pm \varepsilon ]\}\ with\ z,\varepsilon \in {\mathcal {R}}^{\#Z} \\&\quad {\hat{B}}(t) = \sum _{k \in {\mathcal {B}}(z_t)}\frac{|\hat{y_k}-y_t|}{\# {\mathcal {B}}(z_t)}={\hat{r}}_t\\&\quad \ {\hat{\sigma }}(t) = \sqrt{\sum _{k \in {\mathcal {B}}(z_t)}\frac{((\hat{y_k}-y_t)-{\hat{B}}(t))^2}{\# {\mathcal {B}}(z_t)}} \end{aligned}$$

  This issue can be avoided with a variational neural network $M_{\theta }^{var}$ based on an explicit (variational layer) or implicit (variational dropout) random drawing by generating a virtual sampling that self-defines the neighborhood in Z space.

  $$\begin{aligned} \theta = {{\,\mathrm{argmin}\,}}_{\theta } \sum | (M_{\theta }^{var}(X)-y)|^2 \quad \sum _{m}\frac{M_{\theta }^{var}(x_t)}{m}= \hat{y_t} \end{aligned}$$

  The stochastic projections of model $M_{\theta }^{var}$ transform the latent representations $z_t$ into a collection of probabilistic points. We can access the probabilistic clouds of predictions for an element by making many predictions. This gives us a virtual contextual subsampling from which we can estimate the mean and variance.

  $$\begin{aligned} \hat{{\mathcal {B}}}=0 \quad {\hat{\sigma }}(t) =\sqrt{\sum _{m}\frac{(M_{\theta }^{var}(x_t)-{\hat{y}}_t)^2}{m}} \end{aligned}$$

1.2 Forecasting models

1.2.1 Encoding cyclical features

Cyclical encoding aims to encode continuous cyclic attributes by preserving their cyclic structure. Instead of having a large one-hot vector per feature, the sine and cosine encodings project each attribute on a two-dimensional plane. However, some contextual information contains more than one cyclical structure. Using several pairs of sines and cosines with different frequencies can allow us to better express meaningful and compact structures. For instance, we can express several pieces of periodical information (weekly, monthly, and seasonal) by encoding the position of the day in the year with several pairs of sines and cosines with well-chosen frequencies, i.e., 1/53 for a weekly structure, 1/31 for monthly, and 1/1 for a yearly structure.

This technique yields compact and meaningful attributes, in contrast to the bulky and sparse hot encoding.

1.2.2 Random forest training

The forecasting and bias-variance estimation models are optimized through a mean-square error (MSE) optimization loss. We use a standard scikit-learn random forest regressor [40]. The random forest parameters are tuned through a random search using cross-validation combined with early stopping to control the number, size and depth of trees to avoid overfitting.

1.2.3 LSTM EP : architecture and training

In our previous research [11], we proposed an LSTM encoder predictor for ridership forecasting by using both long-term and short-term attributes.

Fig. 11

LSTM-EP architecture with the layer size for the real data

The model was designed to manage the structural variability in the data induced by the transport plan. A simplified version of the model (Fig. 11) is applied in the current work thanks to the regular structure of the data.

First, the long-term features are synthesized through a multilayer perceptron neural network. Then, a pair of encoder–predictor LSTM layers attempt to capture the contextual influence and infer the short-term dynamics of the multivariate time series. Finally, another multilayer perceptron attempts to interpret the prediction embedding $Z^p$ to produce a prediction ${\hat{y}}$. This model takes as input the contextual attribute $X_t$ and past horizon value $y^p_t=[y_{t-p},\ldots ,y_t)$ and aims to forecast a future horizon $[y_{t},y_{t+f})]$. Such a model reconstructs the time step t and then infers the temporal evolution on a future horizon $[t+1,t+f]$. Dropout layers are placed in almost every layer to avoid overfitting and to allow variational dropout. The size is manually chosen through a compromise between three components of size, performance and overfitting.

The encoder predictor model is implemented based on the TensorFlow [41] environment with Keras [42] as a library and a high-level neural network API. Training is performed through 3 training loops with gradient reduction and early stopping. We use an adaptive gradient (ADAM), and we reduce the batch size between each loop. The first training loop is a type of initialization in which we keep only the reconstruction task in the learning loss. Then, we add multistep forecasting with a higher weight to the t+1 prediction loss.