# A Fast Framework for Abrupt Change Detection Based on Binary Search Trees and Kolmogorov Statistic

Authors: Jin-Peng Qi; Jie Qi; Qing Zhang
Journal: Computational Intelligence and Neuroscience
 (2016)

---

## Abstract

Change-Point (CP) detection has attracted considerable attention in the fields of data mining and statistics; it is very meaningful to discuss how to quickly and efficiently detect abrupt change from large-scale bioelectric signals. Currently, most of the existing methods, like Kolmogorov-Smirnov (KS) statistic and so forth, are time-consuming, especially for large-scale datasets. In this paper, we propose a fast framework for abrupt change detection based on binary search trees (BSTs) and a modified KS statistic, named BSTKS (binary search trees and Kolmogorov statistic). In this method, first, two binary search trees, termed as BSTcA and BSTcD, are constructed by multilevel Haar Wavelet Transform (HWT); second, three search criteria are introduced in terms of the statistic and variance fluctuations in the diagnosed time series; last, an optimal search path is detected from the root to leaf nodes of two BSTs. The studies on both the synthetic time series samples and the real electroencephalograph (EEG) recordings indicate that the proposed BSTKS can detect abrupt change more quickly and efficiently than KS,t-statistic (t), and Singular-Spectrum Analyses (SSA) methods, with the shortest computation time, the highest hit rate, the smallest error, and the highest accuracy out of four methods. This study suggests that the proposed BSTKS is very helpful for useful information inspection on all kinds of bioelectric time series signals.


---

## Body

## 1. Introduction

Abrupt change detection is to identify abrupt changes in the statistical properties of a signal series, which occur at unknown instants [1–3]. These changes are interesting because they are indicative of qualitative transitions in the data generation mechanism (DGM) underlying the signals. Currently, CP detection has attracted considerable attention in the fields of data mining and statistics, and it has been widely studied in many real-world problems, such as atmospheric and financial analyses [1], fault detection in engineering system [4, 5], climate change detection [6], genetic time series analyses [7], signal segmentation [8, 9], and intrusion detection in computer network [4].In community of statistics, some nonparametric approaches for CP detection have been widely explored. For example, KS statistic quantifies a distance between the empirical distribution function of the sample and the cumulative distribution function of the reference distribution or between the empirical distribution function of two samples [10, 11]. Also, KS statistic and its modified versions are broadly investigated on many application fields, for example, testing hypotheses regarding activation in blood oxygenation level-dependent functional MRI data [12], modeling the cumulative distribution function of rub-induced AE signals, quantifying the goodness of fit to offer a suitable signal feature for diagnosis [13], as well as abrupt change detecting on EEG signals [14], and gene expression time series [15]. Meanwhile, as for the model-related statistic approaches, some modified cumulative sum (CUSUM) methods provide the asymptotic distributions of test statistics and the consistency of procedures and behave better in finite samples and have a higher stability with respect to the time of change than ordinary CUSUM procedures [16]. The CUSUM method and its revised versions have been widely applied to detect the structural breaks in the parameters of stochastic models, as well as the abrupt changes in the regression parameters of multiple time series regression models, such as multiple CP detection in biological sequences [17], abrupt change detection in the regression parameters of a set of capital asset pricing data related to the Fama-French extension of the CAPM [16], and abrupt change detection in a shape-restricted regression model [18].On the other hand, SSA is a powerful technique for time series analyses. SSA is nonparametric and requires no prior knowledge on the properties of time series signal [19]. The main idea of SSA is applied in the principal component analyses on the trajectory matrix with subsequent reconstruction of the original time series. SSA has been proved to be very successful and has already become a standard tool in the analyses of climatic [10], meteorological, and geophysical time series [11, 19]. Currently, SSA has been successfully applied in the real time series recordings, for example, abrupt change analyses on EMG-onset detection [12] and CP detection in time series [13]. Although SSA is a model-free method, it is not scalable to large-scale datasets, because it is time-consuming and sometimes invalid for time series analyses with less significant data fluctuation.In addition, Wavelet Transform (WT) is another important tool for time series analyses [14, 15, 20–23]. WT has been widely applied in anomaly detection, time series prediction, image processing, and noise reduction [15, 23–25]. WT can represent general function at different scales and positions in a versatile and sophisticated manner, so that the data distribution features can be easily extracted from different time or space scales [25, 26]. As a simple WT, Haar Wavelet (HW) owns some attractive features including fast implementation and ability to analyze the local features. HW is very useful to find abrupt changes of discontinuity and high frequency in time series, so it is a potential candidate in modern electrical and computer engineering applications, such as signal and image compression, eye detection [27], abnormality detection on time series [28, 29], and abrupt change detection on autoregressive conditional heteroscedastic processes [30].However, all of these methods above are time-consuming and sometime invalid for abrupt change detection near the left or the right boundary, especially for insignificant data fluctuation in large-scale time series. To resolve these problems, we propose a fast framework for CP detection based on binary search trees and a modified KS statistic, termed BSTKS for short. In this novel method, first, two BSTs are derived from a diagnosed time series. Second, three search criteria are introduced in terms of the statistic and variance fluctuations between two adjacent time series segments, and then an optimal search path is detected from the root to leaf nodes of two BSTs. Last, the proposed BSTKS and other KS,t, and SSA methods are tested on both the synthetic time series and real EEG recordings and evaluated in terms of computation time, hit rate, error, accuracy, and area under curve (AUC) of Receiver Operating Characteristic (ROC) curve analyses.In general, for a certain bioelectric signal, an abrupt change means an important transition of biological functions or health states before and after a strong attack or an acute perturbation from internal or external environment. Therefore, it is very necessary to not only discern abrupt change from all kinds of physiological and psychological time series signals, but also inspect the significant fluctuation between adjacent time series segments with different scales. The following sections focused on not only presenting the framework of the proposed BSTKS method through theoretical foundation, simulation, and evaluation, but also discussing how it can more quickly and efficiently detect abrupt change on both synthetic and real bioelectric EEG signals than other existing KS,t, and SSA methods. The rest of this paper is organized as follows. Section 2 gives the preliminary of abrupt change by introducing the statistic and variance fluctuations between two adjacent time series segments. Section 3 implements the integrated framework of the BSTKS method in terms of three search criteria in detail. Section 4 provides some representative experiments by using the synthetic time series and real EEG recordings and then analyzes the performance of BSTKS by comparing with other KS, t, and SSA methods. Section 5 gives summary and conclusion from previous sections.


## 2. Preliminary

### 2.1. Statistic Fluctuation

KS statistic is sensitive to differences in both location and shape of the cumulative distribution functions (c.d.f) of two samples. The null distribution of KS statistic is calculated under the null hypothesis that the two samples are drawn from the same distribution or one sample is drawn from the reference distribution. To detect an abrupt change from a diagnosed time seriesZ, we define the statistic fluctuation between two adjacent segments within Z by means of KS statistic as follows [1, 4, 19].Definition 1.
Supposing a time series sample,Z
=
{
z
1
,
…
,
z
N
}, one observes (1)
Z
=
f
i
n
+
X
,
i
=
1
,
…
,
N
,where X
=
{
x
i
}
i
=
1
,
…
,
N is a set of the discrete and centred i.i.d random variables and f is a noisy mean signal with unknown distribution. The statistic fluctuation between two adjacent segments Z
L
=
{
z
a
,
…
,
z
c
} and Z
R
=
{
z
c
+
1
,
…
,
z
b
} is defined as(2)
S
m
n
x
=
m
n
m
+
n
1
/
2
sup
x
∈
R
⁡
F
m
x
-
G
n
x
,in which F
m
(
x
) and G
n
(
x
) are the c.d.f of Z
L and Z
R, respectively; m
=
c
-
a, n
=
b
-
c
-
1, and m
+
n
≤
N. Supposing the hypothesized F
m
(
x
) and G
n
(
x
) in (2) are not available, we can derive the empirical cumulative distribution functions (e.c.d.f) of F
m
(
x
) and G
n
(
x
) from Z
L and Z
R. Then, F
m
(
x
) and G
n
(
x
) can be redefined as(3)
F
m
x
=
P
m
Z
L
≤
x
=
1
m
∑
i
=
a
c
I
z
i
≤
x
,
G
n
x
=
P
n
Z
R
≤
x
=
1
n
∑
j
=
c
+
1
N
I
z
j
≤
x
,where F
m
(
x
) and G
n
(
x
) count the proportion of the sample points below level x.Hypothesis 1.
In order to discern an abrupt change onZ in terms of statistic fluctuation defined above, we introduce KS test for two adjacent segments Z
L and Z
R in Z as(
H
0
)
ifS
m
n
(
z
c
)
≤
δ, no abrupt change occurs in Z;
(
H
1
)
ifS
m
n
(
z
c
)
>
δ, abrupt change occurs in Z,
in whichδ
∈
R is a threshold of the statistic fluctuation within Z belonging to an identical distribution. Then, we test (
H
0
) against (
H
1
) from observations. If an abrupt change c occurs in Z, there exists a value c satisfying S
m
n
(
z
c
)
>
δ, z
c
∈
z
1
,
z
N, and δ
∈
R. In this hypothesis, we assume that the number, the location, and the size of the function f in (1) are unknown, and the upper bound of the statistic fluctuation δ is supposed to be known.


### 2.2. Variance Fluctuation

Provided the statistic fluctuation defined in (2) is insignificant enough, it is difficult to detect abrupt change near the left or the right boundary within Z, especially when sample size N gets smaller. Therefore, we need to introduce another variable to calculate the variance fluctuation between two adjacent parts within a time series sample.Definition 2.
Supposing two adjacent segmentsZ
L
=
{
z
a
,
…
,
z
c
} and Z
R
=
z
c
+
1
,
…
,
z
b in Z
=
{
z
1
,
…
,
z
N
}, the variance fluctuation between Z
L and Z
R is defined as(4)
D
m
n
c
=
s
u
p
1
≤
L
,
R
≤
N
1
m
∑
L
=
a
c
z
L
-
1
n
∑
R
=
c
+
1
b
z
R
,in which m
=
c
-
a, n
=
b
-
c
-
1, and m
+
n
≤
N.Hypothesis 2.
(
H
0
) If D
m
n
(
c
)
≤
β, no abrupt change occurs at c in Z; (
H
1
) if D
m
n
(
c
)
>
β, abrupt change occurs at c in Z.
Here,β
∈
R is a variance threshold of time series Z which obeys an identical distribution. If there exists a value c satisfying D
m
n
(
c
)
>
β
,
z
c
∈
[
z
1
,
z
N
], then an abrupt change occurs at c in Z.


## 2.1. Statistic Fluctuation

KS statistic is sensitive to differences in both location and shape of the cumulative distribution functions (c.d.f) of two samples. The null distribution of KS statistic is calculated under the null hypothesis that the two samples are drawn from the same distribution or one sample is drawn from the reference distribution. To detect an abrupt change from a diagnosed time seriesZ, we define the statistic fluctuation between two adjacent segments within Z by means of KS statistic as follows [1, 4, 19].Definition 1.
Supposing a time series sample,Z
=
{
z
1
,
…
,
z
N
}, one observes (1)
Z
=
f
i
n
+
X
,
i
=
1
,
…
,
N
,where X
=
{
x
i
}
i
=
1
,
…
,
N is a set of the discrete and centred i.i.d random variables and f is a noisy mean signal with unknown distribution. The statistic fluctuation between two adjacent segments Z
L
=
{
z
a
,
…
,
z
c
} and Z
R
=
{
z
c
+
1
,
…
,
z
b
} is defined as(2)
S
m
n
x
=
m
n
m
+
n
1
/
2
sup
x
∈
R
⁡
F
m
x
-
G
n
x
,in which F
m
(
x
) and G
n
(
x
) are the c.d.f of Z
L and Z
R, respectively; m
=
c
-
a, n
=
b
-
c
-
1, and m
+
n
≤
N. Supposing the hypothesized F
m
(
x
) and G
n
(
x
) in (2) are not available, we can derive the empirical cumulative distribution functions (e.c.d.f) of F
m
(
x
) and G
n
(
x
) from Z
L and Z
R. Then, F
m
(
x
) and G
n
(
x
) can be redefined as(3)
F
m
x
=
P
m
Z
L
≤
x
=
1
m
∑
i
=
a
c
I
z
i
≤
x
,
G
n
x
=
P
n
Z
R
≤
x
=
1
n
∑
j
=
c
+
1
N
I
z
j
≤
x
,where F
m
(
x
) and G
n
(
x
) count the proportion of the sample points below level x.Hypothesis 1.
In order to discern an abrupt change onZ in terms of statistic fluctuation defined above, we introduce KS test for two adjacent segments Z
L and Z
R in Z as(
H
0
)
ifS
m
n
(
z
c
)
≤
δ, no abrupt change occurs in Z;
(
H
1
)
ifS
m
n
(
z
c
)
>
δ, abrupt change occurs in Z,
in whichδ
∈
R is a threshold of the statistic fluctuation within Z belonging to an identical distribution. Then, we test (
H
0
) against (
H
1
) from observations. If an abrupt change c occurs in Z, there exists a value c satisfying S
m
n
(
z
c
)
>
δ, z
c
∈
z
1
,
z
N, and δ
∈
R. In this hypothesis, we assume that the number, the location, and the size of the function f in (1) are unknown, and the upper bound of the statistic fluctuation δ is supposed to be known.


## 2.2. Variance Fluctuation

Provided the statistic fluctuation defined in (2) is insignificant enough, it is difficult to detect abrupt change near the left or the right boundary within Z, especially when sample size N gets smaller. Therefore, we need to introduce another variable to calculate the variance fluctuation between two adjacent parts within a time series sample.Definition 2.
Supposing two adjacent segmentsZ
L
=
{
z
a
,
…
,
z
c
} and Z
R
=
z
c
+
1
,
…
,
z
b in Z
=
{
z
1
,
…
,
z
N
}, the variance fluctuation between Z
L and Z
R is defined as(4)
D
m
n
c
=
s
u
p
1
≤
L
,
R
≤
N
1
m
∑
L
=
a
c
z
L
-
1
n
∑
R
=
c
+
1
b
z
R
,in which m
=
c
-
a, n
=
b
-
c
-
1, and m
+
n
≤
N.Hypothesis 2.
(
H
0
) If D
m
n
(
c
)
≤
β, no abrupt change occurs at c in Z; (
H
1
) if D
m
n
(
c
)
>
β, abrupt change occurs at c in Z.
Here,β
∈
R is a variance threshold of time series Z which obeys an identical distribution. If there exists a value c satisfying D
m
n
(
c
)
>
β
,
z
c
∈
[
z
1
,
z
N
], then an abrupt change occurs at c in Z.


## 3. Method

### 3.1. Two BSTs’ Construction

In the first part of the proposed BSTKS method, two BSTs, that is, BSTcA and BSTcD, are constructed from a time series sampleZ, by using multilevel HWT. Generally, as shown in Figure 1, a discrete time series signal Z
=
{
z
1
,
z
2
,
…
,
z
N
} can be decomposed into the kth-level trend c
A
k and k-level fluctuations, that is, c
D
1
,
c
D
2
,
…
,
c
D
k, k
=
1,2
,
…
,
l
o
g
2
N. The k-level HWT is the mapping H
k defined as [13](5)
Z
→
H
k
c
A
k
∣
c
D
k
∣
c
D
k
-
1
∣
⋯
∣
c
D
2
∣
c
D
1
,and then, the mapping H
k can be represented by the approximation and detail coefficient matrices, termed McA and McD as follows: (6)
M
c
A
=
c
A
1,1
⋯
c
A
1
,
N
⋮
c
A
k
,
j
⋮
c
A
M
,
1
0
0
,
M
c
D
=
c
D
1,1
⋯
c
D
1
,
N
⋮
c
D
k
,
j
⋮
c
D
M
,
1
0
0
,where 0
≤
k
≤
M
=
l
o
g
2
N and 1
≤
j
≤
N
/
2
k.Figure 1
The diagram of a discrete time seriesZ decomposition by k-level HWT, which is composed of k-level cA and cD, that is, the average and difference coefficient vectors.Supposing the size of a diagnosedZ is divisible k times by 2, the jth element c
A
k
,
j in c
A
k and the jth element c
D
k
,
j in c
D
k can be denoted as(7)
c
A
k
,
j
=
1
2
∧
k
∑
i
=
a
b
z
i
,
c
D
k
,
j
=
1
2
∧
k
∑
L
=
a
c
z
L
-
∑
R
=
c
+
1
b
z
R
,where 1
≤
k
≤
l
o
g
2
N and 2
k
(
j
-
1
)
+
1
≤
i
≤
j
∗
2
k; a
=
2
k
(
j
-
1
)
+
1, c
=
2
k
(
j
-
1
)
+
2
(
k
-
1
), and b
=
2
k
∗
j.During two BSTs’ construction, as shown in Figure2, the non-leaf nodes in BSTcA and BSTcD are assembled by the k-level coefficient vectors of McA and McD, respectively; and then the leaf nodes are derived directly from the original time series Z. Therefore, the features of abrupt change in Z can be reflected and distributed into the different non-leaf nodes of BSTcA and BSTcD, in accordance with the k level coefficient vectors in McA and McD.Figure 2
The diagrams of two binary trees, BSTcA and BSTcD, which are constructed by McA and McD, as well as the original time seriesZ.
(a)
(b)


### 3.2. CP Detection Based on Three Search Criteria

To find an optimal path towards the potential CP within a given time seriesZ quickly and efficiently, some search criteria need to be introduced, and then the data exceptions can be detected from the root to leaf nodes of two BSTs. As for the statistic fluctuation within BSTcA, first, a new variable z
k
,
j is defined according to a current non-leaf node c
A
k
,
j in BSTcA,(8)
z
k
,
j
=
1
2
k
∑
i
=
a
b
z
i
=
2
∧
k
2
k
c
A
k
,
j
,where 1
≤
k
≤
l
o
g
2
N, 1
≤
j
≤
N
/
2
k; a
=
2
k
(
j
-
1
)
+
1, b
=
2
k
∗
j, and a
≤
i
≤
b. Then, the statistic fluctuation between two adjacent segments Z
L
=
{
z
a
,
…
,
z
c
} and Z
R
=
{
z
c
+
1
,
…
,
z
b
} can be defined by a modified KS statistic as(9)
S
m
n
′
k
,
j
=
n
m
n
+
m
1
/
2
1
n
∑
i
L
=
a
c
I
z
i
L
≤
z
k
,
j
-
1
m
∑
i
R
=
c
+
1
b
I
z
i
R
≤
z
k
,
j
,in which z
k
,
j is a new element defined in (8); m and n stand for the sizes of Z
L and Z
R, respectively; 1
≤
k
≤
l
o
g
2
N, 1
≤
j
≤
N
/
2
k
;
a
=
2
k
(
i
-
1
)
+
1, b
=
2
k
j, and c
=
2
k
(
j
-
1
)
+
2
(
k
-
1
). S
m
n
′
(
k
,
j
) measures the e.c.d.f difference between Z
L and Z
R, and the larger S
m
n
′
(
k
,
j
) means the more significant statistic fluctuation between Z
L and Z
R. Therefore, a potential abrupt change might occur at c in Z with more probability.Definition 3.
For a current non-leaf nodec
A
k
,
j in BSTcA, with its left and right-child nodes c
A
k
-
1,2
j
-
1 and c
A
k
-
1,2
j, the distance of e.c.d.f, S
k
,
j
;
L, and S
k
,
j
;
R can be defined as(10)
S
k
,
j
;
L
=
S
m
n
′
k
,
j
;
k
-
1,2
j
-
1
=
n
m
n
+
m
1
/
2
1
n
∑
i
=
a
b
I
z
i
≤
z
k
,
j
-
1
m
∑
i
L
=
a
c
I
z
i
L
≤
z
k
,
j
=
n
m
n
+
m
1
/
2
W
1
n
∑
i
=
a
b
I
z
i
≤
c
A
k
,
j
-
1
m
∑
i
L
=
a
c
I
z
i
L
≤
c
A
k
,
j
,
S
k
,
j
;
R
=
S
m
n
′
k
,
j
;
k
-
1,2
j
=
n
m
n
+
m
1
/
2
1
n
∑
i
=
a
b
I
z
i
≤
z
k
,
j
-
1
m
∑
i
R
=
c
+
1
b
I
z
i
R
≤
z
k
,
j
=
n
m
n
+
m
1
/
2
W
1
n
∑
i
=
a
b
I
z
i
≤
c
A
k
,
j
-
1
m
∑
i
R
=
c
+
1
b
I
z
i
R
≤
c
A
k
,
j
,where 2
≤
k
≤
l
o
g
2
N, 1
≤
j
≤
N
/
2
k; a
=
2
k
(
j
-
1
)
+
1, b
=
2
k
j, c
=
2
k
(
j
-
1
)
+
2
(
k
-
1
); n
=
2
k, m
=
2
k
-
1; and W
=
(
2
)
∧
k
/
2
k. To estimate an optimal path towards the potential change position within Z, without loss of generality, the first search criterion is introduced based on the statistic fluctuations S
k
,
j
;
L and S
k
,
j
;
R.Criterion 1.
Given two statistic fluctuation variablesS
k
,
j
;
L and S
k
,
j
;
R in accordance with two non-leaf child nodes c
A
k
-
1,2
j
-
1 and c
A
k
-
1,2
j of the current selected node c
A
k
,
j in BSTcA, and 2
≤
k
≤
l
o
g
2
N, (a)
if(
S
k
,
j
;
L
>
S
k
,
j
;
R
)
∧
(
S
k
,
j
;
L
>
C
(
α
)
) holds true, then the left-child node c
A
k
-
1,2
j
-
1 is selected and involved into the current search path; meanwhile, the right-child c
A
k
-
1,2
j is discarded;
(b)
if(
S
k
,
j
;
R
>
S
k
,
j
;
L
)
∧
(
S
k
,
j
;
R
>
C
(
α
)
) holds true, then the right-child node c
A
k
-
1,2
j is selected and involved into the current search path; meanwhile, the left-child c
A
k
-
1,2
j
-
1 is discarded.Proof.
For a selected non-leaf nodec
A
k
,
j in BSTcA, as shown in Figure 3, the original time series Z is divided equally into two adjacent segments Z
L and Z
R, which are covered by two non-leaf child nodes c
A
k
-
1,2
j
-
1 and c
A
k
-
1,2
j, respectively. According to the definitions of S
k
,
j
;
L and S
k
,
j
;
R in (10), the satisfied S
k
,
j
;
L
>
S
k
,
j
;
R indicates that the statistic fluctuation within Z
L is more significant than that one within Z
R; that is, a potential abrupt change might be contained in Z
L with more probability than in Z
R, and vice versa. Furthermore, if S
k
,
j
;
L
>
C
(
α
) holds true, then (
H
1
) of Hypothesis 1 is satisfied; that is, abrupt change occurs in Z
L, and vice versa, where C
(
α
) is the critical value predefined in an identical distribution and α is the significance level. Therefore, one of the two child nodes c
A
k
-
1,2
j
-
1 and c
A
k
-
1,2
j is selected and involved into the current search path; meanwhile, the remaining one is discarded. Once the statistic fluctuation is significant enough, an optimal search path can be detected by Criterion 1 from the top to the last non-leaf level in BSTcA. However, the search procedure is probably forced to cease because the statistic fluctuation is so insignificant that Criterion 1 is invalid for detecting it, especially for the left or the right boundary when sample Z is with smaller size N. Therefore, it is necessary to introduce another search criterion based on the variance fluctuations within BSTcD.Figure 3
The scheme of Criterion1 based on the statistic fluctuations within BSTcA. In terms of this criterion, the left or right-child node, that is, c
A
k
,
2
j
-
1 or c
A
k
,
2
j, might be selected to be involved in the current search path; meanwhile the remaining one is discarded. Thereafter, an optimal path towards the potential abrupt change in Z is expected to be obtained from BSTcA, after log
2
⁡
N binary search steps.Definition 4.
For a current non-leaf nodec
D
k
,
j in BSTcD, with its left and right-child nodes c
D
k
-
1,2
j
-
1 and c
D
k
-
1,2
j, respectively, the variance fluctuations D
k
,
j
;
L and D
k
,
j
;
R are defined in terms of (4) as(11)
D
k
,
j
;
L
=
D
m
n
′
k
,
j
;
k
-
1,2
j
-
1
=
n
m
n
+
m
1
/
2
1
n
∑
i
L
=
a
c
z
i
L
-
∑
i
R
=
c
+
1
b
z
i
R
-
1
m
∑
L
a
=
a
l
c
z
L
a
-
∑
L
b
=
l
c
+
1
c
z
L
b
=
n
m
n
+
m
1
/
2
N
′
c
D
k
,
j
-
M
′
c
D
k
-
1,2
j
-
1
,
D
k
,
j
;
R
=
D
m
n
′
k
,
j
;
k
-
1,2
j
=
n
m
n
+
m
1
/
2
1
n
∑
i
L
=
a
c
z
i
L
-
∑
i
R
=
c
+
1
b
z
i
R
-
1
m
∑
R
a
=
c
+
1
r
c
z
R
a
-
∑
R
b
=
r
c
+
1
b
z
R
b
=
n
m
n
+
m
1
/
2
N
′
c
D
k
,
j
-
M
′
c
D
k
-
1,2
j
,where 2
≤
k
≤
l
o
g
2
N, 1
≤
j
≤
N
/
2
k; a
=
2
k
(
j
-
1
)
+
1, b
=
2
k
j, c
=
2
k
(
j
-
1
)
+
2
(
k
-
1
); l
c
=
2
k
(
j
-
1
)
+
2
(
k
-
2
), r
c
=
c
+
2
(
k
-
2
); n
=
2
k, m
=
2
k
-
1; and N
′
=
(
2
)
∧
k
/
2
k, M
′
=
(
2
)
∧
(
k
-
1
)
/
2
(
k
-
1
).
Suppose Criterion1 is invalid as S
k
,
j
;
L
=
S
k
,
j
;
R
∥
(
m
a
x
(
S
k
,
j
;
L
,
S
k
,
j
;
R
)
≤
C
(
α
)
) holds true; the second search criterion needs to be introduced in terms of the two variance fluctuation variables D
k
,
j
;
L and D
k
,
j
;
R as follows.Criterion 2.
Given two variance fluctuation variablesD
k
,
j
;
L and D
k
,
j
;
R according to the two non-leaf child nodes c
D
k
-
1,2
j
-
1 and c
D
k
-
1,2
j of the selected node c
D
k
,
j in BSTcD, and 2
≤
k
≤
l
o
g
2
N, (a)
if(
D
k
,
j
;
L
>
D
k
,
j
;
R
)
∧
(
D
k
,
j
;
L
>
C
(
β
)
) holds true, then the left-child node c
A
k
-
1,2
j
-
1 in BSTcA is accordingly selected and involved into the current search path; meanwhile the right one is ignored;
(b)
if(
D
k
,
j
;
L
<
D
k
,
j
;
R
)
∧
(
D
k
,
j
;
R
>
C
(
β
)
) holds true, then the right-child node c
A
k
-
1,2
j in BSTcA is accordingly selected and involved into the current search path; meanwhile the left one is ignored.Proof.
Similarly, as illustrated in Figure4, the satisfied (
D
k
,
j
;
L
>
D
k
,
j
;
R
) in Criterion 2 means that the variance fluctuations within Z
L are stronger than that one within Z
R, in terms of the definitions of D
k
,
j
;
L and D
k
,
j
;
R in (11). That is, a potential abrupt change might exist in Z
L with more probability than in Z
R, and vice versa. Meanwhile, if D
k
,
j
;
L
>
C
(
β
) holds true, then (
H
1
) in Hypothesis 2 is satisfied; that is, abrupt change occurs in Z
L, and vice versa, where C
(
β
) is the critical value predefined in an identical distribution and β is the significance level. As a result, one of the two non-leaf child nodes c
A
k
-
1,2
j
-
1 and c
A
k
-
1,2
j in BSTcA can be accordingly selected, and the other one is neglected. Therefore, if Criterion 1 is invalid for less significant statistic fluctuation within BSTcA, Criterion 2 ensures that the search procedure can keep going forward to the potential abrupt change in Z, especially when abrupt change occurs near the left or the right boundary of Z with smaller size N.Figure 4
The scheme of Criterion2 based on the variance fluctuations within BSTcD. Supposing Criterion 1 is invalid for insignificant statistic fluctuation within BSTcA, Criterion 2 ensures that one of the two non-leaf child nodes c
A
k
-
1,2
j
-
1 and c
A
k
-
1,2
j can also be selected from BSTcA, in accordance with the variance fluctuation within BSTcD. Therefore, the search procedure can keep going forward, to the potential abrupt change in Z.Based on Criterions1 and 2 above, a search path can be obtained from the top root to the last non-leaf levels of BSTcA. In order to estimate an abrupt change from the original elements of Z, another criterion needs to be introduced to discern which one can be selected from two adjacent leaf nodes in BSTcA.Definition 5.
Supposing the current nodec
A
k
,
j is selected in the last non-leaf level of BSTcA, k
=
1, with two child leaf nodes z
2
j
-
1 and z
2
j, two statistic fluctuation variables D
L and D
R are defined based on KS test as(12)
D
L
=
D
m
n
z
L
=
m
n
m
+
n
1
/
2
F
m
z
L
-
G
n
z
L
=
m
n
m
+
n
1
/
2
1
m
∑
i
=
1
2
j
-
1
I
z
i
≤
z
L
-
1
n
∑
h
=
2
j
N
I
z
h
≤
z
L
,
D
R
=
D
m
n
z
R
=
m
n
m
+
n
1
/
2
F
m
z
R
-
G
n
z
R
=
m
n
m
+
n
1
/
2
1
m
∑
i
=
1
2
j
I
z
i
≤
z
R
-
1
n
∑
h
=
2
j
+
1
N
I
z
h
≤
z
R
,where z
L
=
z
2
j
-
1 and z
R
=
z
2
j; F
m
(
z
) and G
n
(
z
) refer to the e.c.d.f of Z
L
=
{
z
1
,
…
,
z
m
} and Z
R
=
{
z
m
+
1
,
…
,
z
N
}, respectively; m
=
2
j
-
1 or 2
j and n
=
N
-
m.
Consider that the largest statistic fluctuation betweenF
m
(
x
) and G
n
(
x
) is achieved either before or after one of the jumps, that is,(13)
s
u
p
x
∈
R
G
n
x
-
F
m
x
=
m
a
x
1
≤
i
≤
n
F
m
z
i
-
-
G
n
z
i
-
b
e
f
o
r
e
t
h
e
i
t
h
j
u
m
p
F
m
z
i
-
G
n
z
i
a
f
t
e
r
t
h
e
i
t
h
j
u
m
p
.Then, another two variables D
L
- and D
R
- are defined as(14)
D
L
-
=
D
m
n
z
L
-
=
m
n
m
+
n
1
/
2
F
m
z
L
-
-
G
n
z
L
-
=
m
n
m
+
n
1
/
2
1
m
∑
i
=
1
2
j
-
1
I
z
i
<
z
L
-
1
n
∑
h
=
2
j
N
I
z
h
<
z
L
,
D
R
-
=
D
m
n
z
R
-
=
m
n
m
+
n
1
/
2
F
m
z
R
-
-
G
n
z
R
-
=
m
n
m
+
n
1
/
2
1
m
∑
i
=
1
2
j
-
1
I
z
i
<
z
R
-
1
n
∑
h
=
2
j
N
I
z
h
<
z
R
.Therefore, the maximal statistic fluctuations D
L
′ and D
R
′ can be selected from D
L
- and D
L, as well as D
R
- and D
R. Then, the third search criterion is introduced in terms of D
L
′ and D
R
′ as follows.Criterion 3.
GivenD
L
′ and D
R
′ in accordance with two child leaf nodes z
2
j
-
1 and z
2
j of the selected non-leaf node c
A
k
,
j in BSTcA, k
=
1, (a)
if(
m
a
x
(
D
L
′
,
D
R
′
)
=
D
L
′
)
∧
(
D
L
′
>
C
(
γ
)
) holds true, then the left leaf node z
2
j
-
1 in Z is taken as the estimated CP, and the right one is neglected;
(b)
if(
m
a
x
(
D
L
′
,
D
R
′
)
=
D
R
′
)
∧
(
D
R
′
>
C
(
γ
)
) holds true, then the right leaf node z
2
j in Z is taken as the estimated CP, and the left one is neglected;
(c)
otherwise, no abrupt change is detected fromZ.Proof.
Obviously, ifm
a
x
(
D
L
′
,
D
R
′
)
>
C
(
γ
) is satisfied in Criterion 3, then the statistic fluctuation overtakes the critical value C
(
γ
) which is given in an identical data distribution, and γ is the significance level. Therefore, one of the two leaf nodes z
2
j
-
1 and z
2
j is taken as the estimated CP within Z.Supposing a non-leaf nodec
A
k
,
j is selected in BSTcA, the statistic and variance fluctuations are accordingly calculated between two adjacent segments Z
L and Z
R. Meanwhile, the search procedure is implemented from the root to non-leaf nodes in the last second level of BSTcA, in terms of Criterions 1 and 2. Then, the estimated CP can be obtained from the leaf nodes in BSTcA, by using Criterion 3. Thereafter, an optimal path towards a potential CP within Z is detected from BSTcA, after about l
o
g
2
N binary search steps.


### 3.3. Methods Compared with BSTKS

There are many methods proposed for abrupt change detection in time series, and the following are some typical methods, to evaluate the proposed BSTKS framework.KS Statistic (see [31]). In this method, a diagnosed time series Z is divided into two adjacent segments Z
L
=
{
z
1
,
z
2
,
…
,
z
m
} and Z
R
=
{
z
m
+
1
,
z
m
+
2
,
…
,
z
N
}, and then KS statistic is applied to calculate the statistic distance between Z
L and Z
R as(15)
D
m
n
x
=
m
n
N
1
/
2
s
u
p
x
∈
R
F
n
x
-
G
m
x
=
m
n
N
1
/
2
s
u
p
x
∈
R
∑
R
=
m
+
1
n
I
z
R
<
x
-
∑
L
=
1
m
I
z
L
<
x
,where F
n
(
x
) and G
m
(
x
) stand for the e.c.d.f of Z
L, and Z
R, respectively; N
=
m
+
n, N is the total length of Z, and m refers to a current test position within Z.t-Statistic (see [32]). t also known as Welch’s t-test is used only when the two population variances are assumed different (the two sample sizes may or may not be equal) and hence must be estimated separately. Suppose a diagnosed Z is divided into Z
L
=
{
z
1
,
z
2
,
…
,
z
m
} and Z
R
=
{
z
m
+
1
,
z
m
+
2
,
…
,
z
N
}. Then, t-statistic is calculated as(16)
t
=
Z
-
L
-
Z
-
R
S
Z
-
L
-
Z
-
R
,
S
Z
-
L
-
Z
-
R
=
S
1
2
m
+
S
2
2
n
,where Z
-
L and Z
-
R are the sample means of Z
L and Z
R, respectively; S is an unbiased estimator of the standard deviation, N
=
m
+
n, and m and n are the sizes of two segments Z
L and Z
R, respectively.SSA (see [12, 13, 33]). In SSA method, a windowed portion is chosen within a time series Z
=
{
z
1
,
z
2
,
…
,
z
N
}, where N is large enough and a window width m and the lag parameter M are set such that M
=
m
/
2, K
=
m
-
M
+
1. For each n
=
0,1
,
…
,
N
-
m
-
M, this method takes an interval of the time series [
n
+
1
,
n
+
m
] and then defines the M
×
K trajectory matrix X
n and describes the structure of the windowed portion as an L-dimensional subspace. If the structure changes further, it will not be well described by the computed subspace. Then, the distance between this subspace and the new trajectory vectors will increase; therefore, this increase will signal that an abrupt change occurs in Z.


## 3.1. Two BSTs’ Construction

In the first part of the proposed BSTKS method, two BSTs, that is, BSTcA and BSTcD, are constructed from a time series sampleZ, by using multilevel HWT. Generally, as shown in Figure 1, a discrete time series signal Z
=
{
z
1
,
z
2
,
…
,
z
N
} can be decomposed into the kth-level trend c
A
k and k-level fluctuations, that is, c
D
1
,
c
D
2
,
…
,
c
D
k, k
=
1,2
,
…
,
l
o
g
2
N. The k-level HWT is the mapping H
k defined as [13](5)
Z
→
H
k
c
A
k
∣
c
D
k
∣
c
D
k
-
1
∣
⋯
∣
c
D
2
∣
c
D
1
,and then, the mapping H
k can be represented by the approximation and detail coefficient matrices, termed McA and McD as follows: (6)
M
c
A
=
c
A
1,1
⋯
c
A
1
,
N
⋮
c
A
k
,
j
⋮
c
A
M
,
1
0
0
,
M
c
D
=
c
D
1,1
⋯
c
D
1
,
N
⋮
c
D
k
,
j
⋮
c
D
M
,
1
0
0
,where 0
≤
k
≤
M
=
l
o
g
2
N and 1
≤
j
≤
N
/
2
k.Figure 1
The diagram of a discrete time seriesZ decomposition by k-level HWT, which is composed of k-level cA and cD, that is, the average and difference coefficient vectors.Supposing the size of a diagnosedZ is divisible k times by 2, the jth element c
A
k
,
j in c
A
k and the jth element c
D
k
,
j in c
D
k can be denoted as(7)
c
A
k
,
j
=
1
2
∧
k
∑
i
=
a
b
z
i
,
c
D
k
,
j
=
1
2
∧
k
∑
L
=
a
c
z
L
-
∑
R
=
c
+
1
b
z
R
,where 1
≤
k
≤
l
o
g
2
N and 2
k
(
j
-
1
)
+
1
≤
i
≤
j
∗
2
k; a
=
2
k
(
j
-
1
)
+
1, c
=
2
k
(
j
-
1
)
+
2
(
k
-
1
), and b
=
2
k
∗
j.During two BSTs’ construction, as shown in Figure2, the non-leaf nodes in BSTcA and BSTcD are assembled by the k-level coefficient vectors of McA and McD, respectively; and then the leaf nodes are derived directly from the original time series Z. Therefore, the features of abrupt change in Z can be reflected and distributed into the different non-leaf nodes of BSTcA and BSTcD, in accordance with the k level coefficient vectors in McA and McD.Figure 2
The diagrams of two binary trees, BSTcA and BSTcD, which are constructed by McA and McD, as well as the original time seriesZ.
(a)
(b)


## 3.2. CP Detection Based on Three Search Criteria

To find an optimal path towards the potential CP within a given time seriesZ quickly and efficiently, some search criteria need to be introduced, and then the data exceptions can be detected from the root to leaf nodes of two BSTs. As for the statistic fluctuation within BSTcA, first, a new variable z
k
,
j is defined according to a current non-leaf node c
A
k
,
j in BSTcA,(8)
z
k
,
j
=
1
2
k
∑
i
=
a
b
z
i
=
2
∧
k
2
k
c
A
k
,
j
,where 1
≤
k
≤
l
o
g
2
N, 1
≤
j
≤
N
/
2
k; a
=
2
k
(
j
-
1
)
+
1, b
=
2
k
∗
j, and a
≤
i
≤
b. Then, the statistic fluctuation between two adjacent segments Z
L
=
{
z
a
,
…
,
z
c
} and Z
R
=
{
z
c
+
1
,
…
,
z
b
} can be defined by a modified KS statistic as(9)
S
m
n
′
k
,
j
=
n
m
n
+
m
1
/
2
1
n
∑
i
L
=
a
c
I
z
i
L
≤
z
k
,
j
-
1
m
∑
i
R
=
c
+
1
b
I
z
i
R
≤
z
k
,
j
,in which z
k
,
j is a new element defined in (8); m and n stand for the sizes of Z
L and Z
R, respectively; 1
≤
k
≤
l
o
g
2
N, 1
≤
j
≤
N
/
2
k
;
a
=
2
k
(
i
-
1
)
+
1, b
=
2
k
j, and c
=
2
k
(
j
-
1
)
+
2
(
k
-
1
). S
m
n
′
(
k
,
j
) measures the e.c.d.f difference between Z
L and Z
R, and the larger S
m
n
′
(
k
,
j
) means the more significant statistic fluctuation between Z
L and Z
R. Therefore, a potential abrupt change might occur at c in Z with more probability.Definition 3.
For a current non-leaf nodec
A
k
,
j in BSTcA, with its left and right-child nodes c
A
k
-
1,2
j
-
1 and c
A
k
-
1,2
j, the distance of e.c.d.f, S
k
,
j
;
L, and S
k
,
j
;
R can be defined as(10)
S
k
,
j
;
L
=
S
m
n
′
k
,
j
;
k
-
1,2
j
-
1
=
n
m
n
+
m
1
/
2
1
n
∑
i
=
a
b
I
z
i
≤
z
k
,
j
-
1
m
∑
i
L
=
a
c
I
z
i
L
≤
z
k
,
j
=
n
m
n
+
m
1
/
2
W
1
n
∑
i
=
a
b
I
z
i
≤
c
A
k
,
j
-
1
m
∑
i
L
=
a
c
I
z
i
L
≤
c
A
k
,
j
,
S
k
,
j
;
R
=
S
m
n
′
k
,
j
;
k
-
1,2
j
=
n
m
n
+
m
1
/
2
1
n
∑
i
=
a
b
I
z
i
≤
z
k
,
j
-
1
m
∑
i
R
=
c
+
1
b
I
z
i
R
≤
z
k
,
j
=
n
m
n
+
m
1
/
2
W
1
n
∑
i
=
a
b
I
z
i
≤
c
A
k
,
j
-
1
m
∑
i
R
=
c
+
1
b
I
z
i
R
≤
c
A
k
,
j
,where 2
≤
k
≤
l
o
g
2
N, 1
≤
j
≤
N
/
2
k; a
=
2
k
(
j
-
1
)
+
1, b
=
2
k
j, c
=
2
k
(
j
-
1
)
+
2
(
k
-
1
); n
=
2
k, m
=
2
k
-
1; and W
=
(
2
)
∧
k
/
2
k. To estimate an optimal path towards the potential change position within Z, without loss of generality, the first search criterion is introduced based on the statistic fluctuations S
k
,
j
;
L and S
k
,
j
;
R.Criterion 1.
Given two statistic fluctuation variablesS
k
,
j
;
L and S
k
,
j
;
R in accordance with two non-leaf child nodes c
A
k
-
1,2
j
-
1 and c
A
k
-
1,2
j of the current selected node c
A
k
,
j in BSTcA, and 2
≤
k
≤
l
o
g
2
N, (a)
if(
S
k
,
j
;
L
>
S
k
,
j
;
R
)
∧
(
S
k
,
j
;
L
>
C
(
α
)
) holds true, then the left-child node c
A
k
-
1,2
j
-
1 is selected and involved into the current search path; meanwhile, the right-child c
A
k
-
1,2
j is discarded;
(b)
if(
S
k
,
j
;
R
>
S
k
,
j
;
L
)
∧
(
S
k
,
j
;
R
>
C
(
α
)
) holds true, then the right-child node c
A
k
-
1,2
j is selected and involved into the current search path; meanwhile, the left-child c
A
k
-
1,2
j
-
1 is discarded.Proof.
For a selected non-leaf nodec
A
k
,
j in BSTcA, as shown in Figure 3, the original time series Z is divided equally into two adjacent segments Z
L and Z
R, which are covered by two non-leaf child nodes c
A
k
-
1,2
j
-
1 and c
A
k
-
1,2
j, respectively. According to the definitions of S
k
,
j
;
L and S
k
,
j
;
R in (10), the satisfied S
k
,
j
;
L
>
S
k
,
j
;
R indicates that the statistic fluctuation within Z
L is more significant than that one within Z
R; that is, a potential abrupt change might be contained in Z
L with more probability than in Z
R, and vice versa. Furthermore, if S
k
,
j
;
L
>
C
(
α
) holds true, then (
H
1
) of Hypothesis 1 is satisfied; that is, abrupt change occurs in Z
L, and vice versa, where C
(
α
) is the critical value predefined in an identical distribution and α is the significance level. Therefore, one of the two child nodes c
A
k
-
1,2
j
-
1 and c
A
k
-
1,2
j is selected and involved into the current search path; meanwhile, the remaining one is discarded. Once the statistic fluctuation is significant enough, an optimal search path can be detected by Criterion 1 from the top to the last non-leaf level in BSTcA. However, the search procedure is probably forced to cease because the statistic fluctuation is so insignificant that Criterion 1 is invalid for detecting it, especially for the left or the right boundary when sample Z is with smaller size N. Therefore, it is necessary to introduce another search criterion based on the variance fluctuations within BSTcD.Figure 3
The scheme of Criterion1 based on the statistic fluctuations within BSTcA. In terms of this criterion, the left or right-child node, that is, c
A
k
,
2
j
-
1 or c
A
k
,
2
j, might be selected to be involved in the current search path; meanwhile the remaining one is discarded. Thereafter, an optimal path towards the potential abrupt change in Z is expected to be obtained from BSTcA, after log
2
⁡
N binary search steps.Definition 4.
For a current non-leaf nodec
D
k
,
j in BSTcD, with its left and right-child nodes c
D
k
-
1,2
j
-
1 and c
D
k
-
1,2
j, respectively, the variance fluctuations D
k
,
j
;
L and D
k
,
j
;
R are defined in terms of (4) as(11)
D
k
,
j
;
L
=
D
m
n
′
k
,
j
;
k
-
1,2
j
-
1
=
n
m
n
+
m
1
/
2
1
n
∑
i
L
=
a
c
z
i
L
-
∑
i
R
=
c
+
1
b
z
i
R
-
1
m
∑
L
a
=
a
l
c
z
L
a
-
∑
L
b
=
l
c
+
1
c
z
L
b
=
n
m
n
+
m
1
/
2
N
′
c
D
k
,
j
-
M
′
c
D
k
-
1,2
j
-
1
,
D
k
,
j
;
R
=
D
m
n
′
k
,
j
;
k
-
1,2
j
=
n
m
n
+
m
1
/
2
1
n
∑
i
L
=
a
c
z
i
L
-
∑
i
R
=
c
+
1
b
z
i
R
-
1
m
∑
R
a
=
c
+
1
r
c
z
R
a
-
∑
R
b
=
r
c
+
1
b
z
R
b
=
n
m
n
+
m
1
/
2
N
′
c
D
k
,
j
-
M
′
c
D
k
-
1,2
j
,where 2
≤
k
≤
l
o
g
2
N, 1
≤
j
≤
N
/
2
k; a
=
2
k
(
j
-
1
)
+
1, b
=
2
k
j, c
=
2
k
(
j
-
1
)
+
2
(
k
-
1
); l
c
=
2
k
(
j
-
1
)
+
2
(
k
-
2
), r
c
=
c
+
2
(
k
-
2
); n
=
2
k, m
=
2
k
-
1; and N
′
=
(
2
)
∧
k
/
2
k, M
′
=
(
2
)
∧
(
k
-
1
)
/
2
(
k
-
1
).
Suppose Criterion1 is invalid as S
k
,
j
;
L
=
S
k
,
j
;
R
∥
(
m
a
x
(
S
k
,
j
;
L
,
S
k
,
j
;
R
)
≤
C
(
α
)
) holds true; the second search criterion needs to be introduced in terms of the two variance fluctuation variables D
k
,
j
;
L and D
k
,
j
;
R as follows.Criterion 2.
Given two variance fluctuation variablesD
k
,
j
;
L and D
k
,
j
;
R according to the two non-leaf child nodes c
D
k
-
1,2
j
-
1 and c
D
k
-
1,2
j of the selected node c
D
k
,
j in BSTcD, and 2
≤
k
≤
l
o
g
2
N, (a)
if(
D
k
,
j
;
L
>
D
k
,
j
;
R
)
∧
(
D
k
,
j
;
L
>
C
(
β
)
) holds true, then the left-child node c
A
k
-
1,2
j
-
1 in BSTcA is accordingly selected and involved into the current search path; meanwhile the right one is ignored;
(b)
if(
D
k
,
j
;
L
<
D
k
,
j
;
R
)
∧
(
D
k
,
j
;
R
>
C
(
β
)
) holds true, then the right-child node c
A
k
-
1,2
j in BSTcA is accordingly selected and involved into the current search path; meanwhile the left one is ignored.Proof.
Similarly, as illustrated in Figure4, the satisfied (
D
k
,
j
;
L
>
D
k
,
j
;
R
) in Criterion 2 means that the variance fluctuations within Z
L are stronger than that one within Z
R, in terms of the definitions of D
k
,
j
;
L and D
k
,
j
;
R in (11). That is, a potential abrupt change might exist in Z
L with more probability than in Z
R, and vice versa. Meanwhile, if D
k
,
j
;
L
>
C
(
β
) holds true, then (
H
1
) in Hypothesis 2 is satisfied; that is, abrupt change occurs in Z
L, and vice versa, where C
(
β
) is the critical value predefined in an identical distribution and β is the significance level. As a result, one of the two non-leaf child nodes c
A
k
-
1,2
j
-
1 and c
A
k
-
1,2
j in BSTcA can be accordingly selected, and the other one is neglected. Therefore, if Criterion 1 is invalid for less significant statistic fluctuation within BSTcA, Criterion 2 ensures that the search procedure can keep going forward to the potential abrupt change in Z, especially when abrupt change occurs near the left or the right boundary of Z with smaller size N.Figure 4
The scheme of Criterion2 based on the variance fluctuations within BSTcD. Supposing Criterion 1 is invalid for insignificant statistic fluctuation within BSTcA, Criterion 2 ensures that one of the two non-leaf child nodes c
A
k
-
1,2
j
-
1 and c
A
k
-
1,2
j can also be selected from BSTcA, in accordance with the variance fluctuation within BSTcD. Therefore, the search procedure can keep going forward, to the potential abrupt change in Z.Based on Criterions1 and 2 above, a search path can be obtained from the top root to the last non-leaf levels of BSTcA. In order to estimate an abrupt change from the original elements of Z, another criterion needs to be introduced to discern which one can be selected from two adjacent leaf nodes in BSTcA.Definition 5.
Supposing the current nodec
A
k
,
j is selected in the last non-leaf level of BSTcA, k
=
1, with two child leaf nodes z
2
j
-
1 and z
2
j, two statistic fluctuation variables D
L and D
R are defined based on KS test as(12)
D
L
=
D
m
n
z
L
=
m
n
m
+
n
1
/
2
F
m
z
L
-
G
n
z
L
=
m
n
m
+
n
1
/
2
1
m
∑
i
=
1
2
j
-
1
I
z
i
≤
z
L
-
1
n
∑
h
=
2
j
N
I
z
h
≤
z
L
,
D
R
=
D
m
n
z
R
=
m
n
m
+
n
1
/
2
F
m
z
R
-
G
n
z
R
=
m
n
m
+
n
1
/
2
1
m
∑
i
=
1
2
j
I
z
i
≤
z
R
-
1
n
∑
h
=
2
j
+
1
N
I
z
h
≤
z
R
,where z
L
=
z
2
j
-
1 and z
R
=
z
2
j; F
m
(
z
) and G
n
(
z
) refer to the e.c.d.f of Z
L
=
{
z
1
,
…
,
z
m
} and Z
R
=
{
z
m
+
1
,
…
,
z
N
}, respectively; m
=
2
j
-
1 or 2
j and n
=
N
-
m.
Consider that the largest statistic fluctuation betweenF
m
(
x
) and G
n
(
x
) is achieved either before or after one of the jumps, that is,(13)
s
u
p
x
∈
R
G
n
x
-
F
m
x
=
m
a
x
1
≤
i
≤
n
F
m
z
i
-
-
G
n
z
i
-
b
e
f
o
r
e
t
h
e
i
t
h
j
u
m
p
F
m
z
i
-
G
n
z
i
a
f
t
e
r
t
h
e
i
t
h
j
u
m
p
.Then, another two variables D
L
- and D
R
- are defined as(14)
D
L
-
=
D
m
n
z
L
-
=
m
n
m
+
n
1
/
2
F
m
z
L
-
-
G
n
z
L
-
=
m
n
m
+
n
1
/
2
1
m
∑
i
=
1
2
j
-
1
I
z
i
<
z
L
-
1
n
∑
h
=
2
j
N
I
z
h
<
z
L
,
D
R
-
=
D
m
n
z
R
-
=
m
n
m
+
n
1
/
2
F
m
z
R
-
-
G
n
z
R
-
=
m
n
m
+
n
1
/
2
1
m
∑
i
=
1
2
j
-
1
I
z
i
<
z
R
-
1
n
∑
h
=
2
j
N
I
z
h
<
z
R
.Therefore, the maximal statistic fluctuations D
L
′ and D
R
′ can be selected from D
L
- and D
L, as well as D
R
- and D
R. Then, the third search criterion is introduced in terms of D
L
′ and D
R
′ as follows.Criterion 3.
GivenD
L
′ and D
R
′ in accordance with two child leaf nodes z
2
j
-
1 and z
2
j of the selected non-leaf node c
A
k
,
j in BSTcA, k
=
1, (a)
if(
m
a
x
(
D
L
′
,
D
R
′
)
=
D
L
′
)
∧
(
D
L
′
>
C
(
γ
)
) holds true, then the left leaf node z
2
j
-
1 in Z is taken as the estimated CP, and the right one is neglected;
(b)
if(
m
a
x
(
D
L
′
,
D
R
′
)
=
D
R
′
)
∧
(
D
R
′
>
C
(
γ
)
) holds true, then the right leaf node z
2
j in Z is taken as the estimated CP, and the left one is neglected;
(c)
otherwise, no abrupt change is detected fromZ.Proof.
Obviously, ifm
a
x
(
D
L
′
,
D
R
′
)
>
C
(
γ
) is satisfied in Criterion 3, then the statistic fluctuation overtakes the critical value C
(
γ
) which is given in an identical data distribution, and γ is the significance level. Therefore, one of the two leaf nodes z
2
j
-
1 and z
2
j is taken as the estimated CP within Z.Supposing a non-leaf nodec
A
k
,
j is selected in BSTcA, the statistic and variance fluctuations are accordingly calculated between two adjacent segments Z
L and Z
R. Meanwhile, the search procedure is implemented from the root to non-leaf nodes in the last second level of BSTcA, in terms of Criterions 1 and 2. Then, the estimated CP can be obtained from the leaf nodes in BSTcA, by using Criterion 3. Thereafter, an optimal path towards a potential CP within Z is detected from BSTcA, after about l
o
g
2
N binary search steps.


## 3.3. Methods Compared with BSTKS

There are many methods proposed for abrupt change detection in time series, and the following are some typical methods, to evaluate the proposed BSTKS framework.KS Statistic (see [31]). In this method, a diagnosed time series Z is divided into two adjacent segments Z
L
=
{
z
1
,
z
2
,
…
,
z
m
} and Z
R
=
{
z
m
+
1
,
z
m
+
2
,
…
,
z
N
}, and then KS statistic is applied to calculate the statistic distance between Z
L and Z
R as(15)
D
m
n
x
=
m
n
N
1
/
2
s
u
p
x
∈
R
F
n
x
-
G
m
x
=
m
n
N
1
/
2
s
u
p
x
∈
R
∑
R
=
m
+
1
n
I
z
R
<
x
-
∑
L
=
1
m
I
z
L
<
x
,where F
n
(
x
) and G
m
(
x
) stand for the e.c.d.f of Z
L, and Z
R, respectively; N
=
m
+
n, N is the total length of Z, and m refers to a current test position within Z.t-Statistic (see [32]). t also known as Welch’s t-test is used only when the two population variances are assumed different (the two sample sizes may or may not be equal) and hence must be estimated separately. Suppose a diagnosed Z is divided into Z
L
=
{
z
1
,
z
2
,
…
,
z
m
} and Z
R
=
{
z
m
+
1
,
z
m
+
2
,
…
,
z
N
}. Then, t-statistic is calculated as(16)
t
=
Z
-
L
-
Z
-
R
S
Z
-
L
-
Z
-
R
,
S
Z
-
L
-
Z
-
R
=
S
1
2
m
+
S
2
2
n
,where Z
-
L and Z
-
R are the sample means of Z
L and Z
R, respectively; S is an unbiased estimator of the standard deviation, N
=
m
+
n, and m and n are the sizes of two segments Z
L and Z
R, respectively.SSA (see [12, 13, 33]). In SSA method, a windowed portion is chosen within a time series Z
=
{
z
1
,
z
2
,
…
,
z
N
}, where N is large enough and a window width m and the lag parameter M are set such that M
=
m
/
2, K
=
m
-
M
+
1. For each n
=
0,1
,
…
,
N
-
m
-
M, this method takes an interval of the time series [
n
+
1
,
n
+
m
] and then defines the M
×
K trajectory matrix X
n and describes the structure of the windowed portion as an L-dimensional subspace. If the structure changes further, it will not be well described by the computed subspace. Then, the distance between this subspace and the new trajectory vectors will increase; therefore, this increase will signal that an abrupt change occurs in Z.


## 4. Results and Discussion

In this section, the proposed BSTKS is evaluated on the synthetic time series and real EEG recordings with different sizeN. By comparing with existing KS, t, and SSA methods, the efficiency, sensitivity, and performance are analyzed in terms of the computation time, error and accuracy, hit rate, and AUC of ROC analyses. Furthermore, the novelty of our algorithm and necessity for real application are discussed in the following paragraphs.


### 4.1. CP Detection on Synthetic Time Series

In our simulations, some typical time series samples were derived from the normally distributed datasets (mean,u
=
0, and standard deviation, sd = 1). Each diagnosed sample of size N is composed of a normal segment of size k and an adjacent segment of size N
-
k, in which the abnormal part is simulated by adding a constant variation v into the random numbers of size N
-
k. The proposed BSTKS and other three methods, namely, KS, t, and SSA, were tested, respectively, on 200 samples which were derived from each time series group G
i with N
i
=
2
∧
(
4
+
i
), i
=
1,2
,
…
,
7, and v
i
=
d
(
1
+
l
o
g
2
(
k
-
4
)
), where k
=
l
o
g
2
(
N
i
) and d
=
1.0. For each sample in G
i, a series of test positions were arranged by C
P
K
j
=
j
∗
(
2
∧
(
k
-
4
)
), k
≥
5, and j
=
1,2
,
…
,
15.First, simulations were carried out according to different value of sample sizeN
i and test position C
P
K
j. The average analyses on four methods were listed in Table 1, and the results of simulations on datasets G
1
–
G
7 were illustrated in Figure 5. In general, our BSTKS is the most promising with the shortest computation time, the highest hit rate, the smallest error, and the highest accuracy out of all four methods. Particularly, as sample size N increases from N
1 to N
7, all four methods take longer time for bigger N, and BSTKS is always the fastest one. Meanwhile, BSTKS owns the highest level of hit rate against the low tracks of other three methods; and BSTKS is much more efficient with the smallest error and the highest accuracy, though all four methods tend to be better with N increasing. However, BSTKS has smaller AUC of ROC analyses, that is, bigger search space, than SSA and KS.Table 1
The averaged results on four methods with datasetsG
1 to G
7.
Time
Hit rate
Error
Accuracy
AUC
BSTKS
.0063
.4797
38.5268
.9018
.8922
KS
.3537
.0841
38.7321
.8804
.8984
t
1.0068
.0168
56.3036
.8878
.7960
SSA
1.5218
.0583
41.9464
.8762
.9941Figure 5
The simulations onG
1 to G
7 with size N from 25 to 211. (a) The results in terms of computation time, hit rate, error and accuracy, and AUC of ROC analyses. (b) The average analyses on BSTKS and other three methods. In the histograms, “1,” “2,” “3,” and “4” stand for BSTKS, KS, t, and SSA, respectively.
(a)
Performance analysis on differentN
(b)
Averaged analysisSecond, simulations were carried out based on the datasetsG
1, G
4, and G
7. The proposed BSTKS and other three methods were tested according to the different value of variance v
=
d
(
1
+
log
⁡
2
(
k
-
4
)
), k
=
5,8
,
11, and d
=
0.5,1.0,2.0,3.0, respectively. The average results of four methods on G
1, G
4, and G
7 were summarized in Table 2, and the typical simulations were selected on G
4 and represented in Figure 6. Generally, when v gets larger, all four methods get better hit rate, accuracy, and AUC of ROC analysis, except for longer computation time for bigger size N. Compared with other three methods, the proposed BSTKS is more encouraging because of the shortest computation time, especially when N gets bigger, as well as the highest hit rate and accuracy, especially when N gets smaller. Moreover, the simulations on G
4 with different variance v (Figure 6) explicitly illustrate that BSTKS has the best performance when v gets larger, in terms of the shortest time and the biggest increase of the hit rate out of four methods. For the accuracy and AUC, both BSTKS and KS keep higher sensitivity than t and SSA, as v increases from 0.5 to 3.0. Moreover, the simulations on G
1 and G
7 were omitted, because similar results can be obtained like G
4 above.Table 2
The summary of simulations according to different variances inG
1, G
4, and G
7.
Items
Methods
N
=
2
5
N
=
2
8
N
=
2
11
BSTKS
KS
t
SSA
BSTKS
KS
t
SSA
BSTKS
KS
t
SSA
d
=
0.5
Time
.018
.035
.227
.116
.029
.335
1.85
2.81
.060
6.96
17.2
24.6
Hit rate
.046
.005
.010
.038
.093
.093
.005
.025
.106
.056
.006
.034
Accuracy
.792
.515
.792
.748
.984
.995
.905
.899
.999
.999
.944
.884
AUC
.694
.694
.644
.997
.954
.951
.978
.971
.983
1.00
.999
1.00
d
=
1.0
Time
.018
.034
.223
.113
.031
.356
1.97
2.98
.061
7.09
18.1
26.5
Hit rate
.086
.013
.007
.035
.041
.001
.099
.142
.135
.045
.000
.035
Accuracy
.846
.552
.839
.756
.998
.998
.939
.974
.999
.999
.940
.986
AUC
.695
.695
.851
.997
.993
.992
.997
.991
.992
1.00
.998
1.00
d
=
2.0
Time
.018
.035
.229
.115
.031
.345
1.91
2.88
.065
7.60
19.7
29.1
Hit rate
.165
.061
.007
.049
.181
.090
.000
.049
.167
.028
.000
.053
Accuracy
.927
.737
.958
.765
.998
.997
.971
.983
.999
.999
.984
.997
AUC
.754
.754
.908
.997
.998
.999
.996
1.00
.995
1.00
.999
1.00
d
=
3.0
Time
.019
.037
.245
.125
.034
.382
2.08
3.28
.067
8.08
20.5
31.4
Hit rate
.225
.086
.002
.037
.189
.084
0.00
.046
.169
.035
0.00
.045
Accuracy
.942
.818
.952
.773
.997
.996
.986
.983
.999
.999
.991
.998
AUC
.857
.938
.655
.997
1.00
.996
.728
.100
1.00
.999
.467
1.00Figure 6
The simulations on 200 samples inG
4 with different variances. Under different variances v from 0.5 to 3.0, (a) the computation time, (b) the hit rate, (c) the error, (d) the accuracy, and (e) the AUC of ROC analysis, for BSTKS, KS, t, and SSA, respectively. In all “mean” histograms, “1,” “2,” “3,” and “4” in x-axis stand for BSTKS, KS, t, and SSA methods, respectively.
(a)
Time
(b)
Hit rate
(c)
Error
(d)
Accuracy
(e)
AUCThird, simulations were implemented based on different CP test positions withinG
1 and G
4. The proposed BSTKS and other three methods were analyzed according to the different value of test positionCPK and variance v. The results of simulations on G
1 and G
4 were illustrated in Figure 7, and the results near the left and right boundaries in G
1 and G
4 were summarized in Table 3. In general, all four methods tend to be better, when N increases under a fixed v, or when v increases under a fixed N. Meanwhile, for test positionCPK near the left and right boundaries, the proposed BSTKS produces better performance than other three methods, because of the highest hit rate, the smallest error in all four methods, and higher accuracy and AUC than t and SSA. Moreover, the simulations on G
1 and G
4 near the left and right boundaries were illustrated in Figure 8 in detail. In terms of the distribution of estimated CP (e-CP), PDF of e-CP, and AUC of ROC analysis, these simulations indicate that BSTKS is more sensitive for both left and right boundaries than other three methods, especially when sample size N and variance v get smaller.Table 3
The summary of simulations onG
1 and G
4 near the left and right boundaries according to different variance d.
Items
Methods
N
=
2
5, C
P
K
=
8
N
=
2
5, C
P
K
=
24
N
=
2
8, C
P
K
=
16
N
=
2
8, C
P
K
=
240
BSTKS
KS
t
SSA
BSTKS
KS
t
SSA
BSTKS
KS
t
SSA
BSTKS
KS
t
SSA
d
=
0.5
Hit rate
.200
.055
.036
.160
.215
.060
.010
0.0
.230
.160
0.0
0.0
.265
.165
0.0
.065
Error
2
20
8
3
3
7
8
12
24
1
134
49
28
21
102
33
Accuracy
.937
.375
.750
.906
.906
.781
.750
.625
.960
.996
.476
.808
.891
.918
.601
.871
AUC
.635
.963
.641
.978
.656
.987
.599
.978
.599
.946
.780
.797
.750
.884
.564
.797
d
=
1.0
Hit rate
.515
.160
.040
.295
.485
.190
.005
0.0
.490
.195
0.0
0.0
.535
.175
.025
.100
Error
0
10
4
2
0
3
8
12
1
0
112
10
0
1
90
1
Accuracy
1.0
.687
.875
.937
1.0
.906
.750
.625
.996
1.0
.562
.960
1.00
.996
.648
.996
AUC
.831
.883
.641
.978
.922
.927
.599
.978
.864
.986
.780
.829
.999
.988
.657
.911
d
=
2.0
Hit rate
.655
.240
0.0
.220
.725
.175
0.0
0.0
.510
.160
0.0
0.0
.530
.150
.005
.065
Error
0
6
1
2
0
1
4
12
0
0
107
5
0
1
36
4
Accuracy
1.00
.812
.968
.937
1.0
.968
.875
.625
1.0
1.0
.582
.980
1.0
.996
.859
.984
AUC
.976
.979
.770
.978
.938
.875
.808
.978
.978
.985
.780
.999
1.0
.990
.752
.995
d
=
3.0
Hit rate
.715
210
0.0
.265
.730
.215
0.0
0.0
.530
.195
0.0
0.0
.545
.145
0.0
.060
Error
0
6
1
2
0
1
1
13
0
0
119
5
0
2
11
4
Accuracy
1.0
.812
.968
.937
1.0
.968
.968
.593
1.0
1.0
.535
.980
1.0
.992
.957
.984
AUC
.999
.960
.770
.978
.996
.822
.808
.978
.999
.940
.780
.999
1.0
.990
.752
.998Figure 7
The simulations onG
1 and G
4 according to the different variance v and test positionCPK. The results were shown in (a) G
1 with N
=
32 and v
=
0.5, (b) G
4 with N
=
256 and v
=
0.5, (c) G
1 with N
=
32 and v
=
2.0, and (d) G
4 with N
=
256 and v
=
2.0, in terms of e-CP, hit rate, error, accuracy, and AUC, respectively.
(a)
Results of simulation onG
1 with N
=
32 and v
=
0.5
(b)
Results of simulation onG
4 with N
=
256 and v
=
0.5
(c)
Results of simulation withN
=
32 and v
=
2.0
(d)
Results of simulation onG
4 with N
=
256 and v
=
2.0Figure 8
The results of CP detection on the assembled EEG samplesZ
1
–
Z
4 with different value of sample size N and test positionCPK. For Z
1
–
Z
4 with different N from 27 to 210, (a1–a4) the assembled EEG samples Z
1
–
Z
4 with the assigned test positionCPK, (b1–b4) the e-CP, (c1–c4) the error of e-CP, (d1–d4) the accuracy of e-CP, and (e1–e4) the computation time. In the x-axis of (b–e), the methods “1,” “2,” “3,” and “4” stand for BSTKS, KS, t, and SSA, respectively.Therefore, all simulation results above suggest that our proposed BSTKS is an encouraging and efficient method for abrupt change detection from the synthetic time series datasets, because of the shortest computation time, the highest hit rate, and accuracy out of four methods, especially for less significant statistic fluctuation whenN gets smaller, as well as for less significant variance fluctuation when N gets bigger, and v gets smaller.


### 4.2. Abrupt Change Analyses on EEG Recordings

To verify the proposed method further, we take some representative samples from the CHBMIT Scalp EEG Database. In the PhysioBank platform, the CHBMIT Scalp EEG Database (CHBMIT) was collected at the Children’s Hospital Boston; it consists of EEG recordings from pediatric subjects with intractable seizures [34, 35]. In this CHBMIT EEG database, some subjects were monitored up to several days after withdrawal of antiseizure medication in order to characterize their seizures and assess their candidacy for surgical intervention. Based on these EEG recordings in the CHBMIT EEG database, as well as some existing experiments in [36–39], the proposed BSTKS and other three methods were tested according to different value of test positionCPK and sample size N.First, a diagnosed EEG sampleZ
=
[
Z
L
,
Z
R
] was assembled from two significantly different segments, in which Z
L
=
{
z
1
,
…
,
z
C
P
K
} and Z
R
=
{
z
C
P
K
+
1
,
…
,
z
N
} were derived from chb01_04_edfm and chb01_05_edfm, respectively. Then, BSTKS and other three methods were tested on the assembled EEG recordings Z
1
–
Z
8, respectively, according to the different value of assigned test positionCPK and sample size N. The results of abrupt change detection on these assembled EEG samples were illustrated in Figure 8 and summarized in Table 4. Generally speaking, all four methods can roughly estimate the assigned test position from each assembled EEG recording and then divide it into two adjacent segments Z
L and Z
R. It is worth stressing that the proposed BSTKS can discern the different EEG segments accurately with the smallest error and the highest accuracy out of four methods. Also, BSTKS is the most efficient and encouraging with the shortest time in all four methods.Table 4
The summary of abrupt change detection onZ
1
–
Z
8.
M
Z
Mean
N
27
28
29
210
27
28
29
210
C
P
K
25
50
100
200
100
200
400
900
e-CP
BSTKS
25
50
100
276
100
200
376
900
NA
KS
29
36
95
206
34
92
301
795
NA
t
24
255
31
1023
31
199
33
1023
NA
SSA
32
55
398
1007
106
208
500
907
NA
Err
BSTKS
0
0
0
76
0
0
24
0
12.5
KS
4
14
5
6
66
108
99
105
50.9
t
1
205
69
823
69
1
367
123
207.3
SSA
7
5
298
807
6
8
100
7
154.8
Acc
BSTKS
1.0
1.0
1.0
.93
1.0
1.0
.95
1.0
.98
KS
.97
.94
.99
.99
.48
.57
.81
.89
.83
t
.99
.20
.86
.20
.46
.99
.28
.88
.61
SSA
.94
.97
.42
.21
.94
.97
.80
.99
.78
Time
BSTKS
.023
.031
0.034
.036
.028
.033
.035
.038
.032
KS
.019
.021
.038
.049
.020
.029
.039
.052
.033
t
.03
.063
.088
.170
.031
.050
.081
.174
.086
SSA
.037
.071
.126
.239
.035
.065
.118
.245
.117Moreover, forCPK near the left and right boundaries in Z
1
–
Z
8, BSTKS has much better sensitivity than other KS, t, and SSA methods because of the smallest error and the highest accuracy, especially for less statistic fluctuation when N gets smaller, as well as less significant variance fluctuation when N gets bigger. Supposing the assembled EEG sample indicates that a sharp transition of one’s mental situation occurs before and after a sudden attack or acute stimulation, it is meaningful to estimate the location of the abrupt change and the maximal difference of data distribution exists between two adjacent EEG segments. These experiments above suggest that the proposed BSTKS can successfully detect the change position where a sudden change occurs under a potential mental shock, more quickly and efficiently than KS, t, and SSA methods.Second, the original EEG samplesZ
1
–
Z
6 were selected directly from different recordings in the chb01_05_edfm; then the proposed BSTKS and other three methods were tested according to different sample size N. Because the distance of e.c.d.f (V.e.c.d.f) can partly reflect the data fluctuation between two adjacent EEG segments, we use this V.e.c.d.f variable to distinguish different performance of BSTKS and other three methods. The results of abrupt change analyses on these original EEG recordings were shown in Figure 9 and summarized in Table 5. For all methods above, they can estimate an abrupt change from each of these original EEG samples Z
1
–
Z
6 and then divide it into two adjacent EEG segments. Compared with other three methods, the proposed BSTKS is encouraging for the shortest time out of four methods. Moreover, BSTKS has bigger V.e.c.d.f than t and SSA, which means that it can more reasonably distinguish two adjacent EEG segments with different state of mental health. Although KS has the biggest V.e.c.d.f in all four methods, it takes much more search time than BSTKS, especially when sample size N is getting larger. In addition, t needs the longest search time out of four methods, and it is invalid occasionally, for example, for Z
1, Z
5, and Z
6.Table 5
The summary of CP detection from the original EEG samplesZ
1
–
Z
6.
M
N
Mean
29
210
211
212
213
214
e-CP
BSTKS
328
316
1286
2633
4352
6224
NA
KS
348
317
1342
2252
4673
5947
NA
t
511
314
17
4095
10
16383
NA
SSA
426
854
90
2634
408
11271
NA
V.e.c.d.f
BSTKS
.0649
.2608
.2822
.0997
.1318
.0388
.1464
KS
.4603
.3829
.4407
.3050
.3325
.2234
.3574
t
0
.1257
.5384
0
0
0
.1106
SSA
.1368
.0850
.1260
.0745
.0212
.0012
.0741
Time
BSTKS
.020
.020
.024
.030
.019
.0320
.0241
KS
.016
.041
.112
.466
1.461
5.638
1.289
t
.072
.137
.281
.913
1.726
4.709
1.306
SSA
.107
.209
.415
1.103
1.769
3.548
1.192Figure 9
The analyses of abrupt change on the original EEG samples, by BSTKS, KS,t, and SSA, respectively. (a–f) The results of CP detection from the original EEG recordings Z
1
–
Z
6, with N from 29 to 214, respectively.For these original EEG recordings with intractable seizures, it is of great concern to predict when and where a significant change happens from these EEG signals. This abrupt change probably indicates that a patient encounters a vertical transition from a previous mental status, and it is very important and helpful for diagnosing the patients with intractable seizures. These experiments on original EEG samples above indicate that the proposed BSTKS can not only accurately detect the change position, but also estimate the maximal difference of data distribution existing between two adjacent EEG segments, more quickly and efficiently than existing KS,t, and SSA methods.


## 4.1. CP Detection on Synthetic Time Series

In our simulations, some typical time series samples were derived from the normally distributed datasets (mean,u
=
0, and standard deviation, sd = 1). Each diagnosed sample of size N is composed of a normal segment of size k and an adjacent segment of size N
-
k, in which the abnormal part is simulated by adding a constant variation v into the random numbers of size N
-
k. The proposed BSTKS and other three methods, namely, KS, t, and SSA, were tested, respectively, on 200 samples which were derived from each time series group G
i with N
i
=
2
∧
(
4
+
i
), i
=
1,2
,
…
,
7, and v
i
=
d
(
1
+
l
o
g
2
(
k
-
4
)
), where k
=
l
o
g
2
(
N
i
) and d
=
1.0. For each sample in G
i, a series of test positions were arranged by C
P
K
j
=
j
∗
(
2
∧
(
k
-
4
)
), k
≥
5, and j
=
1,2
,
…
,
15.First, simulations were carried out according to different value of sample sizeN
i and test position C
P
K
j. The average analyses on four methods were listed in Table 1, and the results of simulations on datasets G
1
–
G
7 were illustrated in Figure 5. In general, our BSTKS is the most promising with the shortest computation time, the highest hit rate, the smallest error, and the highest accuracy out of all four methods. Particularly, as sample size N increases from N
1 to N
7, all four methods take longer time for bigger N, and BSTKS is always the fastest one. Meanwhile, BSTKS owns the highest level of hit rate against the low tracks of other three methods; and BSTKS is much more efficient with the smallest error and the highest accuracy, though all four methods tend to be better with N increasing. However, BSTKS has smaller AUC of ROC analyses, that is, bigger search space, than SSA and KS.Table 1
The averaged results on four methods with datasetsG
1 to G
7.
Time
Hit rate
Error
Accuracy
AUC
BSTKS
.0063
.4797
38.5268
.9018
.8922
KS
.3537
.0841
38.7321
.8804
.8984
t
1.0068
.0168
56.3036
.8878
.7960
SSA
1.5218
.0583
41.9464
.8762
.9941Figure 5
The simulations onG
1 to G
7 with size N from 25 to 211. (a) The results in terms of computation time, hit rate, error and accuracy, and AUC of ROC analyses. (b) The average analyses on BSTKS and other three methods. In the histograms, “1,” “2,” “3,” and “4” stand for BSTKS, KS, t, and SSA, respectively.
(a)
Performance analysis on differentN
(b)
Averaged analysisSecond, simulations were carried out based on the datasetsG
1, G
4, and G
7. The proposed BSTKS and other three methods were tested according to the different value of variance v
=
d
(
1
+
log
⁡
2
(
k
-
4
)
), k
=
5,8
,
11, and d
=
0.5,1.0,2.0,3.0, respectively. The average results of four methods on G
1, G
4, and G
7 were summarized in Table 2, and the typical simulations were selected on G
4 and represented in Figure 6. Generally, when v gets larger, all four methods get better hit rate, accuracy, and AUC of ROC analysis, except for longer computation time for bigger size N. Compared with other three methods, the proposed BSTKS is more encouraging because of the shortest computation time, especially when N gets bigger, as well as the highest hit rate and accuracy, especially when N gets smaller. Moreover, the simulations on G
4 with different variance v (Figure 6) explicitly illustrate that BSTKS has the best performance when v gets larger, in terms of the shortest time and the biggest increase of the hit rate out of four methods. For the accuracy and AUC, both BSTKS and KS keep higher sensitivity than t and SSA, as v increases from 0.5 to 3.0. Moreover, the simulations on G
1 and G
7 were omitted, because similar results can be obtained like G
4 above.Table 2
The summary of simulations according to different variances inG
1, G
4, and G
7.
Items
Methods
N
=
2
5
N
=
2
8
N
=
2
11
BSTKS
KS
t
SSA
BSTKS
KS
t
SSA
BSTKS
KS
t
SSA
d
=
0.5
Time
.018
.035
.227
.116
.029
.335
1.85
2.81
.060
6.96
17.2
24.6
Hit rate
.046
.005
.010
.038
.093
.093
.005
.025
.106
.056
.006
.034
Accuracy
.792
.515
.792
.748
.984
.995
.905
.899
.999
.999
.944
.884
AUC
.694
.694
.644
.997
.954
.951
.978
.971
.983
1.00
.999
1.00
d
=
1.0
Time
.018
.034
.223
.113
.031
.356
1.97
2.98
.061
7.09
18.1
26.5
Hit rate
.086
.013
.007
.035
.041
.001
.099
.142
.135
.045
.000
.035
Accuracy
.846
.552
.839
.756
.998
.998
.939
.974
.999
.999
.940
.986
AUC
.695
.695
.851
.997
.993
.992
.997
.991
.992
1.00
.998
1.00
d
=
2.0
Time
.018
.035
.229
.115
.031
.345
1.91
2.88
.065
7.60
19.7
29.1
Hit rate
.165
.061
.007
.049
.181
.090
.000
.049
.167
.028
.000
.053
Accuracy
.927
.737
.958
.765
.998
.997
.971
.983
.999
.999
.984
.997
AUC
.754
.754
.908
.997
.998
.999
.996
1.00
.995
1.00
.999
1.00
d
=
3.0
Time
.019
.037
.245
.125
.034
.382
2.08
3.28
.067
8.08
20.5
31.4
Hit rate
.225
.086
.002
.037
.189
.084
0.00
.046
.169
.035
0.00
.045
Accuracy
.942
.818
.952
.773
.997
.996
.986
.983
.999
.999
.991
.998
AUC
.857
.938
.655
.997
1.00
.996
.728
.100
1.00
.999
.467
1.00Figure 6
The simulations on 200 samples inG
4 with different variances. Under different variances v from 0.5 to 3.0, (a) the computation time, (b) the hit rate, (c) the error, (d) the accuracy, and (e) the AUC of ROC analysis, for BSTKS, KS, t, and SSA, respectively. In all “mean” histograms, “1,” “2,” “3,” and “4” in x-axis stand for BSTKS, KS, t, and SSA methods, respectively.
(a)
Time
(b)
Hit rate
(c)
Error
(d)
Accuracy
(e)
AUCThird, simulations were implemented based on different CP test positions withinG
1 and G
4. The proposed BSTKS and other three methods were analyzed according to the different value of test positionCPK and variance v. The results of simulations on G
1 and G
4 were illustrated in Figure 7, and the results near the left and right boundaries in G
1 and G
4 were summarized in Table 3. In general, all four methods tend to be better, when N increases under a fixed v, or when v increases under a fixed N. Meanwhile, for test positionCPK near the left and right boundaries, the proposed BSTKS produces better performance than other three methods, because of the highest hit rate, the smallest error in all four methods, and higher accuracy and AUC than t and SSA. Moreover, the simulations on G
1 and G
4 near the left and right boundaries were illustrated in Figure 8 in detail. In terms of the distribution of estimated CP (e-CP), PDF of e-CP, and AUC of ROC analysis, these simulations indicate that BSTKS is more sensitive for both left and right boundaries than other three methods, especially when sample size N and variance v get smaller.Table 3
The summary of simulations onG
1 and G
4 near the left and right boundaries according to different variance d.
Items
Methods
N
=
2
5, C
P
K
=
8
N
=
2
5, C
P
K
=
24
N
=
2
8, C
P
K
=
16
N
=
2
8, C
P
K
=
240
BSTKS
KS
t
SSA
BSTKS
KS
t
SSA
BSTKS
KS
t
SSA
BSTKS
KS
t
SSA
d
=
0.5
Hit rate
.200
.055
.036
.160
.215
.060
.010
0.0
.230
.160
0.0
0.0
.265
.165
0.0
.065
Error
2
20
8
3
3
7
8
12
24
1
134
49
28
21
102
33
Accuracy
.937
.375
.750
.906
.906
.781
.750
.625
.960
.996
.476
.808
.891
.918
.601
.871
AUC
.635
.963
.641
.978
.656
.987
.599
.978
.599
.946
.780
.797
.750
.884
.564
.797
d
=
1.0
Hit rate
.515
.160
.040
.295
.485
.190
.005
0.0
.490
.195
0.0
0.0
.535
.175
.025
.100
Error
0
10
4
2
0
3
8
12
1
0
112
10
0
1
90
1
Accuracy
1.0
.687
.875
.937
1.0
.906
.750
.625
.996
1.0
.562
.960
1.00
.996
.648
.996
AUC
.831
.883
.641
.978
.922
.927
.599
.978
.864
.986
.780
.829
.999
.988
.657
.911
d
=
2.0
Hit rate
.655
.240
0.0
.220
.725
.175
0.0
0.0
.510
.160
0.0
0.0
.530
.150
.005
.065
Error
0
6
1
2
0
1
4
12
0
0
107
5
0
1
36
4
Accuracy
1.00
.812
.968
.937
1.0
.968
.875
.625
1.0
1.0
.582
.980
1.0
.996
.859
.984
AUC
.976
.979
.770
.978
.938
.875
.808
.978
.978
.985
.780
.999
1.0
.990
.752
.995
d
=
3.0
Hit rate
.715
210
0.0
.265
.730
.215
0.0
0.0
.530
.195
0.0
0.0
.545
.145
0.0
.060
Error
0
6
1
2
0
1
1
13
0
0
119
5
0
2
11
4
Accuracy
1.0
.812
.968
.937
1.0
.968
.968
.593
1.0
1.0
.535
.980
1.0
.992
.957
.984
AUC
.999
.960
.770
.978
.996
.822
.808
.978
.999
.940
.780
.999
1.0
.990
.752
.998Figure 7
The simulations onG
1 and G
4 according to the different variance v and test positionCPK. The results were shown in (a) G
1 with N
=
32 and v
=
0.5, (b) G
4 with N
=
256 and v
=
0.5, (c) G
1 with N
=
32 and v
=
2.0, and (d) G
4 with N
=
256 and v
=
2.0, in terms of e-CP, hit rate, error, accuracy, and AUC, respectively.
(a)
Results of simulation onG
1 with N
=
32 and v
=
0.5
(b)
Results of simulation onG
4 with N
=
256 and v
=
0.5
(c)
Results of simulation withN
=
32 and v
=
2.0
(d)
Results of simulation onG
4 with N
=
256 and v
=
2.0Figure 8
The results of CP detection on the assembled EEG samplesZ
1
–
Z
4 with different value of sample size N and test positionCPK. For Z
1
–
Z
4 with different N from 27 to 210, (a1–a4) the assembled EEG samples Z
1
–
Z
4 with the assigned test positionCPK, (b1–b4) the e-CP, (c1–c4) the error of e-CP, (d1–d4) the accuracy of e-CP, and (e1–e4) the computation time. In the x-axis of (b–e), the methods “1,” “2,” “3,” and “4” stand for BSTKS, KS, t, and SSA, respectively.Therefore, all simulation results above suggest that our proposed BSTKS is an encouraging and efficient method for abrupt change detection from the synthetic time series datasets, because of the shortest computation time, the highest hit rate, and accuracy out of four methods, especially for less significant statistic fluctuation whenN gets smaller, as well as for less significant variance fluctuation when N gets bigger, and v gets smaller.


## 4.2. Abrupt Change Analyses on EEG Recordings

To verify the proposed method further, we take some representative samples from the CHBMIT Scalp EEG Database. In the PhysioBank platform, the CHBMIT Scalp EEG Database (CHBMIT) was collected at the Children’s Hospital Boston; it consists of EEG recordings from pediatric subjects with intractable seizures [34, 35]. In this CHBMIT EEG database, some subjects were monitored up to several days after withdrawal of antiseizure medication in order to characterize their seizures and assess their candidacy for surgical intervention. Based on these EEG recordings in the CHBMIT EEG database, as well as some existing experiments in [36–39], the proposed BSTKS and other three methods were tested according to different value of test positionCPK and sample size N.First, a diagnosed EEG sampleZ
=
[
Z
L
,
Z
R
] was assembled from two significantly different segments, in which Z
L
=
{
z
1
,
…
,
z
C
P
K
} and Z
R
=
{
z
C
P
K
+
1
,
…
,
z
N
} were derived from chb01_04_edfm and chb01_05_edfm, respectively. Then, BSTKS and other three methods were tested on the assembled EEG recordings Z
1
–
Z
8, respectively, according to the different value of assigned test positionCPK and sample size N. The results of abrupt change detection on these assembled EEG samples were illustrated in Figure 8 and summarized in Table 4. Generally speaking, all four methods can roughly estimate the assigned test position from each assembled EEG recording and then divide it into two adjacent segments Z
L and Z
R. It is worth stressing that the proposed BSTKS can discern the different EEG segments accurately with the smallest error and the highest accuracy out of four methods. Also, BSTKS is the most efficient and encouraging with the shortest time in all four methods.Table 4
The summary of abrupt change detection onZ
1
–
Z
8.
M
Z
Mean
N
27
28
29
210
27
28
29
210
C
P
K
25
50
100
200
100
200
400
900
e-CP
BSTKS
25
50
100
276
100
200
376
900
NA
KS
29
36
95
206
34
92
301
795
NA
t
24
255
31
1023
31
199
33
1023
NA
SSA
32
55
398
1007
106
208
500
907
NA
Err
BSTKS
0
0
0
76
0
0
24
0
12.5
KS
4
14
5
6
66
108
99
105
50.9
t
1
205
69
823
69
1
367
123
207.3
SSA
7
5
298
807
6
8
100
7
154.8
Acc
BSTKS
1.0
1.0
1.0
.93
1.0
1.0
.95
1.0
.98
KS
.97
.94
.99
.99
.48
.57
.81
.89
.83
t
.99
.20
.86
.20
.46
.99
.28
.88
.61
SSA
.94
.97
.42
.21
.94
.97
.80
.99
.78
Time
BSTKS
.023
.031
0.034
.036
.028
.033
.035
.038
.032
KS
.019
.021
.038
.049
.020
.029
.039
.052
.033
t
.03
.063
.088
.170
.031
.050
.081
.174
.086
SSA
.037
.071
.126
.239
.035
.065
.118
.245
.117Moreover, forCPK near the left and right boundaries in Z
1
–
Z
8, BSTKS has much better sensitivity than other KS, t, and SSA methods because of the smallest error and the highest accuracy, especially for less statistic fluctuation when N gets smaller, as well as less significant variance fluctuation when N gets bigger. Supposing the assembled EEG sample indicates that a sharp transition of one’s mental situation occurs before and after a sudden attack or acute stimulation, it is meaningful to estimate the location of the abrupt change and the maximal difference of data distribution exists between two adjacent EEG segments. These experiments above suggest that the proposed BSTKS can successfully detect the change position where a sudden change occurs under a potential mental shock, more quickly and efficiently than KS, t, and SSA methods.Second, the original EEG samplesZ
1
–
Z
6 were selected directly from different recordings in the chb01_05_edfm; then the proposed BSTKS and other three methods were tested according to different sample size N. Because the distance of e.c.d.f (V.e.c.d.f) can partly reflect the data fluctuation between two adjacent EEG segments, we use this V.e.c.d.f variable to distinguish different performance of BSTKS and other three methods. The results of abrupt change analyses on these original EEG recordings were shown in Figure 9 and summarized in Table 5. For all methods above, they can estimate an abrupt change from each of these original EEG samples Z
1
–
Z
6 and then divide it into two adjacent EEG segments. Compared with other three methods, the proposed BSTKS is encouraging for the shortest time out of four methods. Moreover, BSTKS has bigger V.e.c.d.f than t and SSA, which means that it can more reasonably distinguish two adjacent EEG segments with different state of mental health. Although KS has the biggest V.e.c.d.f in all four methods, it takes much more search time than BSTKS, especially when sample size N is getting larger. In addition, t needs the longest search time out of four methods, and it is invalid occasionally, for example, for Z
1, Z
5, and Z
6.Table 5
The summary of CP detection from the original EEG samplesZ
1
–
Z
6.
M
N
Mean
29
210
211
212
213
214
e-CP
BSTKS
328
316
1286
2633
4352
6224
NA
KS
348
317
1342
2252
4673
5947
NA
t
511
314
17
4095
10
16383
NA
SSA
426
854
90
2634
408
11271
NA
V.e.c.d.f
BSTKS
.0649
.2608
.2822
.0997
.1318
.0388
.1464
KS
.4603
.3829
.4407
.3050
.3325
.2234
.3574
t
0
.1257
.5384
0
0
0
.1106
SSA
.1368
.0850
.1260
.0745
.0212
.0012
.0741
Time
BSTKS
.020
.020
.024
.030
.019
.0320
.0241
KS
.016
.041
.112
.466
1.461
5.638
1.289
t
.072
.137
.281
.913
1.726
4.709
1.306
SSA
.107
.209
.415
1.103
1.769
3.548
1.192Figure 9
The analyses of abrupt change on the original EEG samples, by BSTKS, KS,t, and SSA, respectively. (a–f) The results of CP detection from the original EEG recordings Z
1
–
Z
6, with N from 29 to 214, respectively.For these original EEG recordings with intractable seizures, it is of great concern to predict when and where a significant change happens from these EEG signals. This abrupt change probably indicates that a patient encounters a vertical transition from a previous mental status, and it is very important and helpful for diagnosing the patients with intractable seizures. These experiments on original EEG samples above indicate that the proposed BSTKS can not only accurately detect the change position, but also estimate the maximal difference of data distribution existing between two adjacent EEG segments, more quickly and efficiently than existing KS,t, and SSA methods.


## 5. Conclusion

In this paper, a novel BSTKS method is proposed based on binary search trees and a modified KS statistic. In this method, two BSTs were constructed from a diagnosed time series by multilevel HWT, and then an optimal search path is detected from the root to leaf nodes of two BSTs in terms of three search criteria. The novelty of the proposed method is addressed by comparing with other KS,t, and SSA methods, and simulations on the synthetic time series indicate that the proposed BSTKS is more efficient due to the shortest time, the highest hit rate, and the smallest error and highest accuracy out of four methods. Meanwhile, BSTKS has better sensitivity than KS near the left and right boundaries, because of shorter search time, higher hit rate, and bigger AUC, especially when sample size N gets smaller with less significant statistic fluctuation. In addition, the necessity of the proposed method in the real domain is analyzed on real EEG recordings, and the results indicate that the proposed method can successfully discern an abrupt change and then obviously distinguish two adjacent EEG segments from the real EEG recordings. Through inspecting the significant fluctuation between adjacent segments signals, it is encouraging further for useful information inspection on all kinds of physiological and psychological time series signals. In a word, our BSTKS is a novel, efficient, and promising method for abrupt change analysis, and it is very helpful for useful information inspection on all kinds of real time series with different scales.



---

*Source: 8343187-2016-06-16.xml*