Statistical Anomaly Detection via Composite Hypothesis Testing for Markov Models

TL;DR

This paper introduces a new threshold estimator for the Hoeffding test under Markov models, improving anomaly detection accuracy in communication and transportation networks by better controlling false alarms.

Contribution

It develops a novel estimator based on a CLT for the empirical measure, providing theoretical guarantees and demonstrating improved performance over existing methods.

Findings

Our estimator better controls false alarms.

It maintains high detection probabilities.

Effective in cyber security and transportation applications.

Abstract

Under Markovian assumptions, we leverage a Central Limit Theorem (CLT) for the empirical measure in the test statistic of the composite hypothesis Hoeffding test so as to establish weak convergence results for the test statistic, and, thereby, derive a new estimator for the threshold needed by the test. We first show the advantages of our estimator over an existing estimator by conducting extensive numerical experiments. We find that our estimator controls better for false alarms while maintaining satisfactory detection probabilities. We then apply the Hoeffding test with our threshold estimator to detecting anomalies in two distinct applications domains: one in communication networks and the other in transportation networks. The former application seeks to enhance cyber security and the latter aims at building smarter transportation systems in cities.

Figures

Tables

Table 1. TABLE I : ROC points vs. target FPR ( N = 4 𝑁 4 N=4 , n = 50 𝑛 50 n=50 ).

target FPR $\beta$	HTWC-1		HTWC-2		HTSV
	FPR	TPR	FPR	TPR	FPR	TPR
0.001	0.002	0.885	0.0	0.816	0.402	0.999
0.01	0.011	0.965	0.002	0.888	0.752	1.0
0.02	0.018	0.983	0.003	0.943	0.844	1.0
0.03	0.025	0.99	0.01	0.96	0.898	1.0
0.04	0.038	0.99	0.018	0.971	0.927	1.0
0.05	0.047	0.991	0.029	0.981	0.945	1.0

Table 2. TABLE II : ROC points vs. target FPR ( N = 6 𝑁 6 N=6 , n = 100 𝑛 100 n=100 ).

target FPR $\beta$	HTWC-1		HTWC-2		HTSV
	FPR	TPR	FPR	TPR	FPR	TPR
0.001	0.001	1.0	0.0	1.0	0.997	1.0
0.01	0.008	1.0	0.003	1.0	1.0	1.0
0.02	0.017	1.0	0.005	1.0	1.0	1.0
0.03	0.028	1.0	0.017	1.0	1.0	1.0
0.04	0.037	1.0	0.017	1.0	1.0	1.0
0.05	0.055	1.0	0.019	1.0	1.0	1.0

Table 3. TABLE III : Key features of the detected abnormal jams.

index	start time	detected time	latitude	longitude	delay (in seconds)	length (in meters)	alert type
788	12:25:0.302	12:30:0.0	42.361951	-71.117963	232.0	3568.0	heavy traffic
1502	15:35:0.072	15:40:0.0	42.356275	-71.119852	585.0	844.0	heavy traffic
2412	19:25:0.365	19:30:0.0	42.342549	-71.085011	643.0	3568.0	heavy traffic
3005	21:25:0.238	21:30:0.0	42.349125	-71.10778	168.0	1962.0	weather hazard
3094	21:35:0.267	21:40:0.0	42.373336	-71.097731	509.0	897.0	road closed
3126	21:35:0.326	21:40:0.0	42.355048	-71.110335	528.0	1293.0	heavy traffic

Equations

$${\Gamma }_n({\theta }_{ij})=\frac{1}{n}\sum _{l=1}^n\text{mathbbm}1\{Z_l={\theta }_{ij}\},$$

$$\mathbf{\pi }=({\pi }_{ij};i,j=1,\dots ,N)=({\tilde{\pi }}_k;k=1,\dots ,N^2),$$

$$p({\theta }_{ij}|{\theta }_{kl})=\text{mathbbm}1\{i=l\}{q}_{ij},k,l,i,j=1,\dots ,N,$$

$$D({\left.\bm{\Gamma}_{n}\right\|{\bm{\pi}}})=\sum_{i=1}^{N}\sum_{j=1}^{N}\Gamma_{n}(\theta_{ij})\log\frac{\Gamma_{n}(\theta_{ij})/\big{(}\sum\nolimits_{t=1}^{N}\Gamma_{n}(\theta_{it})\big{)}}{\pi_{ij}/\big{(}\sum\nolimits_{t=1}^{N}\pi_{it}\big{)}},$$

$${\mathbf{U}}_n=\sqrt{n}({\mathbf{\Gamma }}_n-\mathbf{\pi }),$$

$${\mathbf{\Gamma }}_n=({\Gamma }_n({\theta }_{11}),\dots ,{\Gamma }_n({\theta }_{1N}),\dots ,{\Gamma }_n({\theta }_{N1}),\dots ,{\Gamma }_n({\theta }_{NN})).$$

$${\mathcal{H}}_n\text{rejects}\mathcal{H}\text{if and only if}D({\mathbf{\Gamma }}_n\Vert \mathbf{\pi })>\eta ,$$

$$\beta =\text{mathbbm}{P}_{\mathcal{H}}(D({\mathbf{\Gamma }}_n\Vert \mathbf{\pi })>\eta ),$$

$${\eta }_{n,\beta }^{\text{sv}}\approx -(1/n)\log (\beta ),$$

$$\frac{{\pi }_{ij}}{{\sum }_{t=1}^N{\pi }_{it}}=\frac{{\pi }_{ij}}{{\sum }_{t=1}^N{\pi }_{ti}}=q_{ij},i,j=1,\dots ,N.$$

$${\mathbf{Q}}=\left[{\begin{array}[]{*{20}{c}}{0.1}&{0.2}&{0.7}\\ 0&{0.2}&{0.8}\\ {0.6}&{0.15}&{0.25}\end{array}}\right],$$

$${\mathbf{P}}=\left[{\begin{array}[]{*{20}{c}}{0.1}&{0.2}&{0.7}&0&0&0&0&0&0\\ 0&0&0&0&{0.2}&{0.8}&0&0&0\\ 0&0&0&0&0&0&{0.6}&{0.15}&{0.25}\\ {0.1}&{0.2}&{0.7}&0&0&0&0&0&0\\ 0&0&0&0&{0.2}&{0.8}&0&0&0\\ 0&0&0&0&0&0&{0.6}&{0.15}&{0.25}\\ {0.1}&{0.2}&{0.7}&0&0&0&0&0&0\\ 0&0&0&0&{0.2}&{0.8}&0&0&0\\ 0&0&0&0&0&0&{0.6}&{0.15}&{0.25}\end{array}}\right],$$

$$U_{n,k}=\sqrt{n}({\Gamma }_n({\tilde{\theta }}_k)-{\tilde{\pi }}_k),k=1,\dots ,N^2.$$

$$f_k(Z)=\text{mathbbm}1\{Z={\tilde{\theta }}_k\}.$$

$${\mathbf{U}}_n\underset{n\to \infty }{\overset{\text{d}}{\to }}\mathcal{N}(0,\mathbf{\Lambda }),$$

$${\Lambda }_{ij}={\tilde{\pi }}_i({\mathbf{I}}_{ij}-{\tilde{\pi }}_j)+\sum _{m=1}^{\infty }[{\tilde{\pi }}_i({\mathbf{P}}_{ij}^m-{\tilde{\pi }}_j)+{\tilde{\pi }}_j({\mathbf{P}}_{ji}^m-{\tilde{\pi }}_i)],$$

$$h(\mathbf{\nu })=D(\mathbf{\nu }\Vert \mathbf{\pi })=\sum _{i=1}^N\sum _{j=1}^N{\nu }_{ij}\log \frac{\frac{{\nu }_{ij}}{{\sum }_{t=1}^N{\nu }_{it}}}{\frac{{\pi }_{ij}}{{\sum }_{t=1}^N{\pi }_{it}}}.$$

$$D({\mathbf{\Gamma }}_n\Vert \mathbf{\pi })$$

$$\frac{{\partial h({\bm{\nu}})}}{{\partial{\nu_{ij}}}}=\log{\nu_{ij}}-\log\bigg{(}{\sum\limits_{t=1}^{N}{{\nu_{it}}}}\bigg{)}-\log{\pi_{ij}}+\log\bigg{(}{\sum_{t=1}^{N}{{\pi_{it}}}}\bigg{)},$$

$$\nabla h(\mathbf{\pi })=0.$$

$$\frac{{\partial }^{2}h(\mathbf{\nu })}{\partial {\nu }_{ij}\partial {\nu }_{kl}}=\{\begin{array}{c}0,\text{ if }k\ne i,\text{hfill}\\ \frac{1}{{\nu }_{ij}}-\frac{1}{{\sum }_{t=1}^N{\nu }_{it}},\text{ if }k=i\text{ and }l=j,\text{hfill}\\ -\frac{1}{{\sum }_{t=1}^N{\nu }_{it}},\text{ if }k=i\text{ and }l\ne j.\text{hfill}\end{array}$$

$$2nD\left({\mathbf{\Gamma }}_n\Vert \mathbf{\pi }\right)$$

$$=n{\left({\mathbf{\Gamma }}_n-\mathbf{\pi }\right)}^{\prime }{\nabla }^{2}h({\tilde{\mathbf{\Gamma }}}_n)\left({\mathbf{\Gamma }}_n-\mathbf{\pi }\right)$$

$$={\mathbf{U}}_n^{\prime }{\nabla }^{2}h({\tilde{\mathbf{\Gamma }}}_n){\mathbf{U}}_n,$$

$${\nabla }^{2}h({\tilde{\mathbf{\Gamma }}}_n)\underset{n\to \infty }{\overset{\text{w.p.1}}{\to }}{\nabla }^{2}h(\mathbf{\pi }).$$

$$D({\mathbf{\Gamma }}_n\Vert \mathbf{\pi })=\frac{1}{2n}{\mathbf{U}}_n^{\prime }{\nabla }^{2}h({\tilde{\mathbf{\Gamma }}}_n){\mathbf{U}}_n\underset{n\to \infty }{\overset{\text{d}}{\to }}\frac{1}{2n}{\mathbf{U}}^{\prime }{\nabla }^{2}h(\mathbf{\pi })\mathbf{U}.$$

$$D({\mathbf{\Gamma }}_n\Vert \mathbf{\pi })=\frac{1}{2n}{\mathbf{U}}_n^{\prime }{\nabla }^{2}h({\tilde{\mathbf{\Gamma }}}_n){\mathbf{U}}_n\underset{n\to \infty }{\overset{\text{d}}{\to }}\frac{1}{2n}\sum _{k=1}^{N^2}{\rho }_k{\chi }_{1k}^2,$$

$${\eta }_{n,\beta }^{\text{wc}}\approx F_{\text{em}}^{-1}(1-\beta ;n),$$

$${{\hat{\tilde{\pi}}}_{k}}=\max\bigg{\{}{\frac{1}{n_{0}}{\sum\limits_{i=1}^{{n_{0}}}{\mathbbm{1}{\big{\{}{Z_{i}^{\scriptscriptstyle{(0)}}={{{\tilde{\theta}}_{k}}}}\big{\}}}}},~{}\varepsilon}\bigg{\}},\;k=1,\ldots,N^{2},$$

$$\hat{p}({\theta }_{ij}|{\theta }_{kl})=\text{mathbbm}1\{i=l\}{\hat{q}}_{ij},k,l,i,j=1,\dots ,N,$$

$${{\hat{\Lambda}}_{ij}}={{\hat{\tilde{\pi}}}_{i}}({{{\mathbf{I}}_{ij}}-{{\hat{\tilde{\pi}}}_{j}}})+\sum\limits_{m=1}^{{m_{0}}}\Big{[}{{\hat{\tilde{\pi}}}_{i}}({\hat{\mathbf{P}}_{ij}^{m}-{{\hat{\tilde{\pi}}}_{j}}})+{{\hat{\tilde{\pi}}}_{j}}({\hat{\mathbf{P}}_{ji}^{m}-{{\hat{\tilde{\pi}}}_{i}}})\Big{]},$$

Full text

I Introduction

For a given system, Statistical Anomaly Detection (SAD) involves learning from data the normal behavior of the system and identifying/reporting time instances corresponding to atypical system behavior. SAD has vast applications. For instance, motivated by the importance of enhancing cyber security, recent literature has seen applications in communication networks; see, e.g., [1, 2, 3, 4]. The behavior of the system is typically represented as a time series of real vectors and, in its most general version, anomaly detection is done through some Composite Hypothesis Test (CHT).

Specifically, a CHT aims to test the hypothesis that a given sequence of observations is drawn from a known Probability Law (PL) (i.e., probability distribution) defined on a finite alphabet [5]. Among numerous such tests, the one proposed by Hoeffding [6] has been well known for decades. When implementing the Hoeffding test in the context of SAD, one must set appropriately a threshold $\eta$ so as to ensure a low false alarm rate while maintaining a reasonably high detection rate. In the existing literature, this threshold is typically estimated by using Sanov’s theorem [7] – a large deviations result. Note that such an estimator (let us denote it by $\eta^{\text{sv}}$) is valid only in the asymptotic sense. In practice, however, only a finite number of observations are available, and it can be observed in simulations that $\eta^{\text{sv}}$ is not accurate enough, especially for relatively small sample sizes.

Our contributions in this paper include:

Under Markovian assumptions, we leverage a Central Limit Theorem (CLT) for a selected empirical measure related to the test statistic of the Hoeffding test, so as to establish weak convergence results for the test statistic, and derive a threshold estimator $\eta^{\text{wc}}$ therefrom, thus, extending the work of [5] which tackles the problem under independent and identically distributed (i.i.d.) assumptions. 2. 2.

We propose algorithms to calculate the threshold estimator $\eta^{\text{wc}}$ obtained above for the ordinary and a robust version of the Hoeffding test, respectively. We assess the advantages of our estimator over earlier work through numerical experiments. 3. 3.

We apply the Hoeffding test with our threshold estimator to two types of systems for the purpose of anomaly detection: $(i)$ a communication network with flow data simulated by the software package SADIT [8]; and $(ii)$ a real transportation network with jam data reported by Waze, a smartphone GPS navigation application. To the best of our knowledge, the latter is a novel application of anomaly detection.

A preliminary conference version of this work appeared in [9]. The present paper includes detailed technical arguments, derives results for the robust version of the Hoeffding test, expands the numerical comparisons with earlier work, and develops the traffic jam anomaly detection application.

The rest of this paper is organized as follows. In Sec. II we review related work. We formulate the threshold estimation problem in Sec. III and derive theoretical results in Sec. IV. Sec. V contains experimental results. Concluding remarks are in Sec. VI and a number of proofs appear in the Appendix.

Notational conventions: All vectors are column vectors. For economy of space, we write ${\mathbf{x}}=(x_{1},\ldots,x_{\text{dim}({\mathbf{x}})})$ to denote the column vector ${\mathbf{x}}$, where $\text{dim}({\mathbf{x}})$ is its dimension. We use prime to denote the transpose of a matrix or vector. Denote by $\mathbb{N_{+}}$ the set of all nonnegative integers. $\|{\mathbf{x}}\|$ denotes the $\ell_{2}$-norm of a vector ${\mathbf{x}}$, $\lfloor x\rfloor$ the integer part of a positive number $x$, $|\mathcal{A}|$ the cardinality of a set $\mathcal{A}$, $\log$ the natural logarithm, $\mathbbm{P}(A)$ the probability of an event $A$, $\mathbbm{E}[X]$ the expectation of a random variable $X$, and $\text{Cov}(X_{1},X_{2})$ the covariance between two random variables $X_{1}$ and $X_{2}$. We use $\mathcal{N}({\mathbf{0}},{\bm{\Sigma}})$ to denote a Gaussian distribution with zero mean and covariance matrix ${\bm{\Sigma}}$. $X_{1}\simeq X_{2}$ indicates that the two random variables $X_{1}$ and $X_{2}$ have approximately the same distribution. $\mathbbm{1}\{\cdot\}$ denotes the indicator function and $\xrightarrow[n\to\infty]{\textit{d}}$ (resp., $\xrightarrow[n\to\infty]{{\textit{w.p.1}}}$) denotes convergence in distribution (resp., with probability one) as $n$ approaches infinity.

II Related Work

Modeling network traffic as stationary in time, [1] applies two methods: one assumes the traffic to be an i.i.d. sequence and the other assumes observations of system activity follow a finite-state Markov chain. Both methods are extended in [4] to the case where system activity is time-varying. When implementing the Hoeffding test, however, both [1] and [4] use the large deviations estimator $\eta^{\text{sv}}$ to calculate the detection threshold in a finite sample-size setting, thus controlling the false alarm rate not well enough.

To derive a more accurate threshold estimator, [5, 10] use a procedure commonly used by statisticians: deriving results based on Weak Convergence (WC) of the test statistic in order to approximate the error probabilities of the Hoeffding test. Under i.i.d. assumptions, [5] (see also [10, 11]) proposes an alternative estimator for $\eta$ (let us denote it by $\eta^{\text{wc}}$), which is typically more accurate than $\eta^{\text{sv}}$, especially when not that many samples are available.

There has also been work on obtaining a tighter approximation of $\eta$ by refining Sanov’s theorem [12]. However, such refinements of large deviation results are typically faced with computational difficulty; for instance, as noted in [10], using the results of [12] requires the computation of a surface integral.

Several alternative anomaly detection approaches have been proposed, using for instance change detection methods [13]. We refer the reader for a comprehensive review of alternative methods to [13] and [1].

III Problem Formulation

To model the statistical properties of a general system, we introduce a few more notational conventions and some definitions. Let $\bm{\Xi}=\{{\xi_{i}};~{}i=1,\ldots,N\}$ be a finite alphabet containing $N$ symbols $\xi_{1},\ldots,\xi_{N}$, and ${\mathbf{Y}}=\{{Y_{l}};~{}l=0,1,2,\ldots\}$ a time series of observations. Define the null hypothesis $\cal H$ as: ${\mathbf{Y}}$ is drawn according to a Markov chain with state set $\bm{\Xi}$ and transition matrix ${\mathbf{Q}}=\left[{{q_{ij}}}\right]_{i,\,j=1}^{N}$. To further characterize the stochastic properties of ${\mathbf{Y}}$, we define the empirical Probability Law (PL) by

[TABLE]

where ${Z_{l}}=({{Y_{l-1}},~{}{Y_{l}}}),~{}l=1,\ldots,n$, ${\theta_{ij}}=({{\xi_{i}},~{}{\xi_{j}}})\in\bm{\Xi}\times\bm{\Xi}$, $i,\,j=1,\ldots,N$. Denote the transformed alphabet ${\bm{\Theta}}=\{{{\theta_{ij}}};~{}i,\,j=1,\ldots,N\}=\{{{{\tilde{\theta}}_{k}}};~{}k=1,\ldots,N^{2}\}$ and note ${\bm{\Theta}}=\bm{\Xi}\times\bm{\Xi}$ with ${{\tilde{\theta}}_{1}}={\theta_{11}},\ldots,{{\tilde{\theta}}_{N}}={\theta_{1N}},\ldots,{{\tilde{\theta}}_{({N-1})N+1}}={\theta_{N1}},\ldots,{{\tilde{\theta}}_{{N^{2}}}}={\theta_{NN}}$. Let also the set of PLs on ${\bm{\Theta}}$ be $\mathcal{P}({\bm{\Theta}})$.

The transformed observations ${\mathbf{Z}}=\{{{Z_{l}};~{}l=1,2,\ldots}\}$ form a Markov chain evolving on ${\bm{\Theta}}$; denote its transition matrix by ${\mathbf{P}}=[{{p_{ij}}}]_{i,\,j=1}^{{N^{2}}}$ and the stationary distribution by

[TABLE]

where ${{\pi_{ij}}}$ denotes the probability of seeing $\theta_{ij}$, and ${{\tilde{\pi}}_{1}}={{\pi}_{11}},\ldots,{{\tilde{\pi}}_{N}}={{\pi}_{1N}},\ldots,{{\tilde{\pi}}_{\left({N-1}\right)N+1}}={{\pi}_{N1}},\ldots,{{\tilde{\pi}}_{{N^{2}}}}={{\pi}_{NN}}$. We have [7]

[TABLE]

which enables us to obtain ${\mathbf{P}}$ directly from ${\mathbf{Q}}$; see Remark 2 for an example. We can now restate the null hypothesis $\mathcal{H}$ as: the Markov chain ${\mathbf{Z}}=\{{{Z_{l}};~{}l=1,2,\ldots}\}$ is drawn from PL ${\bm{\pi}}$.

To quantify the distance between the empirical PL $\bm{\Gamma}_{n}$ and the actual PL ${{\bm{\pi}}}$, one considers the relative entropy (or divergence) between ${\bm{\Gamma}_{n}}$ and ${{{\bm{\pi}}}}$:

[TABLE]

and the empirical measure:

[TABLE]

where ${\bm{\pi}}$ is defined in (2) and $\bm{\Gamma}_{n}$ is the vector

[TABLE]

Let now ${{\cal{H}}_{n}}$ be the output of a test that decides to accept or to reject the null hypothesis $\cal H$ based on the first $n$ observations in the sequence ${\mathbf{Z}}$. Under Markovian assumptions (Assumption 1 in Sec. IV), the Hoeffding test [7] is given by

[TABLE]

where $D({\left.{{\bm{\Gamma}_{n}}}\right\|{{{\bm{\pi}}}}})$ (cf. (4)) is the test statistic and $\eta$ is a threshold.

It is known that the Hoeffding test (6) satisfies asymptotic Newman-Pearson optimality [1, 4], in the sense that it maximizes the exponential decay rate of the misdetection probability over all tests with a false positive probability with exponential decay rate larger than $\eta$. Thus, an appropriate threshold $\eta$ should enable the test to have a small false positive rate while maintaining a satisfactorily high detection rate.

The theoretical false positive rate [5] of the test (6) is given by

[TABLE]

where the subscript $\cal H$ indicates that the probability is taken under the null hypothesis.

Given a tolerable (target) $\beta$, by conducting an ROC analysis for the Hoeffding test using labeled training data, we could “tune” $\eta$ such that the corresponding discrete test111A discrete test corresponds to a fixed specific value for $\eta$ in (6). [14] has a small false alarm rate and a high detection rate. In particular, we could select an $\eta$ corresponding to a point close to the northwest corner of the ROC graph. However, such tuning is too expensive and depends heavily on the quality and quantity of the training data. We can also, in principle, obtain the corresponding $\eta$ in (7) by directly simulating the samples of the test statistic $D({\left.{{\bm{\Gamma}_{n}}}\right\|{{{\bm{\pi}}}}})$, thus deriving an empirical Cumulative Distribution Function (CDF) and using its $(1-\beta)$-quantile. However, we will note in Remark 5 that this is also computationally too expensive when applied through a so-called “windowing” technique for purposes of anomaly detection. Thus, we seek to estimate $\eta$ without directly simulating the statistic. To that end, existing work involves using Sanov’s theorem [7] to derive an estimator for $\eta$. Specifically, for large enough $n$, by replacing the right hand side in (7) by an exponential we can obtain a minimal $\eta$ that suffices to bring the false positive rate below $\beta$ [1, 4]. Such an $\eta$ is given by

[TABLE]

where we use the $n,\beta$ subscript to denote the dependence of this estimator on $\beta$ and $n$ and the label sv indicates that it is obtained from Sanov’s theorem. We note that the estimator (8) does not contain any direct distributional information of the statistic $D({\left.{{\bm{\Gamma}_{n}}}\right\|{{{\bm{\pi}}}}})$; this might be one of the causes leading to inaccurate estimation of $\eta_{n,\beta}$, especially when the sample size $n$ is relatively small in practice. To see this more clearly, one can consider an extreme scenario where $N=4$, $\beta=10^{-1000}$, and $n=50$ (this is a reasonably small value; comparable to $N^{2}=16$). Then by (8), $\eta_{n,\beta}^{\text{sv}}$ would be way bigger than necessary, tending to yield a test with zero false alarm rate but also zero detection rate for a typical test set. The issue arises because we use an asymptotic large deviations result for a relatively modest value of $n$. Our primary goal in this paper is to derive an alternative threshold estimator, which would hopefully be more accurate than $\eta_{n,\beta}^{\text{sv}}$ for modest values of $n$, in terms of a certain metric that we will introduce in Sec. V.

IV Theoretical Results

We introduce the following assumption.

Assumption 1

${\mathbf{Z}}=\{{{Z_{l}};~{}l=1,2,\ldots}\}$* is an aperiodic, irreducible, and positive recurrent Markov chain ([15]) evolving on ${\bm{\Theta}}$ with transition matrix ${\mathbf{P}}$, stationary distribution ${{\bm{\pi}}}$, and with the same ${{\bm{\pi}}}$ as its initial distribution.*

Remark 1

Since ${\bm{\Theta}}$ is a finite set, under Assumption 1 ${\mathbf{Z}}$ is uniformly ergodic [15]. Assuming ${{\bm{\pi}}}$ as the initial distribution is done for notational simplicity; our results apply for any feasible initial distribution. Note also that, under Assumption 1, ${{\bm{\pi}}}$ must have full support over ${\bm{\Theta}}$; i.e., each entry in ${{\bm{\pi}}}$ is strictly positive.**

Lemma 1

Suppose Assumption 1 holds. Then

[TABLE]

Proof:

See Appendix A. ∎

Remark 2

Under Assumption 1, Remark 1 and Lemma 1 imply that all entries of ${\mathbf{Q}}$ are strictly positive, indicating that any two states of the original chain ${\mathbf{Y}}$ are connected. This is a stringent condition; yet, in practice, if some $\pi_{ij}$ in (9) is zero, we can replace it with a small $\varepsilon>0$, and then normalize the modified vector ${{\bm{\pi}}}$, thus ensuring that Assumption 1 is satisfied.**

Another reason why we set the zero entries in ${{\bm{\pi}}}$ to $\varepsilon>0$ is for convenience of computing the original transition matrix ${\mathbf{Q}}$, hence ${\mathbf{P}}$, via (9) and (3). If we simply eliminate the corresponding states in ${\mathbf{Z}}$, then it is possible that the number of the remaining states is not the square of some integer $N$; this would prevent us from easily recovering ${\mathbf{P}}$ from ${{\bm{\pi}}}$. Consider the following example: Assuming

[TABLE]

then by (3) we have

[TABLE]

and, by direct calculation, we obtain ${\bm{\pi}}=(0.03,0.07,0.23,0,\linebreak[3]0.05,0.14,0.3,0.07,0.11)$. Note that only $8$ entries in ${\bm{\pi}}$ are non-zero and $8$ is not the square of some integer $N$. Thus, if we eliminate the state corresponding to the zero entry in ${\bm{\pi}}$, it will be hard to recover ${\mathbf{Q}}$, hence ${\mathbf{P}}$.**

IV-A Weak Convergence of Empirical Measure

Let us first establish CLT results for one-dimensional empirical measures

[TABLE]

For $k\in\{1,\ldots,N^{2}\}$ define

[TABLE]

Lemma 2

Suppose Assumption 1 holds. Then a Central Limit Theorem (CLT) holds for $U_{n,k}$; that is, $U_{n,k}\xrightarrow[n\to\infty]{{\textit{d}}}\mathcal{N}({0,~{}\sigma_{k}^{2}})$ with $\sigma_{k}^{2}=\operatorname{Cov}({{f_{k}}({{Z_{1}}}),~{}{f_{k}}({{Z_{1}}})})+2\sum_{m=1}^{\infty}{\operatorname{Cov}}({{f_{k}}({{Z_{1}}}),~{}{f_{k}}({{Z_{1+m}}})})<\infty$.

Proof:

See Appendix B. ∎

Now we state the CLT [16, Thm. 3.1] for the multidimensional empirical measure ${\mathbf{U}}_{n}=(U_{n,k};\,k=1,\ldots,N^{2})$ as Lemma 3. Several different proofs for this result are available in [16] and the references therein. For completeness, we provide a proof that leverages the results from [15], in terms of extending Lemma 2.

Lemma 3

([16]) Suppose Assumption 1 holds. Then a multidimensional CLT holds for ${\mathbf{U}}_{n}$; that is,

[TABLE]

with $\bm{\Lambda}=[{{\Lambda_{{ij}}}}]_{i,\,j=1}^{{N^{2}}}$ being an $N^{2}\times N^{2}$ covariance matrix given by

[TABLE]

where ${{\mathbf{I}}}_{ij}$ denotes the $(i,\,j)$-th entry of the identity matrix, and ${{\mathbf{P}}}^{m}_{ij}$ (resp., ${{\mathbf{P}}}^{m}_{ji}$) is the $(i,\,j)$-th (resp., $(j,\,i)$-th) entry of the matrix ${{\mathbf{P}}}^{m}$ (the $m$-th power of ${{\mathbf{P}}}$), $i,\,j=1,\ldots,N^{2}$.

Proof:

See Appendix C. ∎

IV-B Weak Convergence of Test Statistic

In this section, to derive weak convergence results for the test statistic ${D}({{{\bm{\nu}}}\|{{\bm{\pi}}}})$, we will leverage a method commonly-used by statisticians in terms of combining a Taylor’s series expansion for the test statistic and the CLT result for the empirical measure [11]. Recently, under i.i.d. assumptions, such a weak convergence analysis for certain test statistics has been conducted in [10, 5].

To this end, for ${{\bm{\nu}}}\in\mathcal{P}({\bm{\Theta}})$ we consider

[TABLE]

Let ${\mathbf{U}}\sim\mathcal{N}({{\mathbf{0}},\bm{\Lambda}})$ with $\bm{\Lambda}$ given by (13). Now, we are in a position to derive weak convergence results for our test statistic ${D}({{{\bm{\nu}}}\|{{\bm{\pi}}}})$.

Theorem 1

Suppose Assumption 1 holds. Then we have the following weak convergence results:

[TABLE]

where ${\nabla^{2}}h({{\bm{\pi}}})$ is the Hessian of $h({\bm{\nu}})$ evaluated at ${\bm{\nu}}={\bm{\pi}}$, $\rho_{k},\,k=1,\ldots,N^{2}$, are the eigenvalues of the matrix ${\nabla^{2}}h({\bm{\pi}})\bm{\Lambda}$, and $\chi_{1k}^{2},\,k=1,\ldots,N^{2}$, are $N^{2}$ independent $\chi^{2}$ random variables with one degree of freedom.

Proof:

Let us first compute the gradient of $h({\bm{\nu}})$. Expanding the logarithm and after some algebra which leads to cancellations of gradient terms with respect to $\nu_{ij}$ in $\sum_{t=1}^{N}\nu_{it}$, for all $i,j=1,\ldots,N$, we obtain

[TABLE]

which implies

[TABLE]

Further, from (17), we compute the Hessian ${\nabla^{2}}h({{\bm{\nu}}})$ by

[TABLE]

Evaluating all the terms in (23) at ${\bm{\nu}}={{\bm{\pi}}}$ yields ${\nabla^{2}}h({{\bm{\pi}}})$, which will play a crucial role in approximating ${D}({{\bm{\Gamma}_{n}}\|{{{{\bm{\pi}}}}}})$. It is seen that ${\nabla^{2}}h({{\bm{\nu}}})$ is continuous in a neighborhood of ${{\bm{\pi}}}$, and we can utilize the second-order Taylor’s series expansion of $h({{\bm{\nu}}})$ centered at ${{\bm{\pi}}}$ to express ${D}({{\bm{\Gamma}_{n}}\|{{{{\bm{\pi}}}}}})=h(\bm{\Gamma}_{n})-h({{\bm{\pi}}})$. Specifically, by (18) and (5) we have

[TABLE]

where ${{\tilde{\bm{\Gamma}}}_{n}}={\xi_{n}}{\bm{\Gamma}_{n}}+\left({1-{\xi_{n}}}\right){\bm{\pi}}$ is determined with some ${\xi_{n}}\in[0,1]$. From the ergodicity of the chain ${\mathbf{Z}}$ it follows ${\bm{\Gamma}_{n}}\xrightarrow[n\to\infty]{{\textit{w.p.1}}}{{\bm{\pi}}}$, leading to ${\tilde{\bm{\Gamma}}_{n}}\xrightarrow[n\to\infty]{{\textit{w.p.1}}}{{\bm{\pi}}}$. By the continuity of ${\nabla^{2}}h({{\bm{\nu}}})$ we obtain

[TABLE]

Applying Slutsky’s theorem [17], by (12), (24), and (25) we attain

[TABLE]

Finally, by means of a linear transformation [18] on the quadratic form ${{\mathbf{U}}^{\prime}}{\nabla^{2}}h({{\bm{\pi}}}){\mathbf{U}}$, we derive the following alternative asymptotic result:

[TABLE]

where $\rho_{k},\,k=1,\ldots,N^{2}$, are the eigenvalues of the matrix ${\nabla^{2}}h({\bm{\pi}})\bm{\Lambda}$, and $\chi_{1k}^{2},\,k=1,\ldots,N^{2}$, are $N^{2}$ independent $\chi^{2}$ random variables with one degree of freedom. ∎

IV-C Threshold Approximation

We use an empirical Cumulative Distribution Function (CDF) to approximate the actual CDF of $D(\bm{\Gamma}_{n}\|{\bm{\pi}})$. In particular, it is seen from (15) that ${D}({{\bm{\Gamma}_{n}}\|{{\bm{\pi}}}})\simeq(1/(2n)){{\mathbf{U}}^{\prime}}{\nabla^{2}}h({{\bm{\pi}}}){\mathbf{U}}$ for large $n$. Thus, to derive an empirical CDF of ${D}({{\bm{\Gamma}_{n}}\|{{\bm{\pi}}}})$, we can generate a set of Gaussian sample vectors independently according to $\mathcal{N}({{\mathbf{0}},\bm{\Lambda}})$ and then plug each such sample vector into the right-hand side of (15) (i.e., replace ${\mathbf{U}}$), thus, obtaining a set of sample scalars, as a reliable proxy for samples of ${D}({{\bm{\Gamma}_{n}}\|{{\bm{\pi}}}})$.

Once we obtain an empirical CDF of ${D}({{\bm{\Gamma}_{n}}\|{{\bm{\pi}}}})$, say, denoted ${F_{\text{em}}}(\cdot;n)$, then, by (7), we can estimate $\eta_{n,\beta}$ as

[TABLE]

where $F_{\text{em}}^{-1}(\cdot;n)$ is the inverse of $F_{\text{em}}(\cdot;n)$. Note that, the $\eta_{n,\beta}^{\text{wc}}$ derived by (26) depends on the entries of the PL ${\bm{\pi}}$. In practice, if ${{\bm{\pi}}}$ is not directly available, we can replace it by the empirical PL evaluated over a long past sample path. For such cases, we summarize the procedures of estimating the threshold based on our weak convergence analysis as Alg. 1, where $\hat{\bm{\pi}}$ is a good estimate for ${{\bm{\pi}}}$. We note that the length $n_{0}$ of the past sample path should be sufficiently large (e.g., $n_{0}\geq 500N^{2}$) so as to guarantee the validity of taking ${\bm{\pi}}$ to be $\hat{\bm{\pi}}$. In addition, the small positive number $\varepsilon$ (e.g., $\varepsilon\leq 10^{-6}$) introduced in Step 1 is to avoid division by zero, thus ensuring the numerical stability of the algorithm. If, on the other hand, the actual PL ${\bm{\pi}}$ is known, then we can still apply Alg. 1 by replacing the $\hat{\bm{\pi}}$ therein with ${\bm{\pi}}$.

Similar to (26), we can derive another weak convergence based threshold estimator $\bar{\eta}_{n,\beta}^{\text{wc}}$ from (16). However, an easy way of calculating $\bar{\eta}_{n,\beta}^{\text{wc}}$ (also summarized in Alg. 1) still cannot avoid simulations; it is hard to conclude any advantage of $\bar{\eta}_{n,\beta}^{\text{wc}}$ over ${\eta}_{n,\beta}^{\text{wc}}$. As a matter of fact, calculating the eigenvalues of ${\nabla^{2}}h({\bm{\pi}})\bm{\Lambda}$ makes the calculation of $\bar{\eta}_{n,\beta}^{\text{wc}}$ numerically not as stable, compared to the calculation of ${\eta}_{n,\beta}^{\text{wc}}$ via Alg. 1. Other methods for numerically obtaining $\bar{\eta}_{n,\beta}^{\text{wc}}$ can be found, e.g., in [19] and the references therein. Another fact we should point out is that, in [20, p. 30], a slightly different statistic is considered and therefore an even simpler asymptotic distribution can be derived correspondingly. Moreover, some other papers, e.g., [21, 22], also considered similar but different statistics.

We will illustrate by extensive experiments that our weak convergence analysis can empirically produce more accurate estimation of the threshold than Sanov’s theorem for moderate values of $n$; the price we have to pay, however, is a relatively long but still acceptable computation time.

Remark 3

In Alg. 1, due to acceptable numerical errors, the originally estimated $\hat{\bm{\Lambda}}$ (Step 5) could be neither symmetric nor positive semi-definite. Symmetry is imposed by Step 6. Further, to ensure positive semi-definiteness we can diagonalize $\hat{\bm{\Lambda}}$ as

[TABLE]

*where ${\mathbf{O}}$ is an orthogonal matrix and $\text{diag}(\bm{\lambda})$ a diagonal matrix with the elements of $\bm{\lambda}$ in the main diagonal. Due to numerical errors we might encounter cases where some $\lambda_{i}$ are either negative or too small; we can replace them with small positive numbers and recalculate the right-hand side of (27), thus obtaining an updated positive-definite $\hat{\bm{\Lambda}}$. For implementation details, the reader is referred to [23]. ***

IV-D A Robust Hoeffding Test

Many actual systems exhibit time-varying behavior. In this section, we extend our methodology to accommodate such systems and use a set of PLs (instead of a single PL ${\bm{\pi}}$) to model past system activity.

Let the null hypothesis $\mathcal{H}$ be defined as: ${\mathbf{Z}}=\{{{Z_{l}};~{}l=1,2,\ldots}\}$ is drawn according to the set of PLs $\bm{\Pi}=\{{\bm{\pi}}^{\scriptscriptstyle{(1)}},\ldots,{\bm{\pi}}^{\scriptscriptstyle{(L)}}\}\subset\mathcal{P}({\bm{\Theta}})$, i.e., ${\mathbf{Z}}$ is drawn from one of the PLs in $\bm{\Pi}$ but we do not know from which one. Consider a robust version of the Hoeffding test [4, 5, 24] under Markovian assumptions:

[TABLE]

Essentially, the test selects the most likely PL from $\bm{\Pi}$ and uses that to make a decision as in (6). Asymptotic Newman-Pearson optimality of this test is shown in [4].

For $l=1,\ldots,L$, let ${\mathbf{P}}^{{\scriptscriptstyle(l)}}$ denote the transition matrix corresponding to ${\bm{\pi}}^{\scriptscriptstyle{(l)}}$ and, similar to (2), we write

[TABLE]

Assume ${\mathbf{Z}}$ is drawn from PL ${\bm{\pi}}^{\scriptscriptstyle{(l)}}$ which satisfies Assumption 1. Let ${\mathbf{U}}_{n}^{\scriptscriptstyle{(l)}}=\sqrt{n}({{\bm{\Gamma}_{n}}-{{\bm{\pi}}^{\scriptscriptstyle{(l)}}}})$. By Lemma 3, we have

[TABLE]

where \bm{\Lambda}^{\scriptscriptstyle{(l)}}=\big{[}\Lambda^{\scriptscriptstyle{(l)}}_{ij}\big{]}_{i,\,j=1}^{{N^{2}}} is given by

[TABLE]

with ${\mathbf{P}}_{ij}^{\scriptscriptstyle{(l)}m}$ being the $(i,j)$-th entry of the matrix ${\mathbf{P}}^{\scriptscriptstyle{(l)}m}$ (the $m$-th power of ${\mathbf{P}}^{\scriptscriptstyle{(l)}}$). Let ${\mathbf{U}}^{\scriptscriptstyle{(l)}}\sim{\cal{N}}({{\mathbf{0}},{\bm{\Lambda}^{\scriptscriptstyle{(l)}}}})$. Using (15) we obtain

[TABLE]

which leads to an approximation for the infimum term in (28):

[TABLE]

By the right-hand side of (30), we can generate Gaussian samples to compute a reliable proxy for the samples of $\mathop{\inf}_{{\bm{\pi}}\in\bm{\Pi}}D({{\bm{\Gamma}_{n}}\|{\bm{\pi}}})$, thereby, obtaining an empirical CDF, denoted ${F^{\text{rob}}_{\text{em}}}(\cdot;n)$, of $\mathop{\inf}_{{\bm{\pi}}\in\bm{\Pi}}D({{\bm{\Gamma}_{n}}\|{\bm{\pi}}})$. Thus, given a target false positive rate $\beta$, similar to (26), we can estimate the threshold $\eta_{n,\beta}$ as

[TABLE]

where $(F_{\text{em}}^{\text{rob}})^{-1}(\cdot;n)$ denotes the inverse of ${F^{\text{rob}}_{\text{em}}}(\cdot;n)$. Similar to (16), we can also derive a $\chi^{2}$-type asymptotic approximation to the distribution of $\mathop{\inf}\nolimits_{{\bm{\pi}}\in\bm{\Pi}}D({{\bm{\Gamma}_{n}}\|{\bm{\pi}}})$, thus obtaining another WC-based threshold estimator $\bar{\eta}_{n,\beta}^{\text{wc}}$; for economy of space, we omit the details. For the cases where the PLs are not directly available, we summarize the calculation of $\eta_{n,\beta}^{\text{wc}}$ for the robust Hoeffding test as Alg. 2.

V Experimental Results

In this section, we assess the accuracy of our threshold estimator and the performance of the anomaly detection procedure. We start with a numerical evaluation of the threshold’s accuracy and then perform anomaly detection in two application settings using simulated and actual data.

V-A *Numerical Results for Threshold

Approximation*

In this subsection, for simplicity we consider the ordinary (and not the robust) Hoeffding test. We have developed a software package TAM [23] to perform the experiments. We will use ${\bm{\Theta}}=\{1,2,\ldots,N^{2}\}$ to indicate the states and assume the stationary distribution ${{\bm{\pi}}}$ to also be the initial distribution.

In the following numerical examples, we first randomly create a valid (i.e., such that Assumption 1 holds) $N\times N$ transition matrix ${\mathbf{Q}}$, giving rise to an $N^{2}\times N^{2}$ transition matrix ${\mathbf{P}}$, and then generate $T$ test sample paths of the chain ${\mathbf{Z}}$, each with length $n$, denoted ${{\mathbf{Z}}^{\scriptscriptstyle{(t)}}}=\{{Z_{1}^{\scriptscriptstyle{(t)}},\ldots,Z_{n}^{\scriptscriptstyle{(t)}}}\}$, $t=1,\ldots,T$. We use these samples to derive empirical CDF’s. To simulate the case where the PL ${\bm{\pi}}$ is not directly available, we generate one more independent reference sample path ${{\mathbf{Z}}^{\scriptscriptstyle{(0)}}}=\{{Z_{1}^{\scriptscriptstyle{(0)}},\ldots,Z_{n_{0}}^{\scriptscriptstyle{(0)}}}\}$ of length $n_{0}\gg\left|\bm{\Theta}\right|=N^{2}$, thus enabling us to obtain a good estimate of ${\bm{\pi}}$. Note that we do not rely on the test sample paths to estimate the PL ${\bm{\pi}}$. The ground truth ${\bm{\pi}}$ is computed by taking any row of ${\mathbf{P}}^{m_{0}}$ for some sufficiently large $m_{0}$.

Having the ground truth PL ${\bm{\pi}}$ at hand, with the test sample paths ${{\mathbf{Z}}^{\scriptscriptstyle{(t)}}}=\{{Z_{1}^{\scriptscriptstyle{(t)}},\ldots,Z_{n}^{\scriptscriptstyle{(t)}}}\}$, $t=1,\ldots,T$, we can compute $T$ samples of the scalar random variable ${D}({{\bm{\Gamma}_{n}}\|{{{{\bm{\pi}}}}}})$, by (4). Using these samples, we obtain an empirical CDF of ${D}({{\bm{\Gamma}_{n}}\|{{{{\bm{\pi}}}}}})$, denoted ${F}(\cdot;n)$, which can be treated as a dependable proxy of the actual one. The threshold given by (26) with $F_{\text{em}}(\cdot;n)$ replaced by $F(\cdot;n)$ is then taken as a reliable proxy of $\eta_{n,\beta}$. We denote this proxy by $\eta^{*}_{n,\beta}$. To emphasize the dependence on $N$, we write $\eta_{n,\beta}$ (resp., $\eta^{*}_{n,\beta}$) as $\eta_{n,\beta}(N)$ (resp., $\eta^{*}_{n,\beta}(N)$). Next, using the reference sample path ${{\mathbf{Z}}^{\scriptscriptstyle{(0)}}}$ and applying Alg. 1, we obtain $\eta_{n,\beta}^{\text{wc}}(N)$ and $\bar{\eta}_{n,\beta}^{\text{wc}}(N)$.

Let the target false positive rate be $\beta=0.001$. Consider four different scenarios where $N$ is $2$, $4$, $6$, and $8$, respectively. Set $\varepsilon=10^{-10}$, $T=1000$, $m_{0}=1000$, and $n_{0}=1000N^{2}$. Here we note that, in all our experiments, an estimate $\hat{\bm{\pi}}$ for ${\bm{\pi}}$ with $\|\hat{\bm{\pi}}-{\bm{\pi}}\|\leq 10^{-6}$ can be obtained by executing Alg. 1 with parameters $n_{0}\geq 500N^{2}$ and $\varepsilon\leq 10^{-8}$. In Figs. 1(a) through 1(d), the red line plots $\eta^{*}_{n,\beta}(N)$, the blue line $\eta_{n,\beta}^{\text{wc}}(N)$, the magenta line $\bar{\eta}_{n,\beta}^{\text{wc}}(N)$, and the green line $\eta_{n,\beta}^{\text{sv}}(N)$ (cf. (8)), all as a function of the sample size $n$. Setting sample sizes $n$ reasonably small ($n$ should at least be comparable to $N^{2}$), it can be seen that $\eta_{n,\beta}^{\text{wc}}(N)$ and $\bar{\eta}_{n,\beta}^{\text{wc}}(N)$ are more accurate than $\eta_{n,\beta}^{\text{sv}}$, except for the case $N=2$ where all estimators perform approximately equally well. In particular, as $N$ increases, the estimation errors of $\eta_{n,\beta}^{\text{wc}}(N)$ and $\bar{\eta}_{n,\beta}^{\text{wc}}(N)$ are consistently close to zero, while the approximation error of $\eta_{n,\beta}^{\text{sv}}$ increases significantly. Moreover, for the scenarios $N=6,\,8$, $\eta_{n,\beta}^{\text{wc}}(N)$ and $\bar{\eta}_{n,\beta}^{\text{wc}}(N)$ are very close.

Remark 4

In Figs. 1(a)-1(d), the red line representing the “actual” value $\eta^{*}_{n,\beta}$ is not smooth; this is because each time when varying the sample size $n$, we regenerate all the sample paths ${{{\mathbf{Z}}}^{\scriptscriptstyle{(t)}}}=\{{Z_{1}^{\scriptscriptstyle{(t)}},\ldots,Z_{n}^{\scriptscriptstyle{(t)}}}\},~{}t=1,\ldots,T$ from scratch. On the other hand, the blue (resp., magenta) line corresponding to $\eta_{n,\beta}^{\text{wc}}$ (resp., $\bar{\eta}_{n,\beta}^{\text{wc}}$) is smooth because we only need to generate the $T$ Gaussian (resp., $\chi^{2}$-type) sample vectors once. In our experiments, most of the running time is spent generating the sample paths ${{{\mathbf{Z}}}^{\scriptscriptstyle{(t)}}}$ and calculating $\eta^{*}_{n,\beta}$ therefrom. In practice, we will neither generate such samples nor calculate $\eta^{*}_{n,\beta}$, and only need to focus on obtaining $\eta_{n,\beta}^{\text{wc}}$ or $\bar{\eta}_{n,\beta}^{\text{wc}}$, which is computationally not expensive.**

Remark 5

Theoretically speaking, we could use the “actual” threshold $\eta^{*}_{n,\beta}$ as obtained above, but it is of little practical value; the reason is that in statistical anomaly detection applications, we are typically faced with a long series of observations and want to use a so-called windowing technique (see Sec. V-C), which divides the observations into a sequence of detection windows with the same time length. The sample sizes $n$ in different windows may not necessarily be equal, leading to different threshold settings when sliding the windows. If we use the simulated “actual” threshold, then, when varying the detection windows, we will need to regenerate the corresponding samples (for threshold estimation purposes) from scratch, which is computationally too expensive, especially when there are many detection windows. In contrast, to compute our estimator $\eta_{n,\beta}^{\text{wc}}$ (resp., $\bar{\eta}_{n,\beta}^{\text{wc}}$), we only need to generate one set of Gaussian (resp., $\chi^{2}$-type) sample vectors (cf. Remark 4), which can be shared by all the detection windows, thus, saving a lot of computation time. To see this more clearly, let us denote by $\tau_{1}$ the average running time for generating a set of samples with $T$ ($T=1000$ is empirically a good choice) Gaussian (resp., $\chi^{2}$-type) vectors according to (15) (resp., (16)), and $\tau_{2}$ the average running time for calculating a threshold via (26) given the corresponding sample vectors required to derive the empirical CDF. Clearly, we have $\tau_{1}\gg\tau_{2}>0$. Assume we have $W$ detection windows. Then, if we directly simulate the statistic so as to estimate the threshold for each and every detection window, the total running time would be $c_{1}W\tau_{1}+c_{2}W\tau_{2}=(c_{1}\tau_{1}+c_{2}\tau_{2})W$, where $c_{1},c_{2}>0$ are two scaling constants satisfying $c_{1}\tau_{1}\gg c_{2}\tau_{2}$. On the other hand, by simulating Gaussian (resp., $\chi^{2}$-type) samples, the total running time required to estimate all the thresholds for the $W$ detection windows would be $c_{3}\tau_{1}+c_{4}\tau_{2}W$, where $c_{3},c_{4}>0$ are two scaling constants satisfying $c_{4}\approx c_{2}$, leading to $0<c_{4}\tau_{2}\ll c_{1}\tau_{1}+c_{2}\tau_{2}$. Thus, for large $W$ we have $c_{3}\tau_{1}+c_{4}\tau_{2}W\ll(c_{1}\tau_{1}+c_{2}\tau_{2})W$.**

To further investigate the performance of different classes of threshold estimators, we now take the randomness of the transition matrix ${\mathbf{P}}$ into account and define a simulation-based metric $d\left({\hat{\eta},\eta^{*};n,\beta,N,K}\right)$ to quantify the average squared empirical estimation error, specified as follows:

[TABLE]

Recall that $N$ is a parameter representing the number of states in the original chain ${\mathbf{Y}}$. We denote by $\hat{\eta}$ the threshold estimator class (could be $\eta^{\text{sv}}$, $\eta^{\text{wc}}$, or $\bar{\eta}^{\text{wc}}$), and by $\eta^{*}$ a proxy of the actual threshold class (derived by directly simulating the samples of the test statistic). Denote by $K$ the number of independent repetitions of the calculation for $({\hat{\eta}^{\scriptscriptstyle{(k)}}_{n,\beta}({N})-\eta^{\scriptscriptstyle{*(k)}}_{n,\beta}({N})})^{2}$, where $\hat{\eta}^{\scriptscriptstyle{(k)}}_{n,\beta}({N})$ (resp., $\eta^{\scriptscriptstyle{*(k)}}_{n,\beta}({N})$) denotes the class $\hat{\eta}$ (resp., $\eta^{*}$) instantiated under parameters $n$, $\beta$, $N$, and $k\in\{1,\ldots,K\}$.

Setting $\beta=0.001$, $K=200$, $N\in\{2,\,4,\,6,\,8\}$, and $n\in\left\{{\bar{n}=2{N^{2}}+i\times\left\lfloor{0.2{N^{2}}+1}\right\rfloor:\bar{n}<6{N^{2}}+5,i\in\mathbb{N_{+}}}\right\}$, we evaluate $d\left({\hat{\eta},\eta^{*};n,\beta,N,K}\right)$. The results are shown in Fig. 2. Several observations can be made from Figs. 2(a)-2(d): (i) Except for the case $N=2$, both $\eta^{\text{wc}}$ and $\bar{\eta}^{\text{wc}}$ outperform $\eta^{\text{sv}}$, that is, $d\left({\eta^{\text{wc}},\eta^{*};n,\beta,N,K}\right)<d\left({\eta^{\text{sv}},\eta^{*};n,\beta,N,K}\right)$ and $d\left({\bar{\eta}^{\text{wc}},\eta^{*};n,\beta,N,K}\right)<d\left({\eta^{\text{sv}},\eta^{*};n,\beta,N,K}\right)$. (ii) For the cases $N=6,\,8$, $\eta^{\text{wc}}$ and $\bar{\eta}^{\text{wc}}$ perform almost equally well, with both $d\left({\eta^{\text{wc}},\eta^{*};n,\beta,N,K}\right)$ and $d\left({\bar{\eta}^{\text{wc}},\eta^{*};n,\beta,N,K}\right)$ being very close to zero and, for the cases $N=2,\,4$, $\eta^{\text{wc}}$ outperforms $\bar{\eta}^{\text{wc}}$, i.e., $d\left({\eta^{\text{wc}},\eta^{*};n,\beta,N,K}\right)<d\left({\bar{\eta}^{\text{wc}},\eta^{*};n,\beta,N,K}\right)$. (iii) Only for the case $N=2$, $\eta^{\text{sv}}$ performs the best among the three estimators and, $\eta^{\text{wc}}$ performs approximately equally well with $\eta^{\text{sv}}$ in this case. More extensive comparison results can be derived using TAM [23]. We may empirically conclude that, $\eta^{\text{wc}}$ performs consistently the best among the three for almost all scenarios that we have considered and, on the other hand, $\eta^{\text{sv}}$ performs unsatisfactorily when $N>2$, and $\bar{\eta}^{\text{wc}}$ is numerically not as stable as $\eta^{\text{wc}}$, especially for the cases where $N\leq 4$.

V-B *ROC Analysis for the Hoeffding Test with

Different Threshold Estimators*

In this subsection, for simplicity and economy of space, we again only consider the ordinary (and not the robust) Hoeffding test. We note here that similar results can be derived for the robust Hoeffding test. The numerical experiments are conducted using the software package ROM [25].

Let ${\bm{\Theta}}=\{1,2,\ldots,N^{2}\}$ containing $N^{2}$ states. For a given sample size $n$ and a given target false positive rate (FPR) $\beta$, the three thresholds $\eta_{n,\beta}^{\text{wc}}$, $\bar{\eta}_{n,\beta}^{\text{wc}}$, and $\eta_{n,\beta}^{\text{sv}}$, respectively, give rise to three different discrete tests (denote them by “HTWC-1,” “HTWC-2,” and “HTSV,” respectively). To compare their performances, we will conduct the Receiver Operating Characteristic (ROC) [14] analysis (detection rate vs. false alarm rate) using simulated data.

Similar to what we have done in Sec. V-A, we first randomly create a valid $N\times N$ transition matrix ${\mathbf{Q}}$, hence an $N^{2}\times N^{2}$ transition matrix ${\mathbf{P}}$, and then generate $T$ sample paths of the chain ${\mathbf{Z}}$, each with length $n$, denoted ${{\mathbf{Z}}^{\scriptscriptstyle{(t)}}}=\{{Z_{1}^{\scriptscriptstyle{(t)}},\ldots,Z_{n}^{\scriptscriptstyle{(t)}}}\}$, $t=1,\ldots,T$. From ${\mathbf{P}}$ we derive the PL ${\bm{\pi}}$. Next, to simulate anomalies, we create another valid $N\times N$ transition matrix $\bar{{\mathbf{Q}}}$, hence an $N^{2}\times N^{2}$ transition matrix $\bar{{\mathbf{P}}}$, and generate $T$ sample paths of the corresponding chain $\bar{{\mathbf{Z}}}$, each with length $n$, denoted ${{\bar{{\mathbf{Z}}}}^{\scriptscriptstyle{(t)}}}=\{{{\bar{Z}}_{1}^{\scriptscriptstyle{(t)}},\ldots,{\bar{Z}}_{n}^{\scriptscriptstyle{(t)}}}\}$, $t=1,\ldots,T$. Label each sample path of $\bar{{\mathbf{Z}}}$ (resp., ${{\mathbf{Z}}}$) with length $n$ as “positive” (resp., “negative”). Then, $\{{{\mathbf{Z}}^{\scriptscriptstyle{(t)}}}:\,t\in\{1,\ldots,T\}\}\cup\{{\bar{{\mathbf{Z}}}}^{\scriptscriptstyle{(t)}}:\,t\in\{1,\ldots,T\}\}$ will be our test set, which contains $T$ negative (${\mathbf{Z}}^{\scriptscriptstyle{(t)}}$) and $T$ positive (${\bar{{\mathbf{Z}}}}^{\scriptscriptstyle{(t)}}$) sample paths.

Now, by executing Alg. 1 without estimating ${\bm{\pi}}$ (since the ground truth is available), we obtain $\eta_{n,\beta}^{\text{wc}}$ and $\bar{\eta}_{n,\beta}^{\text{wc}}$. Also, by (8) we obtain $\eta_{n,\beta}^{\text{sv}}$. For each sample path in the test set, we compute ${D}({{\bm{\Gamma}_{n}}\|{{{{\bm{\pi}}}}}})$ by (4). Next, using $\eta_{n,\beta}^{\text{wc}}$ (resp., $\bar{\eta}_{n,\beta}^{\text{wc}}$, $\eta_{n,\beta}^{\text{sv}}$), we can apply HTWC-1 (resp., HTWC-2, HTSV) to detect each sample path as positive or negative. Then, we integrate these reports with the ground truth labels so as to calculate the true positive rate (TPR) and FPR, thereby, obtaining a point of the ROC space.

In our experiments, we take $T=1000$. Fig. 3(a) (resp., 3(b)) shows the ROC graphs of HTWC-1, HTWC-2, and HTSV for a scenario corresponding to $N=4,\,n=50$ (resp., $N=6,\,n=100$); different points on the graph are obtained by $\beta$ taking values from a predesignated finite set $\{0.001\}\cup\{0.01,0.02,\ldots,0.19\}$. It is seen from Fig. 3(a) (or Fig. 3(b)) that all TPR values are very close to 1, which is good, but for most cases (each case corresponds to a specific “small” target FPR $\beta$) HTWC-1 and HTWC-2 have much closer FPR values to the target FPR value than HTSV, meaning HTWC-1 and HTWC-2 are able to control for false alarms better than HTSV. To see this more clearly, we show a few specific values of the (TPR, FPR) pair in Tables I and II. It is worth noting that in the $N=6$ scenario, HTSV is almost a random guess for all the target FPR cases that are considered. More extensive experiments show that, as $N$ increases, the performance of HTSV gets worse and worse; in particular, when $N\geq 6$, HTSV is very likely merely a random guess yielding an ROC point close to $(1,1)$. During our experiments, another observation is that, for each fixed $N$ and $\beta$, when $n$ increases, all HTWC-1, HTWC-2, and HTSV perform better and better; this is because with larger sample sizes, all the three estimators $\eta_{n,\beta}^{\text{wc}}$, $\bar{\eta}_{n,\beta}^{\text{wc}}$, and $\eta_{n,\beta}^{\text{sv}}$ approximate the actual $\eta_{n,\beta}$ better. We therefore conclude that HTWC-1 (or HTWC-2) typically outperforms HTSV in the sense that the former has a better capability of controlling the false alarm rate (i.e., FPR) while maintaining a satisfactory detection rate (i.e., TPR).

Remark 6

A natural concern about the ROC analysis above might be the setting of the target FPR ($\beta$) values; one may ask: How about always setting $\beta$ to a “very small” value, say, $10^{-10}$, $10^{-100}$, or even $10^{-1000}$? We have actually already discussed this partly in Sec. III. Setting a too small $\beta$ would typically lead to an unsatisfactory detection rate (TPR). In addition, note that $\eta_{n,\beta}^{\text{wc}}$ (or $\bar{\eta}_{n,\beta}^{\text{wc}}$) is numerically obtained from an empirical CDF, say, $G(x)$, of some scalar random variable; we have $G(x)$ nondecreasing, and $\mathop{\lim}\nolimits_{x\to+\infty}G\left(x\right)=1$, implying that finding an “accurate” $x$ such that $G(x)=1-\beta$ would be hard for a too small $\beta\in(0,1)$. An empirically “good” choice of $\beta$ is 0.001 (see Tabs. I and II), which is what we use in our applications. Because HTWC-1 and HTWC-2 perform almost equally well in our experiments, but HTWC-1 is more stable and less computationally demanding, we will only apply HTWC-1 in the following.**

V-C Simulation Results for Network Anomaly Detection

In this subsection we test our approach in a communication network traffic anomaly detection application. We will use the term traffic and flow interchangeably. We perform the simulations using the software package SADIT [8], which, based on the fs-simulator [26], is capable of efficiently generating flow-level network traffic datasets with annotated anomalies.

As shown in Fig. 4, the simulated network consists of an internal network involving eight normal users (CT1-CT8), a server (SRV) that stores sensitive information, and three Internet nodes (INT1-INT3) that connect to the internal network via a gateway (GATEWAY).

As in [4, Sec. III.A], to characterize the statistical properties of the flow data, we use as features the flow duration and size (bits). We also cluster the source/destination IP addresses and use as features for each flow the assigned cluster ID and the distance of the flow’s IP from the cluster center. For each feature, we quantize its values into discrete symbols so as to obtain a finite alphabet $\bm{\Xi}$, hence ${\bm{\Theta}}$, for our model. Based on the time stamps (the start times) of the flows, we divide the flow data into a series of detection windows, each of which contains a set of flow observations (see [4] for details).

To implement our anomaly detection approach, we first estimate a PL ${\bm{\pi}}$ (resp., a PL set $\bm{\Pi}$) from the stationary (resp., time-varying) normal traffic. Note that, for either case, the reference data should be anomaly-free ideally. However, in our experiments, for the stationary case we use as reference traffic the entire flow sequence with anomalies injected at some time interval; this makes sense because the size of a typical detection window is much smaller than that of the whole flow sequence and the fraction of anomalies is indeed very small, leading to an estimation for the PL with acceptable accuracy. On the other hand, for the time-varying case we generate the reference traffic without anomalies and the test traffic with anomalies separately, sharing all the parameter settings in the statistical model used in SADIT except the ones for introducing anomalies. Note that, estimating a PL for the stationary traffic is relatively easy, while, for the time-varying traffic, we need to make an effort to estimate several different PLs corresponding to certain periods of the day. We apply the two-step procedure proposed in [4]; that is, we first generate a relatively large PL set and then refine the candidate PLs therein by solving a weighted set cover problem. Note also that, if we already know the periodic system activity pattern, then we can directly estimate the PL set period by period; see another anomaly detection application in Sec. V-D for example.

Now, having the reference PL (resp., PL set) at hand, we persistently monitor the test traffic and report an anomaly instantly as long as the relative entropy $D({{{\bm{\Gamma}_{n}}}\|{{\bm{\pi}}}})$ (resp., $\mathop{\inf}_{{\bm{\pi}}\in\bm{\Pi}}D({{\bm{\Gamma}_{n}}\|{\bm{\pi}}})$) exceeds the threshold $\eta_{n,\beta}^{\text{wc}}$ for the current detection window, where $n$ is the number of flow samples within the window. It is worth pointing out that, for the current application, we will not seek to identify which flows belonging to an abnormal detection window contribute mostly to causing the anomaly, but, in some other applications, e.g., the one in Sec. V-D, we will do so.

In the following, we consider two scenarios – one for stationary traffic and the other for time-varying traffic.

V-C1 Stationary Network Traffic – Scenario

V-C-1

We mimic anomalies caused by a large file download [3, Sec. IV.A.2]. The simulation time is $7000$ s. A user increases its mean flow size to $10$ times the usual value between $1000$ s and $1500$ s. The interval between the starting points of two consecutive time windows is taken as $w_{\text{d}}=50$ s, the window-size is set to $w_{\text{s}}=200$ s, and the target false positive rate is set to $\beta=0.001$. The number of user clusters is $k=2$ and the quantization level for flow duration, flow size, and distance to cluster center is set to $n_{1}=1$, $n_{2}=2$, and $n_{3}=2$, respectively. Thus, the original chain has $N=2\times 1\times 2\times 2=8$ states, and we have $N^{2}=64$ states in the transformed chain.

The detection results are shown in Figs. 5(a) and 5(b), both of which depict the relative entropy (divergence) metric defined in (4). The green dashed line in Fig. 5(a) is the threshold estimated using Sanov’s theorem (i.e., $\eta_{n,\beta}^{\text{sv}}$ given by (8), where $n$ is the sample size in each specific detection window). The green dashed line in Fig. 5(b) is the threshold given by our estimator (i.e., $\eta_{n,\beta}^{\text{wc}}$ computed by Alg. 1). The interval during which the divergence curve is above the threshold line (the red segment) corresponds to the time instances reported as abnormal. Fig. 5(a) shows that, if $\eta_{n,\beta}^{\text{sv}}$ is used as the threshold, then the Hoeffding test reports too many false alarms, and, Fig. 5(b) shows that, if, instead, we use $\eta_{n,\beta}^{\text{wc}}$ as the threshold, then the Hoeffding test does not report any false alarm while successfully identifying the true anomalies between $1000$ s and $1500$ s.

V-C2 Time-Varying Network Traffic – Scenario

V-C-2

Consider the case where the network in Fig. 4 is simulated with a day-night traffic pattern in which the flow size follows a log-normal distribution. We use precisely the same scenario as that in [4, Sec. IV.B.2]. The ground truth anomaly (consider an anomaly where node CT2 increases its mean flow size by 30%) is injected beginning at $59$ h and lasting for $80$ minutes.

Using the two-step procedure proposed in [4, Sec. III.C], we first obtain $32$ rough PL candidates. Then, using the PL refinement algorithm given in [4, Sec. III.D] equipped with the cross-entropy threshold parameter $\lambda=0.028$, which is determined by applying Alg. 2, we finally obtain $6$ PLs, being active during morning, afternoon, evening, mid-night, dawn, and the transition time around sunrise, respectively. Note that, since we have obtained the PL set in a different way, in the following, when applying Alg. 2 for each detection window, we can skip the first two steps (lines 2 and 3). In the subsequent detection procedure, the chief difference between our method and the one used in [4] is that we no longer set the threshold universally as a constant; instead, we calculate the threshold $\eta_{n,\beta}^{\text{wc}}$ for each detection window using Alg. 2. Set $k=1$, $n_{1}=1$, $n_{2}=4$, and $n_{3}=1$. Thus, the original chain has $N=1\times 1\times 4\times 1=4$ states, and we have $N^{2}=16$ states in the transformed chain for this case. Take $w_{\text{d}}=1000$ s, $w_{\text{s}}=1000$ s, and $\beta=0.001$. We see from Fig. 6 that the anomaly is successfully detected, without any false alarms.

V-D Anomaly Detection for Waze Jams

V-D1 Dataset Description

The Waze datasets under investigation are kindly provided to us by the Department of Innovation and Technology (DoIT) in the City of Boston. The datasets include three parts: the jam data $\mathcal{J}_{1}$ (traffic slowdown information generated by Waze based on users’ location and speed; note that each jam consists of a set of points), the corresponding point data $\mathcal{J}_{2}$ (latitudes and longitudes of the points within jams), and the alert data $\mathcal{J}_{3}$ (traffic incidents reported by active users; we will call such a user a “Wazer”). For each part, we only list the features that we have used in our algorithms. In particular, each entry (jam) in $\mathcal{J}_{1}$ has the following fields: uuid (unique jam ID), start time, end time, speed (current average speed on jammed segments in meters per second), delay (delay caused by the jam compared to free flow speed, in seconds), and length (jam length in meters). The information for each entry in $\mathcal{J}_{2}$ includes a jam uuid and the locations (latitudes and longitudes) of the points within the jam. The fields of each entry in $\mathcal{J}_{3}$ include: uuid (unique system ID; this is different from the jam ID in $\mathcal{J}_{1}$), location (latitude and longitude per report), type (event type; e.g., accident, weather hazard, road closed, etc.), start time, and end time. It is seen that, by combing $\mathcal{J}_{1}$ and $\mathcal{J}_{2}$, we can denote each jam in $\mathcal{J}_{1}$ as

[TABLE]

where $i$ is the index, uuid is the unique jam ID, “loc” (resp., “startTime”) is the abbreviation for location (resp., start time). Because we are only interested in detecting the abnormal jams in real-time, we will not use the jam end times.

V-D2 Anomaly Description

Typically we can observe lots of jams in certain areas during rush hour, e.g., the AM/PM peaks, and most of them are “normal” except those with extremely atypical features (delay, length, etc.). On the other hand, if a jam was observed outside of rush hours or typical areas, then it would likely be “abnormal.”

V-D3 Description of the Experiments

Treating Waze jams as a counterpart of the network flows in Sec. V-C, we implement the robust Hoeffding test on the quantized jam data in the following experiments.

Consider an area around the Boston University (BU) bridge, whose location is specified by latitude and longitude $(42.351848,-71.110730)$ (see the green marker in Fig. 7). Extract the jam data no farther than 3 kilometers from BU (within the circle in Fig. 7). Note that it is possible for Waze to report several jams at the same time. To assign each jam a unique time stamp, we slightly perturb the start time of the jams that share the same time stamp in the raw data. Such slight adjustments would not alter the original data significantly.

Reference (resp., test) data are taken as jams reported on March 9, 2016 (resp., March 16, 2016). Both dates are Wednesdays, representing typical workdays. There are 3218 jams in the reference data, and 3882 jams in the test data. Note that we have historical data for a relatively long time period (compared to the test data within a detection window); including all the jams reported within the selected reference time period would not hurt the accuracy of the PLs (anomaly-free ideally) to be estimated.

The features that we use for anomaly detection are location, speed, delay, and length. The time stamp of a jam is taken as its start time. To quantize the location, we need to define the distance between two jams. For any valid index $i$, denote the complete location data of jam $i$ by

[TABLE]

where $x$’s and $y$’s denote the latitudes and longitudes, respectively, and $i_{n}$ is the number of points in jam $i$ (typically, $i_{n}$ is greater than $4$). Noting that most of the jams are approximately linear in shape, we simplify (33) by using the 4 vertices of the “smallest” rectangle that covers all the points in the jam and update (33) by

[TABLE]

where ${x_{i,\min}}=\min\{{{x_{i,1}},\ldots,{x_{i,{i_{n}}}}}\}$, ${x_{i,\max}}=\max\{{{x_{i,1}},\ldots,{x_{i,{i_{n}}}}}\}$, ${y_{i,\min}}=\min\{{{y_{i,1}},\ldots,{y_{i,{i_{n}}}}}\}$, and ${y_{i,\max}}=\max\{{{y_{i,1}},\ldots,{y_{i,{i_{n}}}}}\}$. Note that $\text{loc}[i]$ in (34) only contains $4$ points. Denote the point-to-point distance (in meters) yielded by Vincenty’s formula [27] as $d_{V}(\cdot,\cdot)$. Then, for any pair of jams, say, indexed $i$ and $j$, we define the distance between them as

[TABLE]

Using the distance defined above and setting the quantization level for “location” as $3$, we apply the commonly used $K$-means clustering method [28], thus obtaining 3 cluster centers as depicted in Fig. 7 (note that, by (34) each cluster center is represented by $4$ red markers).

In all our experiments, we take the quantization level to be 1 for “speed,” and set the target false alarm rate as $\beta=0.001$. The window size is taken as $w_{\text{s}}=10$ minutes, and the distance between two consecutive windows is $w_{\text{d}}=5$ minutes. To estimate the PLs, we divide a whole day into 4 subintervals: 5:00-10:00 (AM), 10:00-15:00 (MD), 15:00-19:00 (PM), and 19:00-5:00 (NT). So, for each scenario we end up with $4$ PLs, corresponding to the AM peak, the middle day, the PM peak, and the night, respectively. To calculate the threshold $\eta_{n,\beta}^{\text{wc}}$ for each detection window, we use Alg. 2.

V-D4 Detection Results

First, let the quantization level for “delay” be 2 and for “length” be 1. The original sample path has $N=3\times 1\times 2\times 1=6$ states. Thus, we have $N^{2}=36$ states in the transformed chain. We use relatively sparse quantization levels for “delay” and “length” to avoid unnecessary computational overhead in the quantization subroutine for the jam location data. After running our algorithm in the initial step, 910 out of 3882 jams are reported within abnormal detection windows, which correspond to the red segments in Fig. 8(a). We then perform a refinement procedure by selecting jams in these windows with non-typical individual features as follows. For each selected feature, we calculate the sample mean $\mu$ and sample standard deviation $\sigma$ using the reference data. We then label as anomalous any jam with feature value exceeding $\mu+3\sigma$. We first consider the delay feature. Using the $3\sigma$-rule on delay, we obtain an anomaly list $\mathcal{L}_{1}$ containing $4$ jams.

Second, let the quantization level for “delay” be 1 and for “length” be 2. Then, again, the original sample path has $N=3\times 1\times 1\times 2=6$ states, and we have $N^{2}=36$ states in the transformed chain. After rerunning the algorithm in the initial step, 590 out of 3882 jams are reported within abnormal detection windows, which correspond to the red segments in Fig. 8(b), and, after refining by use of the $3\sigma$-rule on the feature “length”, we end up with an anomaly list $\mathcal{L}_{2}$ containing $2$ jams.

Finally, we take $\mathcal{L}=\mathcal{L}_{1}\cup\mathcal{L}_{2}$ as our ultimate anomaly list, which contains $6$ jams in total. By checking the time stamps and the alarm instances, we see that all of these 6 jams would be reported as abnormal by our method within $5$ minutes from their start time; this is satisfactory in a real-time traffic jam anomaly detection application. Note that we can tune $w_{\text{d}}$ and $w_{\text{s}}$ such that the detection becomes even faster while maintaining good accuracy in identifying anomalies. Specifically, smaller $w_{\text{d}}$ leads to faster detection while $w_{\text{s}}$ should be reasonably big (the number of jams in a window should at least be comparable to $N^{2}$). By comparing the locations and time stamps, we map the jams in the final anomaly list to the alert data $\mathcal{J}_{3}$, and find that one of them was reported by Wazers as “road closed,” another as “weather hazard,” and all the others as “jam heavy traffic.” In addition, all of them occurred during non-peak hours. We list the key features of these abnormal jams in Tab. III, where the atypical values of the features “delay” and “length” have been highlighted in bold red. It is worth pointing out that jam 2412 is reported as abnormal based on “delay,” but its length (highlighted in bold magenta) is also above the threshold for refining the detection results based on “length.” Note also that the latitude and longitude in each row of Tab. III represent the closest location of the Wazer who reported the alert for the corresponding jam (extracted from the alert data $\mathcal{J}_{3}$); the shapes of the actual jams have been visualized as colored bold curves in Fig. 7. While in this application we do not have ground truth, it is reassuring that the jams we identify as anomalous have indeed been reported as non-typical by Wazers. Clearly, depending on how such a detection scheme will be used by a City’s transportation department, our approach provides flexibility in setting thresholds to adjust the volume of reported anomalous jams. This volume will largely depend on the resources that City personnel have to further investigate anomalous jams (e.g., using cameras) and intervene.

Remark 7

If we directly apply the $3\sigma$-rule on the whole test data without implementing the Hoeffding test to obtain a potential anomaly list first, then we would very likely end up with too many anomalies, which might include undesirable false alarms. Indeed, when we apply the $3\sigma$-rule on the whole test data for “delay” (resp., “length”), we obtain 38 (resp., 62) “anomalies,” which are much more than those in our final anomaly list (6 only). Thus, including the well-validated Hoeffding test in our method ensures a good control of false alarms.**

VI Conclusions and Future Work

We have established weak convergence results for the relative entropy in the Hoeffding test under Markovian assumptions, which enables us to obtain a tighter estimator (compared to the existing estimator based on Sanov’s theorem) for the threshold needed by the test. We have demonstrated good performance of our estimator by applying the Hoeffding test in extensive numerical experiments for the purpose of statistical anomaly detection. The application scenarios involve not only simulated communication networks, but also real transportation networks. Our work contributes to enhancing cyber security and helping build smarter cities.

As for future work, it is of interest to establish theoretical comparison results concerning the tightness of the threshold estimators. The challenge in this direction arises from associating the finite sample-size setting with the asymptotic properties of the Central Limit Theorem and the large deviations results (Sanov’s theorem). It is also of interest to conduct rigorous analysis relating the computation time of the proposed estimation approach to its accuracy. Also, it is possible to consider additional applications.

Appendix A Proof of Lemma 1

Expanding the first $N$ entries of ${{\bm{\pi}}}{\mathbf{P}}={{\bm{\pi}}}$, we obtain $q_{1i}\sum_{t=1}^{N}{{\pi_{t1}}}={\pi_{1i}}$, $i=1,\ldots,N$. Summing up both sides of these equations, it follows

[TABLE]

Noticing $\sum_{i=1}^{N}q_{1i}=1$, (A.1) implies $\sum_{t=1}^{N}{{\pi_{t1}}}=\sum_{t=1}^{N}{{\pi_{1t}}}$, which, together with $q_{11}\sum_{t=1}^{N}{{\pi_{t1}}}={\pi_{11}}$, yields

[TABLE]

Similarly, we can show (9) holds for all the other $(i,\,j)$’s.

Appendix B Proof of Lemma 2

This can be established by applying [15, Corollary 1]. Noting $f_{k}(\cdot)$ is an indicator function, thus Borel measurable and bounded, and the chain ${\mathbf{Z}}$ is uniformly ergodic, we see that, $\exists B\in({0,\infty})$ s.t. $|{{f_{k}}(Z)}|\leq B,\,\forall Z$, implying that $\mathbbm{E}[{{{|{{f_{k}}(Z)}|}^{3}}}]\leq{B^{3}}<\infty$, and [15, (3)] holds with $M(\cdot)$ bounded, leading to $\mathbbm{E}[M]<\infty$, and $\gamma(n)=t^{n}$ for some $t\in(0,1)$, indicating that $\sum\nolimits_{n}{{{({\gamma(n)})}^{{1\mathord{/{\vphantom{13}}.\kern-1.2pt}3}}}}=\sum\nolimits_{n}{{t^{{n\mathord{/{\vphantom{n3}}.\kern-1.2pt}3}}}}=\sum\nolimits_{n}{{{({{t^{{1\mathord{/{\vphantom{13}}.\kern-1.2pt}3}}}})}^{n}}}<\infty$. Thus, all the conditions needed by [15, Corollary 1] are satisfied.

Appendix C Proof of Lemma 3

We can directly extend Lemma 2 to the multidimensional case (see [29, Chap. 8]). In particular, under Assumption 1, (12) holds with $\bm{\Lambda}$ given by

[TABLE]

where ${\bm{\Lambda}^{\scriptscriptstyle{(0)}}}$ and ${\bm{\Lambda}^{\scriptscriptstyle{(m)}}}$ are specified, respectively, by

[TABLE]

Let the subscript $ij$ denote the $(i,j)$ elements of the matrices $\bm{\Lambda},\bm{\Lambda}^{\scriptscriptstyle{(0)}},\bm{\Lambda}^{\scriptscriptstyle{(m)}}$. By the Markovian properties, after some direct algebra, for $i,\,j=1,\ldots,{N^{2}}$ we obtain $\Lambda^{\scriptscriptstyle{(0)}}_{ij}={{\tilde{\pi}}_{i}}({{{{\mathbf{I}}}_{ij}}-{{\tilde{\pi}}_{j}}})$ and

[TABLE]

Acknowledgments

We thank Jing Wang for his contributions and help in developing the software package SADIT [8]. We also thank the DoIT of the City of Boston, Chris Osgood, Alex Chen, and Connor McKay for supplying the Waze data. We thank Athanasios Tsiligkaridis for his help in deriving Tab. III in Sec. V-D. We finally thank the anonymous reviewers for useful comments on preliminary versions of this paper.