An Improved Composite Hypothesis Test for Markov Models with
Applications in Network Anomaly Detection∗

arXiv:1509.01706v3 [] 16 Sep 2016

Jing Zhang† and Ioannis Ch. Paschalidis‡
Abstract— Recent work has proposed the use of a composite
hypothesis Hoeffding test for statistical anomaly detection.
Setting an appropriate threshold for the test given a desired
false alarm probability involves approximating the false alarm
probability. To that end, a large deviations asymptotic is
typically used which, however, often results in an inaccurate
setting of the threshold, especially for relatively small sample
sizes. This, in turn, results in an anomaly detection test that
does not control well for false alarms. In this paper, we develop
a tighter approximation using the Central Limit Theorem
(CLT) under Markovian assumptions. We apply our result to
a network anomaly detection application and demonstrate its
advantages over earlier work.
Index Terms— Hoeffding test, weak convergence, Markov
chains, network anomaly detection.

I. I NTRODUCTION
During the last decade, the applications of statistical
anomaly detection in communication networks using the
Large Deviations Theory (LDT) [1] have been extensively
explored; see, e.g., [2], [3], [4], among others. Statistical
anomaly detection involves characterizing the “normal behavior” of the system and identifying the time instances
corresponding to abnormal system behavior. Assuming that
the network traffic is stationary in time, [2] applies two
methods to that end. The first method, termed “model-free,”
models traffic as an independent and identically distributed
(i.i.d.) sequence. The second method, termed “model-based,”
models traffic as a finite-state Markov chain. These are then
extended in [4] to the case where the network traffic is
time-varying. Essentially, each of these methods is designed
to tackle a certain Universal Hypothesis Testing Problem
(UHTP).
A UHTP 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].
It is well known that the test proposed by Hoeffding [6]
is optimal in an error exponent sense [1], [2], [4]. When
implementing Hoeffding’s test, one must set a threshold η,
which can be estimated by using Sanov’s theorem [1], a
result within the scope of LDT. Note that such an estimate
(let us denote it by η sv ) is valid only in the asymptotic sense.
* Research partially supported by the NSF under grants CNS-1239021
and IIS-1237022, by the ARO under grants W911NF-11-1-0227 and
W911NF-12-1-0390, by the ONR under grant N00014-10-1-0952, and by
the Cyprus Research Promotion Foundation under Grant New Infrastructure
Project/Strategic/0308/26.
† Division of Systems Engineering, Boston University, jzh@bu.edu.
‡
Department
of
Electrical
and
Computer
Engineering
and Division of Systems Engineering, Boston University, 8
St. Mary’s St., Boston, MA 02215, yannisp@bu.edu,
http://sites.bu.edu/paschalidis.

In practice, however, only a finite number of observations are
available, and it can be shown by experiments (e.g., using
the software package TAHTIID [7]) that η sv is typically not
accurate enough.
To improve the accuracy of the estimation for η, [5] (or,
see [8]) borrows a method typically used by statisticians;
that is, deriving results based on Weak Convergence (WC)
of the test statistic to approximate the error probabilities
of Hoeffding’s test. Under the i.i.d. assumption, by using
the technique combining WC results together with Taylor’s
series expansion, [5] (or, see [8], [9]) gives an alternative
approximation for η (let us denote it by η wc ). Using [7],
one can verify that, in the finite sample-size setting, η wc is
typically much more accurate than η sv .
It is worth pointing out that some researchers have also
tried to obtain a more accurate approximation for η by
refining Sanov’s theorem [10]. However, as noted in [8],
such refinements of large deviation results are typically hard
to calculate.
In this work we establish η wc under a Markovian assumption, thus, extending the work of [5] which considered
the i.i.d. case. We apply the proposed procedure to network
anomaly detection by embedding it into the software package
SADIT [11].
The rest of the paper is organized as follows. We first
formulate the problem in Sec. II and derive the theoretical
results in Sec. III. We then present the experimental results
in Sec. IV and finally provide concluding remarks in Sec. V.
Notation: By convention, all vectors are column vectors.
To save space, we write x = (x1 , . . . , xn ) to denote the column
vector x. We use prime to denote the transpose of a matrix
or vector.
II. P ROBLEM F ORMULATION
Let Ξ = {ξi ; i = 1, . . . , N} be a finite alphabet. The observations Y = {Yl ; l = 0, 1, 2, . . .} form a stochastic process
taking values from Ξ.
In the model-based UHTP [1], [2], [4], under the null
hypothesis H , the observations Y = {Yl ; l = 0, 1, 2, . . .} are
drawn according to a Markov chain with state set Ξ. Let its
transition matrix be Q = (qi j )Ni,j=1 .
Let 1{·} denote the indicator function. Define the empirical
PL by
Γn (θi j ) =

1 n
1
,
n ∑l=1 {Zl =θi j }

(1)

where Zl = (Yl−1 , Yl ) , l = 1, . . . , n, θi j = (ξi , ξ j ) ∈ Ξ ×
Ξ, i, j = 1, . . . , N, and denote the new alphabet Θ =



θi j ; i, j = 1, . . . , N = θ̃k ; k = 1, . . . , N 2 . Note that Θ =
Ξ × Ξ and θ̃1 = θ11 , . . ., θ̃N = θ1N , . . ., θ̃(N−1)N+1 = θN1 , . . .,
θ̃N 2 = θNN . Let also the set of PLs on Θ be P(Θ).
It is seen that the transformed observations Z =
{Zl ; l = 1, 2, . . .} form a Markov chain evolving on Θ. Let its
2
transition matrix be P = (pi j )Ni, j=1 with stationary distribution

π = (πi j ; i, j = 1, . . . , N) = π̃k ; k = 1, . . . , N 2 ,

to be π is for notational simplicity; our results apply for any
feasible initial distribution. Note also that, under (C), π must
have full support over Θ, meaning each entry in π is strictly
positive.

where πi j denotes the probability of seeing θi j , and π̃1 = π11 ,
. . ., π̃N = π1N , . . ., π̃(N−1)N+1 = πN1 , . . ., π̃N 2 = πNN . Let
p (θi j |θkl ) denote the transition probability from θkl to θi j .
Then we have

Proof: Expanding the first N entries of π P = π , we
N
πt1 = π1i , i = 1, . . . , N. Summing up both sides
obtain q1i ∑t=1
of these equations, it follows


N
N
N
(8)
∑i=1 q1i · ∑t=1 πt1 = ∑t=1 π1t .

p (θi j |θkl ) = 1{i=l} qi j ,

k, l, i, j = 1, . . . , N,

which enables us to obtain P directly from Q.
Write π = (π11 , . . . , π1N , . . . , πN1 , . . . , πNN ), and Γ n =
(Γn (θ11 ), . . . , Γn (θ1N ), . . . , Γn (θN1 ), . . . , Γn (θNN )).
Define the relative entropy (i.e., divergence)
.
N N
Γn (θi j ) ∑Nj=1 Γn (θi j )
.
D ( Γ n k π ) = ∑ ∑ Γn (θi j ) log
, (2)
i=1 j=1
πi j ∑Nj=1 πi j
and the empirical measure
√
Γn − π ) .
Un = n (Γ

(3)

For any positive integer n, let Hn be the output of the
test that decides to accept or to reject the null hypothesis H
based on the first n observations in the sequence Z. Under
the Markovian assumption, Hoeffding’s test [1] is given by
Hn = H ⇔ D ( Γ n k π ) ≤ ηn ,

(4)

where D ( Γ n k π ), defined in (2), is the test statistic, and ηn
is a threshold. Note that, when applied in anomaly detection,
Hoeffding’s test will report an anomaly if and only if Hn 6=
H , or, equivalently, D ( Γ n k π ) > ηn .
Under hypothesis H , the false alarm probability [5] of
the test sequence {Hn } as a function of n is defined by
βn = PH {D ( Γ n k π ) > ηn } ,

(5)

where PH {A} denotes the probability of event A under hypothesis H . Sanov’s theorem [1] implies the approximation
for ηn
1
(6)
ηnsv ≈ − log (βn ) .
n
Given βn ∈ (0, 1), our central goal in this paper is to improve
the accuracy of the approximation for ηn .
III. T HEORETICAL R ESULTS
We introduce the following condition:
Condition (C). Z = {Zl ; l = 1, 2, . . .} is an aperiodic, irreducible, and positive recurrent Markov chain [12], [13]
evolving on Θ with stationary distribution π , and with the
same π as its initial distribution.
Remark 1: Since Θ is a finite set, we know that Z is uniformly ergodic [12], [14]. Assuming the initial distribution

Lemma III.1 Assume (C) holds. Then
πi j
πi j
= N
= qi j , i, j = 1, . . . , N.
N
π
∑t=1 it
∑t=1 πti

(7)

N
Noticing ∑Ni=1 q1i = 1, we see that (8) implies ∑t=1
πt1 =
N
N
∑t=1 π1t , which, together with q11 ∑t=1 πt1 = π11 , yields
π11
π11
= N
= q11 .
N
∑t=1 π1t
∑t=1 πt1

Similarly, we can show (7) holds for all the other (i, j)’s.
Remark 2: By Remark 1 and Lemma III.1 we see that,
under Condition (C), all entries in Q are strictly positive,
indicating that any two states of the original chain Y are
connected. This is a strict condition; yet, in engineering
practice, if some πi j in (7) are equal to zero, we can always
replace them with a small positive value, and then normalize
the modified vector π , thus ensuring that Condition (C) can
be satisfied. See Sec. IV-A for an illustration of this trick.
A. Weak convergence of Un
Let us first establish CLT results for one-dimensional
empirical measures


√
Un(k) = n Γn θ̃k − π̃k , k = 1, . . . , N 2 .
(9)
Fix k ∈ {1, . . . , N 2 }. Define
(
fk (Z) = 1{Z=θ̃k } =

1, if Z = θ̃k ,

0, if Z ∈ Θ\ θ̃k .

(10)

Lemma III.2 Assume (C) holds. Then a CLT holds

d.
for Un(k) ; that is, Un(k) −−−→ N 0, σk2 with σk2 =
n→∞
Cov ( fk (Z1 ) , fk (Z1 )) + 2 ∑∞
m=1 Cov ( f k (Z1 ) , f k (Z1+m )) <
∞.
Proof: This can be established by applying [12, Corollary 1]. Noting fk (·) is bounded and the chain Z is uniformly
ergodic, we see that all the conditions needed by [12,
Corollary 1] are satisfied.
Now we can directly extend Lemma III.2 to the multidimensional case (for an informal argument, see, e.g., [15];
cf. [16, Sec. III.6] for the definition of a Gaussian random
variable with a general normal distribution). Due to its role in
the following applications, we state the result as a theorem.
Theorem III.3 Assume (C) holds. Then a multidimensional
CLT holds for Un ; that is,
d.

Un −−−→ N (0, Λ )
n→∞

(11)

with Λ being an N 2 × N 2 covariance matrix given by
∞

Λ = Λ0 + ∑m=1 Λm ,

(12)

where Λ 0 and Λ m are specified, respectively,
2
by
Λ 0 = [Cov ( fi (Z1 ) , f j (Z1 ))]Ni, j=1
and
Λm =
2

[Cov ( fi (Z1 ) , f j (Z1+m )) + Cov ( f j (Z1 ) , fi (Z1+m ))]Ni, j=1 ,
m = 1, 2, . . ..

Further, from (17), we compute the Hessian ∇2 h (νν ) by
considering three cases:
2 ν)
= 0.
(i) If k 6= i, then ∂∂νi h(ν
j ∂ νkl
(ii) If k = i and l = j, then
∂ 2 h (νν )
∂ 2 h (νν )
1
1
.
=
−
=
∂ νi j ∂ νkl
νi j ∑Nj=1 νi j
∂ νi2j
(iii) If k = i and l 6= j, then

The rest of this subsection is devoted to the computation of

N 2
Λ ; we will express Λ = Λ(i, j) i, j=1 in terms of the quantities
determining the probabilistic structure of the chain Z.
In particular, we have

Λ0(i, j) = Cov 1{Z1 =θ̃i } , 1{Z1 =θ̃ j } , i, j = 1, . . . , N 2 , (13)

∂ 2 h (νν )
∂ 2 h (νν )
1
.
=
=− N
∂ νi j ∂ νkl
∂ νi j ∂ νil
∑ j=1 νi j
π ), which will play
All these terms evaluated at π yield ∇2 h (π
Γn kπ
π ).
a crucial role in approximating D (Γ
a.s.
From the ergodicity of the chain Z, it follows Γ n −−−→ π .
n→∞

and
j)
Λ(i,
= Cov 1{Z1 =θ̃i } , 1{Z1+m =θ̃ j }
m




+ Cov 1{Z1 =θ̃ j } , 1{Z1+m =θ̃i } ,

(14)

m = 1, 2, . . . ; i, j = 1, . . . , N 2 .
We first determine (13) as follows. By direct calculation,
we obtain Λ(1,1)
= π̃1 (1 − π̃1 ), Λ(1,2)
= −π̃1 π̃2 , and so on. In
0
0
j)
general, we have Λ(i,
=
π̃
(I
−
π̃
)
, i, j = 1, . . . , N 2 , where
i
i
j
j
0
Ii j denotes the (i, j)-th entry of the identity matrix.
j)
Now we compute Λ(i,
2 by (14). Omitting the details, again

by direct calculation, we arrive at Cov 1{Z1 =θ̃i } , 1{Z3 =θ̃ j } =

π̃i P2i j − π̃ j , where P2i j is the (i, j)-th entry of the matrix P2
(the square of the transition matrix P). Interchanging
the

indexes i, j, we obtain Cov 1{Z1 =θ̃ j } , 1{Z3 =θ̃i } = π̃ j P2ji −



j)
2
2
π̃i . Thus, we have Λ(i,
2 = π̃i Pi j − π̃ j + π̃ j P ji − π̃i .
Computing in a similar way, in general, we derive


j)
m
m
Λ(i,
m = 1, 2, . . . ,
m = π̃i Pi j − π̃ j + π̃ j P ji − π̃i ,
m
where Pm
i j is the (i, j)-th entry of the matrix P (the m-th
power of the transition matrix P).
Thus, for i, j = 1, . . . , N 2 we finally attain



∞
m
Λ(i, j) = π̃i (Ii j − π̃ j ) + ∑m=1 π̃i Pm
i j − π̃ j + π̃ j P ji − π̃i .
(15)

It is seen that both ∇h (νν ) and ∇2 h (νν ) are continuous in a
neighborhood of π , and we can use the second-order Taylor’s
series expansion of h (νν ) centered at π to approximate
Γn kπ
π ) = h(Γ
Γn ) − h (π
π ). To be specific, we have
D (Γ
1
Γn kπ
π ) ≈ (Γ
Γn − π )0 ∇2 h (π
π ) (Γ
Γn − π ) .
D (Γ
2
C. Threshold approximation

We use an empirical CDF to approximate the actual CDF
of D ( Γ n k π ). In particular, by (19) we have
0

√
1 √
Γn kπ
Γn − π ) ∇2 h (π
π ) n (Γ
Γn − π ) ,
π)≈
D (Γ
n (Γ
2n
(20)
√
d.
Γn − π ) −−−→ N (0, Λ ) with Λ
where, by Thm. III.3, n (Γ
n→∞
given by (15). Thus, to derive an empirical CDF, we can
generate a set of Gaussian sample vectors independently
according to N (0, Λ ) and √
then plug each such sample
Γn − π )) and compute the
vector into (20) (i.e., replace n (Γ
right hand side of (20), thus, obtaining a set of sample
Γn kπ
π ).
scalars, as an approximation for samples of D (Γ
Γn kπ
π ), say,
Once we obtain an empirical CDF of D (Γ
denoted Fem (·; n), then, by (5), we can use it to estimate
ηn as
−1
(1 − βn ; n) ,
ηnwc ≈ Fem

Γn kπ
π)
B. Taylor’s series expansion of D (Γ
Inspired by [5] (or, see [8]) wherein the i.i.d. assumption is
required, we use a Taylor’s series expansion to approximate
Γn kπ
π ).
the relative entropy D (Γ
To this end, for ν ∈ P (Θ) let us consider
N

N

π ) = ∑ ∑ νi j log
h (νν ) = D (νν kπ
i=1 j=1

νi j
N ν
∑t=1
it
πi j
N π
∑t=1
it

.

By direct calculation, we derive


∂ h (νν )
N
= log νi j − log ∑t=1 νit − log πi j
∂ νi j


N
+ log ∑t=1 πit , i, j = 1, . . . , N,

(16)

(21)

−1 (·;
Fem

where
n) is the inverse of Fem (·; n).
Clearly, to calculate the right hand side of (20), we need
the parameter π . In practice, we typically do not know the
actual values of the underlying chain parameters, so we need
to estimate them accordingly. This can be done by computing
Γ n over a long past sample path.
IV. E XPERIMENTAL R ESULTS
A. Numerical results for threshold approximation

(17)

which leads to
π ) = 0.
∇h (π

(19)

(18)

The experiments in this subsection are performed using
the software package TAHTMA [17]. For convenience, we
use Θ = {1, 2, . . . , N 2 } to indicate the states, and assume π
(ground truth) is the initial distribution.
In the following numerical examples, we first randomly
create a valid (i.e., such that (C) holds) N × N transition
matrix Q, hence an N 2 × N 2 transition matrix P, and then

generate T samples
of the chain Z, each with length n,

denoted Z(t) = Z1(t) , . . . , Zn(t) , t = 1, . . . , T . These samples
will be used to derive empirical CDF’s. Also, to estimate the
parameters
 of interest, we generate one more training sample
Z(0) = Z1(0) , . . . , Zn(0)0 , where n0 is the length. The stationary
distribution π (estimation) of the chain is computed by taking
any row of Pm0 , where m0 is a sufficiently large integer.
With samples Z(t) = Z1(t) , . . . , Zn(t) , t = 1, . . . , T , we can
Γn kπ
π ),
compute T samples of the scalar random variable D (Γ
Γn kπ
π ), denoted
by (2). Thus, an empirical CDF of D (Γ
F (·; n), can be derived. We will treat F (·; n) as the actual
Γn kπ
π ). The threshold given by (21) with F (·; n)
CDF of D (Γ
(i.e., F(·; n) plays the role of Fem (·; n)) is then taken as the
theoretical (actual) value ηn .
Γn kπ
π ),
Next we derive an estimated empirical CDF of D (Γ
using Taylor’s series expansion together with Thm. III.3.
Let us first estimate the parameters of interest.
Recall π =

(πi j ; i, j = 1, . . . , N) = π̃k ; k = 1, . . . , N 2 . We estimate π̃k
by

 n0 

 ∑i=1 1 Z (0) =k
i
, ε , k = 1, . . . , N 2 , (22)
π̃ˆ k = max


n0
where ε > 0 is a small number. The purpose of introducing
ε is to avoid division
 2 by zero. Then, we obtain an estimated
π = π̃ˆ k ∑Nj=1 π̃ˆ j ; k = 1, . . . , N 2 , or written π̂
π =
π as π̂


(π̂i j ; i, j = 1, . . . , N), where π̂11 = π̃ˆ 1 ŝ, . . ., π̂1N = π̃ˆ N ŝ,


2
. . ., π̂N1 = π̃ˆ (N−1)N+1 ŝ, . . ., π̂NN = π̃ˆ N 2 ŝ, with ŝ = ∑Nj=1 π̃ˆ j
being the normalization constant.
π ), by
We are now ready to calculate the estimate for ∇2 (π
2
π
evaluating ∇ (π̂ ).
Further, by Lemma III.1, we estimate Q by
q̂i j =

π̂i j
N
π̂it
∑t=1

,

i, j = 1, . . . , N.

(23)

Then an estimation of P, denoted P̂, can be directly derived
from Q̂.
Finally, by (15), we estimate the covariance matrix Λ using


m0  ˆ
ˆ
Λ̂(i, j) = π̃ˆ i Ii j − π̃ˆ j + ∑m=1
π̃ i P̂m
i j − π̃ j

+ π̃ˆ j P̂mji − π̃ˆ i , i, j = 1, . . . , N 2 .
(24)
It is worth noting that, to ensure theestimated Λ is at least
0
Λ by (Λ̂
Λ + Λ̂
Λ ) 2. Also, to ensure Λ̂
Λ to
symmetric, we update Λ̂
be positive semi-definite, we use the trick of QR factorization
[18]; for details, the reader is referred to [17].
We now generate T Gaussian sample vectors according to
Λ) and then plug each of them into (20) and compute
N (0, Λ̂
π ) replaced by ∇2 (π̂
π ),
the right hand side of (20) with ∇2 (π
thus, obtaining T scalar samples as approximations for samΓn kπ
π ). Therefore, an estimated empirical CDF of
ples of D (Γ
Γn kπ
π ), denoted Fem (·; n), can be derived accordingly.
D (Γ
ηnwc given by (21) is then the WC approximation for the
threshold.
Let the false alarm probability be βn = β = 0.001. We set
the aforementioned parameters as N = 12, ε = 10−8 , T =
1000, m0 = 1000, and let n0 = 1000N 2 . In Fig. 1, the red line

indicates the theoretical value of the threshold ηn , the blue
line shows the WC approximation ηnwc , and the green line
demonstrates the estimate ηnsv obtained by Sanov’s theorem,
all as a function of the number of samples n. It is seen that
ηnwc is much more accurate than ηnsv .

Fig. 1. Threshold (η) versus number of samples (n); the number of states
in the original Markov chain is N = 12.

Remark 3: In Fig. 1, the red line corresponding to the
actual value ηn is not as smooth; the reason is that each time
when varying the
 number of samples n, we regenerate all the
samples Z(t) = Z1(t) , . . . , Zn(t) , t = 1, . . . , T , from scratch. On
the other hand, the blue line corresponding to ηnwc looks
smooth; this is because we only need to generate the T
Λ) once. Note
Gaussian sample vectors according to N (0, Λ̂
that, in our experiments, most 
of the running time is spent
generating the samples Z(t) = Z1(t) , . . . , Zn(t) , t = 1, . . . , T ,
and then calculating ηn . In practice, when implementing
the approximation procedure proposed in this subsection,
we will only need to focus on obtaining ηnwc , which is
computationally inexpensive.
B. Simulation results for network anomaly detection
We will use the term traffic and flow interchangeably in
this subsection; they mean the same thing. The simulations
are done using the software package SADIT [11], which,
based on the fs-simulator [19], can efficiently generate flowlevel network traffic datasets with annotated anomalies.
The network (see Fig. 2) consists of an internal network
involving 8 normal users (CT1-CT8) and 1 server (SRV)
that stores some sensitive information, and 3 Internet nodes
(INT1-INT3) that visit the internal network via a gateway
(GATEWAY).
When dealing with the data, we use the features “flow
duration,” “flow size,” and “distance to cluster center” (this
is related to a user IP address). The data preprocessing
involves three steps: i) network traffic representation and flow
aggregation, ii) quantization, and iii) window aggregation;
cf. [4, Sec. III.A].
To implement Hoeffding’s test, we first estimate a PL
(resp., a set of PLs) for the stationary (resp., time-varying)
normal traffic. For the former case, we use as the reference
traffic the whole flow sequence with anomalies at some

Fig. 2.

Simulation setting (this figure is from [4] (or, see [20])).

time interval that we seek to detect. For the latter case, we
generate the reference traffic and the traffic with anomalies
separately, sharing all of the parameters except the ones
regarding the anomalies. Using the windowing technique
[20], [4], we then persistently monitor the traffic and report
anomalies instantly whenever the relative entropy is higher
than the threshold corresponding to the detection window.
The key difference between the whole detection procedures for the two types of traffic is that, estimating a PL
for the stationary traffic is straightforward, while, for the
time-varying traffic, we need to make an effort to identify
several different PLs corresponding to certain time intervals.
We apply the two-steps procedure proposed in [4]; that is,
with some rough estimation of the shifting and periodic
parameters of the traffic, we first generate a large set of
PLs, and then refine the large family of PLs by solving a
weighted set cover problem. For further theoretical details
and the implementation, the reader is referred to [4], [11].
Note that, to deal with the time-varying traffic, [4] proposes a generalized Hoeffding’s test, which we will use for
Scenario 2 in the following. However, there is no essential
difference compared to the typical Hoeffding’s test. Note also
that, in the following experiments, we only use the modelbased method [20], [4].
1) Scenario 1 (Stationary Network Traffic): We mimic the
case for anomaly caused by “large access rate;” cf. [20, Sec.
IV.A.3]. The simulation time is 7000 s. A user suspiciously
increases its access rate to 10 times of its normal value
between 4000 s and 4500 s. The interval between the starting
point of two consecutive time windows is taken as 50 s, the
window-size is chosen to be ws = 400 s, and the false alarm
probability is set as β = 0.001.
Set the quantization level for “flow duration,” “flow size,”
and “distance to cluster center” to be n1 = 1, n2 = 2, and
n3 = 2, respectively. Set the number of user clusters as k = 3.
See Fig. 3 and Fig. 4 for the detection results. Both figures
depict the relative entropy (divergence) metric in (2). The
green dashed line in Fig. 3 is the threshold calculated using
Sanov’s theorem (i.e., ηnsv , where n is the number of flows in
a specific detection window). The green dashed line in Fig.
4 is the threshold computed using the WC result derived in
Sec. III (i.e., ηnwc ). The interval during which the entropy
curve is above the threshold line (the red part) is the interval

Fig. 3. Detection result for Scenario 1 with ws = 400 s, n1 = 1, n2 = 2,
n3 = 2, k = 3; threshold is estimated by use of Sanov’s theorem.

Fig. 4. Detection result for Scenario 1 with ws = 400 s, n1 = 1, n2 = 2,
n3 = 2, k = 3; threshold is estimated by use of weak convergence result.

reported as abnormal.
Fig. 3 shows that, if using ηnsv as the threshold, then
the detection method reports an anomaly for all windows
including the ones wherein the behavior of the network traffic
is actually normal; in other words, there are too many false
alarms. Clearly, this is undesirable. Fig. 4 shows that, if,
instead, we use ηnwc as the threshold, then the detection
method will not report any false alarm while successfully
identifying the true anomalies between 4000 s and 4500 s.
The reason for such detection results is that ηnwc is more
accurate than ηnsv .
2) Scenario 2 (Time-Varying Network Traffic): Consider
the network with day-night pattern where the flow size
follows a log-normal distribution. We use exactly the same
scenario as that in [4, Sec. IV.B.2].
Applying the same procedure as in [4, Sec. III.C], we first
obtain 32 rough PLs. Then, using the heuristic PL refinement
algorithm given in [4, Sec. III.D] equipped with the up-bound
nominal cross-entropy parameter λ = 0.027565, which is
determined by the WC approximation, we finally obtain 6
PLs (see Fig. 5). These PLs are active during morning,
afternoon, evening, mid-night, dawn, and the transition time
around sunrise, respectively. In the following detection pro-

cedure, the key difference between our method and the one
used in [4] is that we no longer manually set the threshold
universally as a constant; instead, the threshold ηnwc for each
detection window is automatically calculated by use of the
WC approximation. We set the quantization level for “flow
duration,” “flow size,” and “distance to cluster center” to be
n1 = 1, n2 = 4, and n3 = 1, respectively, and set the number
of user clusters to k = 1. The interval between the starting
point of two consecutive time windows is chosen as 1000 s,
the window-size is taken as ws = 1000 s, and the false alarm
probability is set as β = 0.001.
Noting that the anomaly is simulated beginning at 59 h
and lasting for 80 minutes [4], we see from Fig. 6 that the
anomaly is successfully detected, without any false alarms.

Fig. 5. PLs identification for Scenario 2 with ws = 1000 s, n1 = 1, n2 = 4,
n3 = 1, k = 1.

Fig. 6. Detection result for Scenario 2 with ws = 1000 s, n1 = 1, n2 = 4,
n3 = 1, k = 1.

V. C ONCLUSIONS
We establish a weak convergence result for the relative
entropy in Hoeffding’s test under a Markovian assumption,
which enables us to obtain a more accurate approximation
(compared to the existing approximation given by Sanov’s
theorem) for the threshold needed by the test. We validate
the accuracy of such approximation through numerical experiments and simulations for network anomaly detection.

ACKNOWLEDGEMENTS
The authors would like to thank Dr. Jing Wang for his
contributions and help in developing the software package
SADIT [11].
R EFERENCES
[1] A. Dembo and O. Zeitouni, Large deviations techniques and applications. Springer, 1998.
[2] I. C. Paschalidis and G. Smaragdakis, “Spatio-temporal network
anomaly detection by assessing deviations of empirical measures,”
IEEE/ACM Transactions on Networking (TON), vol. 17, no. 3, pp.
685–697, 2009.
[3] S. Meyn, A. Surana, Y. Lin, and S. Narayanan, “Anomaly detection
using projective Markov models in a distributed sensor network,” in
Proceedings of the 48th IEEE Conference on Decision and Control,
2009 held jointly with the 2009 28th Chinese Control Conference.
CDC/CCC 2009., Dec 2009, pp. 4662–4669.
[4] J. Wang and I. Paschalidis, “Statistical traffic anomaly detection in
time-varying communication networks,” IEEE Transactions on Control
of Network Systems, vol. 2, no. 2, pp. 100–111, 2015.
[5] J. Unnikrishnan and D. Huang, “Weak convergence analysis of asymptotically optimal hypothesis tests,” IEEE Transactions on Information
Theory, submitted, 2013.
[6] W. Hoeffding, “Asymptotically optimal tests for multinomial distributions,” The Annals of Mathematical Statistics, pp. 369–401, 1965.
[7] J. Zhang, “Threshold Approximation for Hoeffding’s Test under I.I.D.
Assumption (TAHTIID),” https:// github.com/ jingzbu/ TAHTIID, 2015.
[8] J. Unnikrishnan, D. Huang, S. P. Meyn, A. Surana, and V. V. Veeravalli, “Universal and composite hypothesis testing via mismatched
divergence,” IEEE Transactions on Information Theory, vol. 57, no. 3,
pp. 1587–1603, 2011.
[9] S. S. Wilks, “The large-sample distribution of the likelihood ratio for
testing composite hypotheses,” The Annals of Mathematical Statistics,
vol. 9, no. 1, pp. 60–62, 1938.
[10] M. Iltis, “Sharp asymptotics of large deviations in Rd ,” Journal of
Theoretical Probability, vol. 8, no. 3, pp. 501–522, 1995.
[11] J. Wang, J. Zhang, and I. C. Paschalidis, “Systematic Anomaly Detection of Internet Traffic (SADIT),” https:// github.com/ hbhzwj/ SADIT,
2014.
[12] G. L. Jones, “On the Markov chain central limit theorem,” Probability
Surveys, vol. 1, pp. 299–320, 2004.
[13] S. P. Meyn and R. L. Tweedie, Markov Chains and Stochastic Stability.
New York: Springer, 1993.
[14] G. O. Roberts and J. S. Rosenthal, “General state space Markov chains
and MCMC algorithms,” Probability Surveys, vol. 1, pp. 20–71, 2004.
[15] F. Pène, “Rate of convergence in the multidimensional central limit
theorem for stationary processes. Application to the Knudsen gas and
to the Sinai billiard,” The Annals of Applied Probability, vol. 15, no. 4,
pp. 2331–2392, 2005.
[16] W. Feller, An introduction to probability theory and its applications.
John Wiley, 1957, vol. 2.
[17] J. Zhang, “Threshold Approximation for Hoeffding’s Test under Markovian Assumption (TAHTMA),” https:// github.com/ jingzbu/
TAHTMA, 2015.
[18] R. A. Horn and C. R. Johnson, Matrix analysis. Cambridge university
press, 2012.
[19] J. Sommers, R. Bowden, B. Eriksson, P. Barford, M. Roughan,
and N. Duffield, “Efficient network-wide flow record generation,” in
INFOCOM, 2011 Proceedings IEEE. IEEE, 2011, pp. 2363–2371.
[20] J. Wang, D. Rossell, C. G. Cassandras, and I. C. Paschalidis, “Network
anomaly detection: A survey and comparative analysis of stochastic
and deterministic methods,” in IEEE 52nd Annual Conference on
Decision and Control (CDC), Dec 2013, pp. 182–187.

