Journal of Machine Learning Research 1 (2017) 1-48

Submitted 4/00; Published 10/00

Scan B-Statistic for Kernel Change-Point Detection
Shuang Li

sli370@gatech.edu

H. Milton Stewart School of Industrial and Systems Engineering
Georgia Institute of Technology
Atlanta, GA 30332, USA

arXiv:1507.01279v4 [cs.LG] 24 Sep 2017

Yao Xie (corresponding author)

yao.xie@isye.gatech.edu

H. Milton Stewart School of Industrial and Systems Engineering
Georgia Institute of Technology
Atlanta, GA 30332, USA

Hanjun Dai

hanjundai@gatech.edu

College of Computing
Georgia Institute of Technology
Atlanta, GA 30332, USA

Le Song

lsong@cc.gatech.edu

College of Computing
Georgia Institute of Technology
Atlanta, GA 30332, USA

Editor:

Abstract
Detecting the emergence of an abrupt change-point is a classic problem in statistics and
machine learning. Kernel-based nonparametric statistics have been used for this task which
enjoy fewer assumptions on the distributions than the parametric approach and can handle
high-dimensional data. In this paper we focus on the scenario when the amount of background data is large, and propose two related computationally efficient kernel-based statistics for change-point detection, which are inspired by the recently developed B-statistics.
A novel theoretical result of the paper is the characterization of the tail probability of these
statistics using the change-of-measure technique, which focuses on characterizing the tail
of the detection statistics rather than obtaining its asymptotic distribution under the null
distribution. Such approximations are crucial to control the false alarm rate, which corresponds to the significance level in offline change-point detection and the average-run-length
in online change-point detection. Our approximations are shown to be highly accurate.
Thus, they provide a convenient way to find detection thresholds for both offline and online cases without the need to resort to the more expensive simulations or bootstrapping.
We show that our methods perform well on both synthetic data and real data.
Keywords: Change-Point Detection, Kernel-Based Tests, Online Algorithm, Tail Probability, False-Alarm Control.

1. Introduction
Given a sequence of samples, x1 , x2 , . . . , xt , from a domain X , we are interested in detecting
a possible change-point τ , such that before the change samples xi are i.i.d. with a null
c 2017 Shuang Li, Yao Xie, Hanjun Dai, Le Song.

Li, Xie, Dai and Song

distribution P , and after the change samples xi are i.i.d. with a distribution Q. Here, the
time horizon t is either fixed, t = T0 , which we call the offline or fixed-sample change-point
detection, or t is not fixed, meaning that one can keep getting new samples, which we call
the online or sequential change-point detection. In the offline setting, our goal is to detect
the existence of a change. In the online setting, our goal is to detect the emergence of
a change as soon as possible after it occurs. Here, we restrict our attention to detecting
one change-point. One such instance is seismic event detection (Ross and Ben-Zion, 2014),
where we would like to either detect the presence of a weak event in retrospect to better
understand the geophysical structure or detect the event as quickly as possible for online
monitoring.
Ideally, the detection algorithm should be free of distributional assumptions to be robust
when applied to real data. To achieve this goal, various kernel based statistics have been
proposed in the machine learning literature (Harchaoui et al., 2008; Enikeeva and Harchaoui,
2014; Zou et al., 2014b; Kifer et al., 2004; Song et al., 2013; Desobry et al., 2005), which are
suitable for vector observations and typically work better in real data. Kernel approaches
are distribution free and more robust as they provide consistent results over larger classes of
data distributions; albeit they can be less powerful in settings where a clear distributional
assumption can be made. However, most kernel based statistics cost O(n2 ) to compute
over n samples. In online change-point detection, the number of samples grows with time
and hence we cannot directly use the naive approach. Recently, Zaremba et al. (2013)
develop the so-called B-test statistic to reduce the computational complexity. The B-test
statistic samples N pairs of blocks of size B from the two-sample data, compute the unbiased
estimates of the kernel-based statistic between each pair, and then taken an average. The
computational complexity of the B-statistics reduces to O(BN 2 ) or equivalently O(N n)
instead of O(n2 ).
In this paper, we present two related scan B-statistics for change-point detection. The
proposed statistics are based on kernel maximum mean discrepancy (MMD) (Gretton
et al., 2012; Harchaoui et al., 2013). Our methods are inspired by the B-statistic but
differ in various ways to tailor to the need of change-point detection. Typically, there is a
small number of post-change samples (for instance, seismic events are relatively rare, and
in online change-point detection, one would like to detect the change quickly). But there is
a large amount of reference data. So when constructing the detection statistic, we reuse the
post-change samples for the test block, but construct multiple and disjoint reference blocks.
This leads to non-negligible dependence between the MMD statistics being averaged over
and we cannot use the approach in existing works to analyze them. Moreover, the scanning
nature of the proposed statistic also introduces non-negligible dependence. We construct
the reference and test blocks in a structured way, so that analytical expressions can be
obtained when analyzing the false alarm rates.
Our main theoretical contribution includes accurate theoretical approximations to the
false-alarm rate of scan B-statistics. Controlling false alarms is a key challenge in changepoint detection. Specifically, this means to quantify the significance level for offline changepoint detection, and the average run length (ARL) for online change-point detection. We
cannot directly use the null property of the B-statistic from the existing work, because the
scan statistics take the maximum of multiple B-statistics computed over highly overlapping
data blocks. Due to such strong correlations, one cannot use the central limit theorem or
2

Scan B-Statistic for Kernel Change-Point Detection

even the martingale central limit theorem as used by existing work. Instead, we adopt a
recently developed change-of-measure technique (see, e.g., Yakir (2013)) for scan statistics,
which are capable of dealing with the more challenging situation here. Such a technique is
quite general, although it has not been used much in the machine learning community.
Our contribution also includes: (1) obtaining a closed-form variance estimator, which
allows easy calculation of the B-statistics; (2) further improving the accuracy of our approximations by taking into account the skewness of the kernel-based statistics. The accuracy of
our approximations is validated by numerical examples. Finally, we demonstrate the good
performance of our method using real speech and human activity data.
1.1 Related work
Classic parametric approaches for change-point detection can be found in Siegmund (1985);
Tartakovsky et al. (2014). There has been an array of nonparametric change-point detection
methods. Notable non-parametric schemes for change-point detection include Gordon and
Pollak (1994); Picard (1985), which are designed for scalar observations and not suitable for
vector observations. Brodsky and Darkhovsky (2013) provide a comprehensive introduction
to the methodologies and applications of nonparametric change-point detection. Bibinger
et al. (2017) construct a nonparametric minimax-optimal test to discriminate continuous
paths with volatility jumps and prove weak convergence of the test statistic to an extreme
value distribution. In the online setting, Kifer et al. (2004) present a meta-algorithm which
compares data in some “reference window” to the data in the current window, using empirical distance measures that are not kernel-based; Desobry et al. (2005) detect abrupt
changes by comparing two sets of descriptors extracted online from the signal at each time
instant: the immediate past set and the immediate future set, and then use a soft margin
single-class support vector machine to build a dissimilarity measure (which is asymptotically equivalent to the Fisher ratio in the Gaussian case) in the feature space between
those sets without estimating densities as an intermediate step; Song et al. (2013) present
a density-ratio estimation method to detect change-points, fitting the density-ratio using a
non-parametric Gaussian kernel model, whose parameters are updated online via stochastic
gradient descent approach. Another important branch of nonparametric change-point detection method is based on Kolmogorov-Smirnov test (Jr, 1951; Lilliefors, 1967), which has
been used in biostatistics (Wilcox, 2005; Wang et al., 2014; Wan et al., 2014). The generalization of Kolmogorov-Smirnov test from the univariate setting to the multi-dimensional
setting is given by the Peacock’s test (Fasano and Franceschini, 1987), which, however, is
less convenient to use that the kernel-based statistic test.
Seminal works by Csörgő and Horváth (1988) studied kernel-based U -statistic for the
change-point detection. They show that the statistic indexed by the assumed change-point
location parameter k, after proper standardization and rescaling of time and magnitude,
converges in distribution to a Gaussian process under the null, and converges to a deterministic path in probability under the alternative, when the number of samples goes to
infinity. These results are useful for bounding the detection statistics under the null with
high-probability (hence, controlling the false detection), and for studying the consistency of
tests. Csörgö and Horváth (1997) and Serfling (2009) contain comprehensive discussions on
asymptotic theory of non-parametric statistics including U -statistics. Our scan B-statistic
3

Li, Xie, Dai and Song

can also be viewed as a form of U -statistic using appropriate definition of the kernel. The
main differences of these classic work from our proposed B-statistic are: (1) our statistic
uses B-test block decomposition and averaging to make the test statistic more computationally efficient; (2) our statistic is more challenging to analyze due to the block structure
and correlation introduced by scanning statistics; (3) our analytical approach is different:
Csörgő and Horváth (1988) leverage invariance principle to establish convergence of the
entire sample path; we focus on characterizing the tail probability of the statistic under the
null and use the change-of-measure technique to achieve good approximation accuracy.
Other existing works that also focus on establishing asymptotic distribution of the detection statistic under the null for controlling the false alarm rate including the following:
Harchaoui et al. (2008) present a maximum kernel Fisher discriminant ratio statistic and
study its asymptotic null distribution; Dehling et al. (2015) investigate the two-sample test
U -statistic for dependent data. Our approach is different from above, in that we focus on
directly approximating the tail of the detection statistic under the null, rather than trying
to obtain its asymptotic distribution. Moreover, traditional analyses are usually done for
offline change-point detection, while our analytical framework based on change-of-measure
can be applied to both offline and online change-point detection.
Change-point detection problems are related to the classical statistical two-sample test.
However, they are usually more challenging than two-sample test because the change-point
location τ is unknown. Hence, when forming the detection statistic, one has to “taking the
maximum” of the detection statistics. The statistic being maxed over are usually highly
correlated since they are computed using overlapping data. Hence, when analyzing the false
detection rate, more sophisticated probabilistic tools have to be used.
Our techniques for approximating false alarm rates differ from large-deviation techniques
(see, e.g., Dembo and Zeitouni (2009)), which establish exponential rate by which the
probability converges to zero. In certain scenarios, the first-order approximation obtained
from large-deviation technique may not be sufficient for choosing threshold. Our method
provides more refined approximations to include polynomial terms and constants.
Finally, there are also works taking different approaches than hypothesis test for changepoint detection. For instance, Harchaoui and Cappé (2007) develop a kernel-based multiplechange-point detection approach, where the optimal location to segment the data is obtained
by dynamic programming; Arlot et al. (2012) treat multiple change-point detection as an
estimation problem, develop a kernelized linear model and provide a non-asymptotic oracle inequality for the estimation error. In the offline setting, Zou et al. (2014b) study
a problem when there are s anomalous sequences out of n sequences to be detected, and
the test statistic is constructed using MMD; Matteson and James (2014) propose a nonparametric approach based on U -statistics, which is capable of consistently estimating an
unknown number of multiple change-point locations; Zou et al. (2014a) propose a nonparametric maximum likelihood approach, with the number of change-points determined
from the Bayesian information criterion (BIC) and the locations of the change-points are
estimated via dynamic programming
Our notations are standard. Let Ik denote the identity matrix of size k-by-k; [1, N ] =
{1, 2, . . . , N }. Let E[A; B] = E[A1B ] denote the expectation conditioned on event B, where
1B represents the indicator function that takes value 1 when the event B happens and takes
value 0, otherwise. Let Var(·) and Cov(·) denote the variance and the covariance. Let 0 and
4

Scan B-Statistic for Kernel Change-Point Detection

e denote vectors of all zeros and all ones, respectively. Let [Σ]ij denote the ij-th element
of a matrix Σ. In Section 4, EB , VarB , and CovB denote the values computed under the
distribution after the change-of-measure. Similarly, in Section 5, Et , Vart , and Covt denote
the values obtained under the distribution after the change-of-measure.

2. Background
We briefly review kernel-based statistics and the maximum mean discrepancy. A reproducing kernel Hilbert space (RKHS) F on X with a kernel k(x, x0 ) is a Hilbert space of
functions f (·) : X 7→ R with inner product h·, ·iF . Its element k(x, ·) satisfies the reproducing property: hf (·), k(x, ·)iF = f (x), and consequently, hk(x, ·), k(x0 , ·)iF = k(x, x0 ),
meaning that we can view the evaluation of a function f at any point x ∈ X as an inner
product. Commonly used RKHS kernel functions include Gaussian radial basis function
(RBF) kernel k(x, x0 ) = exp(− kx − x0 k2 /2σ 2 ) where σ > 0 is the kernel bandwidth, and
polynomial kernel k(x, x0 ) = (hx, x0 i + a)d where a > 0 and d ∈ N (Schölkopf and Smola,
2002). RKHS kernels can also be defined for sequences, graph and other structured object (Schölkopf et al., 2004). In this paper, if not otherwise stated, we will assume that
Gaussian RBF kernel is used.
Assume there are two sets, each containing n samples, each taking value on a general domain X , where X = {x1 , x2 , . . . , xn } are i.i.d. with a distribution P , and Y =
{y1 , y2 , . . . , yn } are i.i.d. with a distribution Q. The maximum mean discrepancy (MMD)
is defined as (Gretton et al., 2012)
MMD[F, P, Q] := sup {EX∼P [f (X)] − EY ∼Q [f (Y )]} .
f ∈F

An unbiased estimator of MMD2 can be obtained using U -statistic (Gretton et al., 2012)
X
1
MMD2u [F, X, Y ] =
h(xi , xj , yi , yj ),
(1)
n(n − 1)
i,j∈[1,n],i6=j

where h(·) is the kernel for U -statistic, defined in terms of an RKHS kernel as
h(xi , xj , yi , yj ) = k(xi , xj ) + k(yi , yj ) − k(xi , yj ) − k(xj , yi ).

(2)

Intuitively, the empirical test statistic MMD2u is expected to be small (close to zero) if
P = Q, and large if P and Q are far apart. The complexity for evaluating MMD2u is O(n2 ),
since we have to form the so-called Gram matrix for the data, which is of size n-by-n. Under
the null hypothesis, P = Q, the U -statistic is degenerate and has the same distribution as
an infinite sum of Chi-square variables.
To improve computational efficiency, recently, an alternative approach to eatimate
MMD2 has been proposed (called the B-test by (Zaremba et al., 2013)), which we refer
to as the B-statistic in the following. The key idea is to partition the n samples from P
and Q into N non-overlapping blocks, X1 , . . . , XN and Y1 , . . . , YN , each of size B. Then
compute MMD2u [F, Xi , Yi ] for each pair of blocks and then take the average:
MMD2B [F, X, Y ] =

N
1 X
MMD2u [F, Xi , Yi ].
N
i=1

5

Li, Xie, Dai and Song

Since B is constant and N is on the order of O(n), the computational complexity of
MMD2B [F, X, Y ] is O(B 2 n), which is significantly lower than the O(n2 ) complexity of
MMD2u [F, X, Y ]. Furthermore, by averaging MMD2u [F, Xi , Yi ] over blocks, when the blocks
are independent, the B-statistic is asymptotically normal under the null due to central limit
theorem. This property allows a simple threshold to be derived for the B-test.

3. Scan B-statistic in change-point detection
Now we present our scan B-statistic based change-point detection schemes. Given a sequence of data {. . . , x−2 , x−1 , x0 , x1 , . . . , xt }, each taking value on a general domain X .
Let {. . . , x−2 , x−1 , x0 } denote the reference data that we know to follow the pre-change
distribution, and there is a large amount of reference data is available.
In offline change-point detection, the number of samples is fixed and our goal is to
detect the existence of a change-point τ , such that before the change-point, the samples
are i.i.d. with a distribution P , and after the change-point, the samples are i.i.d. with a
different distribution Q. The location τ where the change-point occurs is unknown. In
other words, we are concerned with testing the null hypothesis
H0 : xi ∼ P,

i = 1, . . . , t,

against the single change-point alternative
H1 : ∃1 ≤ τ < t

(
Q, i > τ
xi ∼
P, otherwise.

Note that we are interested in the case of a sustained change: before the change, all samples
follow one distribution, and after the change, all samples follow another distribution and
never switch back. In online change-point detection, the number of samples is not fixed,
and we keep getting new samples; the goal is to detect the emergence of a change-point.
In various change-point detection settings, the number of post-change samples is small, but
the number of reference samples is large. Therefore, when constructing MMD statistics over
blocks, we will use a common post-change block, and multiple disjoint pre-change reference
blocks.
3.1 Offline change-point detection
For each possible change location τ , the post-change block Y consists of the most recent
samples indexed from τ to t. Since we do not know the change-point location, we consider
all possible meaningful τ values. This corresponds to considering a range of post-change
block sizes B ranging from 2 (i.e., the most recent 2 samples are post-change samples; we
exclude B = 1 because the corresponding procedure is very unstable) to Bmax .
The detection statistic is constructed as follows, as illustrated in Fig. 1(a). Data are
split into (N + 1) blocks of size Bmax . Treat the last block as the post-change block, and
the remaining N blocks as reference data. Then we sample data from each block to form
smaller sub-blocks of various size B, B ∈ [2, Bmax ]. The reference blocks are denoted as
(B)
Xi , i = 1, . . . , N , and the test block as Y (B) . We compute MMD2u for each reference
6

Scan B-Statistic for Kernel Change-Point Detection

Block containing
potential change
point

Pool of reference data
B
MMD2u XN max , Y (Bmax

……

B

XN max
Bmax

……

Bmax

X3

Bmax

Bmax

X2

Bmax

Bmax

X1

Bmax

Y

MMD2u 𝑋𝑖

……

Bmax

Bmax

B

B

B

, 𝑌 (𝐵0,𝑡

Pool of reference data

2

2

……

B

2

Y (B0,t
𝑡

time

𝐵0
B0 ,t

Pool of reference data

Xi

2

(a): offline

Y (B0,t+1

𝑡+1

sample

2

Block containing
potential change
point

sample

𝑡 time

Xi

B

𝐵0 ,𝑡

B0 ,t+1

(b): sequential

Figure 1: Illustration of (a) offline change-point detection: data are initially split into blocks
of size Bmax ; then one samples from each initial block to form smaller blocks of various size B,
B = 2, . . . , Bmax ; (b) online change-point detection: the most recent B0 samples form the test
block, which is constantly updated with new data; the reference blocks of the same size B0 are
sampled from the reference pool of data.
sub-block with respect to the common post-change block, and take an average:
N
1 X
(B)
ZB =
MMD2u (Xi , Y (B) ).
N

(3)

i=1

Since the estimator MMD2u is unbiased, under the null hypothesis P = Q, E[ZB ] = 0.
Let Var[ZB ] denote the variance of ZB under the null. The variance of ZB depends on the
block size B and the number of blocks N . Hence, we divide each ZB by their standard
deviation
0
ZB
= ZB /(Var[ZB ])1/2 ,
and take the maximum over all B to form the offline B-statistic. The variance Var[ZB ] is
given in Lemma 1, which facilitates its computation. A change-point is detected whenever
the offline B-statistic exceeds a pre-specified threshold b:
max

2≤B≤Bmax

0
ZB
> b. {offline change-point detection}

(4)

3.2 Online change-point detection
In the online setting, new samples arrive in a streaming fashion and we constantly test
whether the incoming samples come from a different distribution. To reduce computational
burden, in the online setting, we fix the block-size and adopt a sliding window approach.
The detection statistic is constructed as follows, as illustrated in Fig. 1(b). At each
time t, we treat the most recent B0 samples as the post-change block. Since in online
change-point detection, we want to detect the change as quickly as possible, and hence
typically there are not many post-change samples. There is large amount of reference data.
To utilize data efficiently, we compare a common test block consisting of the most recent
samples, with N different reference blocks. The reference blocks are formed by taking
7

Li, Xie, Dai and Song

N B0 samples without replacement from the reference pool. Compute MMD2u between each
reference block with respect to the post-change block, and take an average:
ZB0 ,t =

N
1 X
(B ,t)
MMD2u (Xi 0 , Y (B0 ,t) ),
N

(5)

i=1

(B ,t)

where B0 is the fixed block-size, Xi 0 is the i-th reference block at time t, and Y (B0 ,t) is
the the post-change block at time t. When there is a new sample, we append it to the postchange block and purge the oldest sample. We show later that this construction allows for
an explicit characterization of the false-alarm rate. We divide each statistic by its standard
deviation to form the online B-statistic:
0
ZB
= ZB0 ,t /(Var[ZB0 ,t ])1/2 .
0 ,t

The calculation of Var[ZB0 ,t ] can also be achieved using Lemma 1. The online changepoint detection procedure is a stopping time: an alarm is fired whenever the standardized
detection statistic exceeds a pre-specified threshold b > 0:
0
T = inf{t : ZB
> b}.
0 ,t

{online change-point detection}

(6)

The online B-statistic can be computed efficiently. Note that the variance of the ZB0 ,t
only depends on the block size B0 but independent of t. Hence, it can be pre-computed.
Moreover, there is a simple recursive implementation of the online B statistic, with more
details specified in Appendix A.
3.3 Analytical expression for Var[ZB ]
We obtain an analytical expression for Var[ZB ], which is useful when forming the detection
statistic in (4) and (6).
Lemma 1 (Variance of ZB under the null) Given a fixed block size B ≥ 2 and number
of blocks N , under the null hypothesis,
 −1 



B
1
N −1
2
0
0
0
0
00 000
0
Var[ZB ] =
E[h (x, x , y, y )] +
Cov h(x, x , y, y ), h(x , x , y, y ) , (7)
2
N
N
where x, x0 , x00 , x000 , y, and y 0 are i.i.d. random variables with the null distribution P .
The lemma is proved by making a connection between MMD2u and U -statistic (Serfling,
1980) and utilizing the properties of U -statistic. Detail proof is provided in Appendix B.
We verified numerically that Lemma 1 is quite accurate, as shown in Fig. 2.
Lemma 1 provides a much simpler way to estimate Var[ZB ] than direct Monte Carlo
simulation. When using direct Monte Carlo simulation, to generate 10000 instances of ZB ,
one needs a total of 10000(N + 1)B samples. Each ZB needs (N + 1)B samples and is
computationally expensive to evaluate. Using Lemma 1, one only needs to evaluate terms
such as E[h2 (x, x0 , y, y 0 )] and Cov[h(x, x0 , y, y 0 ), h(x00 , x000 , y, y 0 )], which requires much fewer
samples to achieve a good estimation accuracy.
8

B=2

Percent Error

Percent Error

Scan B-Statistic for Kernel Change-Point Detection

0.4
0.2
0
0

2000
4000
6000
8000
Amount of Reference Data

B=200
0.4
0.2
0
0

10000

(e): B = 2, % error for estimated variance

2000
4000
6000
8000
Amount of Reference Data

10000

(f): B = 200, % error for estimated variance

Figure 2: Accuracy of Lemma 1 in estimating the variance of ZB when B = 2 and B = 200.
Figures show the percentage difference between the estimate from Lemma 1 relative to that
obtained from direct Monte Carlo. We form 10000 instances of ZB using data sampled from
the null distribution N (0, I20 ) Here, “Amount of reference data” denotes the size of pool to
sample the six-tuple (x, x0 , x00 , x000 , y, y 00 ) to evaluate the quantity in (7).

3.4 Examples of detection statistics
Below, we present a few examples to demonstrate the versatility of scan B-statistics.
Gaussian to Gaussian mixture. In Fig. 3(a), P = N (0, I2 ), Q is a mixture Gaussians:
0.3N (0, I2 ) + 0.7N (0, 0.1I2 ), and τ = 250. The online procedure stops at time 270, meaning
the change is detected with a small delay of 20 unit time.

After Change

Before Change

30
20

50

100
Time

150

200

0
100

200

Time

300

400

500

(a): Gaussian to GMM, τ = 250

b=3.83

400
600
Time

800

1000

Bandw=Med
Bandw=100Med
Bandw=0.1Med

50
b=3.55

2
0
0

200

100

6
4

0

−100

Statistic

Statistic

Statistic

0
0
8

5 b=3.73

50

−50

10

10

0

100

40

Seismic Signal

Index of Adjacency
Matrix Entry

Sequence of graphs. In Fig. 3(b), we consider detecting the emergence of a community
inside the network modeled using a stochastic block model, which has been considered
in (Maragoni-Simonsen and Xie, 2015). Assume that before the change, each sample is
a realization of an Erdős-Rényi random graph, with the probability of forming an edge
p0 = 0.1 uniform across the graph. After the change, a “community” emerges, which is
a subset of nodes, where the edges are formed in between these nodes with much higher
probability p1 = 0.3. This models a community where the members inside the community
interact more often. Our online procedure stops at time 102, meaning the change is detected
with a small delay of 2 unit times.

0

50

100
Time

150

(b) Graphs, τ = 100

200

200

400
600
Time

800

1000

(c): Real seismic signal

Figure 3: Examples of offline and online B-statistics with N = 5. All thresholds (red lines)
are theoretical values obtained from Theorem 2 (offline) and Theorem 9 (online).
9

Li, Xie, Dai and Song

Real seismic signal; effect of kernel bandwidth. In Fig. 3(c), we consider a segment
of real seismic signal that contains a change-point. Using the seismic signal, we also illustrate the effect of kernel bandwidth. The performance of the B-statistic depends
on the

kernel bandwidth. For Gaussian RBF kernel k(Y, Y 0 ) = exp −kY − Y 0 k2 /2σ 2 , the kernel
bandwidth σ > 0 is typically chosen using a “median trick” (Schölkopf and Smola, 2001),
where σ is set to be the median of the pairwise distances between data points. The accuracy
of the median trick is justified by Ramdas et al. (2015). The offline B-statistic can detect
the seismic event reliably across different choices of the kernel bandwidth.

4. Theoretical approximation of significance level for offline B-statistic
In the offline setting, the choice of the threshold b involves a tradeoff between two standard
performance metrics: (1) significance level (SL), which is the probability that the B-statistic
exceeds the threshold b under the null hypothesis (i.e., when there is no change); and (2)
power, which is the probability of the statistic exceeds the threshold under the alternative
hypothesis.
We present an accurate approximation to the SL of the offline B-statistic, assuming the
detection threshold b tends to infinity and the number of blocks N is fixed. The following
theorem is our main result.
Theorem 2 (SL of offline B-statistic) When b → ∞, and Bmax → ∞, with b/(Bmax )1/2
held a fixed positive constant, the significance level of the offline B-statistic defined in (4)
is given by
s
!


B
max
X
1 2
2B
−
1
(2B
−
1)
0
√
> b = be− 2 b ·
P
max ZB
ν b
· [1 + o(1)] , (8)
2≤B≤Bmax
B(B − 1)
2 2πB(B − 1)
B=2

where the special function
ν(µ) ≈

(2/µ)(Φ(µ/2) − 0.5)
,
(µ/2)Φ(µ/2) + φ(µ/2)

(9)

φ(x) and Φ(x) are the probability density function and the cumulative distribution function
of the standard normal distribution, respectively.
We can show numerically that the approximation (8) is accurate even in the nonasymptotic regime. Hence, it is quite useful in choosing a threshold b given a target SL,
without having to resort to the onerous numerical simulations. First, consider synthetic
data that are i.i.d. normal P = N (0, I20 ). We set Bmax to be 50, 100, 150, and in each
case, N = 5. We compare the thresholds obtained by the theoretical approximation and by
simulation, for a prescribed SL α. To obtain threshold by simulation, we generate Monte
Carlo trials to for offline B-statistics and find the (1 − α)-quantile as the estimated threshold. Table 1 shows that for various choices of Bmax , the thresholds predicted by Theorem 2
match quite well with those obtained by simulation. The accuracy can be further improved
for smaller α values, by skewness correction in Section 7.
Next, we consider real speech data from the CENSREC-1-C dataset (more details in
Section 8). Here, the null distribution P corresponds to the unknown distribution of the
10

Scan B-Statistic for Kernel Change-Point Detection

background signal, and we are interested in detecting the onset of speech signals. This
case is more challenging because the true distribution can be arbitrary. In the dataset,
there are 3000 reference samples. We bootstrap these reference samples to generate 10000
re-samples to estimate the tail of the detection statistic. Table 2 demonstrates that the
thresholds predicted by theory are quite accurate. Note that in this case, the accuracy
improves significantly by skewness correction.
Table 1: Thresholds for the offline B-statistics using synthetic data, obtained by simulation, theory (Theorem 2), and theory with Skewness Correction (Section 7), respectively,
for various SL values α.
α
0.10
0.05
0.01

b (sim)
2.41
2.77
3.54

Bmax = 50
b (theory)
2.38
2.67
3.23

b (SC)
2.57
2.97
3.64

b (sim)
2.43
2.76
3.47

Bmax = 100
b (theory)
2.50
2.78
3.32

b (SC)
2.76
3.17
3.82

b (sim)
2.53
2.97
3.64

Bmax = 150
b (theory)
2.56
2.83
3.37

b (SC)
2.89
3.22
3.89

Table 2: Thresholds for the offline B-statistics using real speech data, obtained by simulation, theory (Theorem 2), and theory with Skewness Correction (Section 7), respectively,
for various SL values α.
α
0.10
0.05
0.01

b (boot)
2.96
3.62
4.85

Bmax = 50
b (theory)
2.38
2.67
3.23

b (SC)
3.23
3.68
4.61

b (boot)
3.16
3.82
5.20

Bmax = 100
b (theory)
2.50
2.78
3.32

b (SC)
3.59
4.06
5.03

b (boot)
3.21
3.86
5.42

Bmax = 150
b (theory)
2.56
2.83
3.37

b (SC)
3.94
4.43
5.45

4.1 Proof of main theorem
The key to our approximation is to quantify the tail probability of the detection statistic
when the threshold is large. In the offline setting, to quantify the significance level, we are
0 > b} under the null distribution. The
interested in find the probability P {max2≤B≤Bmax ZB
0 }
main assumption here is that the collection of random variables {ZB
B=2,...,Bmax , form a
Gaussian random field. This means the finite-dimensional joint distributions of the collection of random variables are all Gaussian, and they are completely specified by its mean
and covariance functions, which we characterize below.
0 ] is zero since they are constructed
Under the null distribution, the expectation E[ZB
using the unbiased MMD estimator. The covariance under the null distribution is given by
the following lemma:
Lemma 3 (Covariance structure of Z-statistic in the offline setting) Under the null
distribution, the covariance between two offline B-statistics is given by
s  


u v
u∨v
ru,v = Cov Zu0 , Zv0 =
, 2 ≤ u, v ≤ Bmax ,
(10)
2
2
2
where u ∨ v = max{u, v}.
11

Li, Xie, Dai and Song

The proof can be found in Appendix B.2.
The proof to Theorem 2 is mainly based on the change-of-measure technique. The idea
is similar to importance sampling, but it is much more sophisticated. Under null distri0 > b} for a large threshold b is rare (so
bution, the boundary exceeding event {maxB ZB
that false alarm remains low). Since quantify such small probability is hard, we considers an alternative probability measure under which the event happens with much higher
probability. Thus, under the new measure, one can use the local central limit theorem to
a obtain an analytical expression. Finally, the original small probability will be related to
the probability under the alternative measure using the Mill’s ratio (Yakir, 2013).
Below, we present the main steps of the proof, including (1) exponential tilting; (2)
change-of-measure by the likelihood identity; (3) establish properties of the local field and
the global term; (4) perform asymptotic approximation using the localization theorem by
showing that the “global” log likelihood and the “local process” are asymptotically independent. Finally, we collect terms together to obtain the result.
4.1.1 Step One: Exponential tilting
We first introduce exponential tilting, which creates a family of distributions that is related
0 be
to the original distribution of B-statistic. Let the log moment generating function of ZB
0

ψ(θ) = log E[eθZB ].

(11)

 0
dPB = exp θZB
− ψ(θ) dP,

(12)

Define a family of new measures

0 under the null distribution P ,
where P represents the original probability measure of ZB
PB is the new measure after the transformation, and θ parameterizes the family of new
measures. Note that the new measures take the form of exponential family, and θ can be
viewed as the natural parameter. It can be verified that PB is indeed a probability measure
since
Z
Z
 0
dPB = exp θZB
(x) − ψ(θ) dP = 1.

It can be shown ψ̇(θ) is equal to the expected value under the tilted measure, since
0

0 eθZB ]
E[ZB
0
0 θZB
ψ̇(θ) =
= E[ZB
e −ψ(θ) ] =
0
θZ
E[e B ]

Z
xe

θx−ψ(θ)

Z
dP =

xdPB .

Similarly, one can show that ψ̈(θ) equals to the variance under the tilted measure.
0 has zero mean and unit variance, since
Recall that, under the null distribution, ZB
it is constructed using unbiased MMD estimators and has been divided by the standard
0 is a standard Gaussian random variable with zero mean and
deviation. Now we assume ZB
unit variance. Under this assumption, the corresponding log moment generating function
is given by ψ(θ) = θ2 /2.
One still has the freedom to select the value of θ. We will set the value of θ such that
the mean under the tilted measure is equal to a given threshold b. This means that the new
measure peaks at the threshold b (this essentially enables us to use the local central limit
12

Scan B-Statistic for Kernel Change-Point Detection

theorem later on). This can be done by choosing θ such that ψ̇(θ) = b. Since ψ(θ) = θ2 /2
for standard normal random variable, this means θ = b. Note that the solution θ does
not depend on B. Hence, we can set the mean under the transformed measure to b, by
uniformly choosing θ = b for any B. Given such a choice, the transformed measure is
0 (x) − b2 /2 dP. We also define, for each B, the log-likelihood ratio
given by dPB = exp bZB
log(dPB /dP) of the form
0
`B = bZB
− b2 /2.

(13)

0 to a likelihood ratio, even if Z 0
This way, we have associated the detection statistic ZB
B
itself does not come out of a likelihood ratio.
0 under the new measure has the same unit variance
The following lemma shows that ZB
and its mean has been shifted to b. This key fact will lead to the desired exponential tail.

Lemma 4 (Mean and variance under tilted measure) Define EB and VarB as the
expectation and variance under the transformed measures
EB [U ] = E[U e`B ],

(14)

VarB [U ] = E[U 2 e`B ] − E2B [U ].

(15)

0 ] = b, and Var [Z 0 ] = 1.
We have EB [ZB
B B
0 ] = ψ̇(b) = b by construction. To show Var [Z 0 ] = 1, note that
Proof First, EB [ZB
B B
0
0
0 )2 ebZB
log E[ebZB ] = b2 /2. Taking the derivative of ψ(θ) with respect to b twice gives E[(ZB
]=
2 /2
2 /2
0 −ψ(b)
b
2
b
0
2
0
2
θZ
2
0
0
2
B
e
+b e
. Hence, EB [(ZB ) ] = E[(ZB ) e
] = 1+b , and VarB [ZB ] = EB [(ZB ) ]−
b2 = 1.

4.1.2 Step Two: Change-of-measure by likelihood ratio identity
Now we are ready to analyze the tail probability. The basic idea is to convert the original
problem of finding the small probability that the maximum of a random field exceeds a
large threshold, to another problem: finding an alternative measure under which the event
happens with a much higher probability.
Here, the alternative measure will be a mixture of simple exponential tilted measures.
Define the maximum and the sum for likelihood ratio differences relative to a particular
parameter value B:
MB =

max

s∈{2,...,Bmax }

e`s −`B ,

SB =

X

e`s −`B .

s∈{2,...,Bmax }

Also define a re-centered likelihood ratio, which we call the global term
0
`˜B = b(ZB
− b).

13

(16)

Li, Xie, Dai and Song

With the definitions above and the log-likelihood ratios `B in (13), we have the following
likelihood ratio identity:

P

max

2≤u≤Bmax



0
Zu > b = E 1;

#
"P

Bmax `B
e
; max Zu0 > b
max Zu0 > b = E PB=2
Bmax `s 2≤u≤B
2≤u≤Bmax
max
e
| s=2
{z }
=1

B
max
X




Bmax
1
(14) X
0
0
=
E P `s ; max Zu > b =
EB P `s ; max Zu > b
2≤u≤Bmax
2≤u≤Bmax
se
se
B=2
B=2


B
max
X
MB −(`˜B +log MB ) ˜
−b2 /2
=e
EB
e
; `B + log MB ≥ 0
SB


e`B

(17)

B=2

where an intermediate step is done by changing the measure to PB , and the last equality
can be verified by simple algebra. Recall our notation EB [A; B] = EB [A1{B}] for a random
quantity A and event B; 1 denotes an indicator function.
In a nutshell, the last equation in (17) converts the tail probability to a product of
2
two terms: a deterministic term e−b /2 associated with the large deviation rate, and a sum
of conditional expectations under the transformed measures. A close examination of the
conditional expectations of the form EB [· · · ; [· · · ] ≥ 0] reveals that it involves a product of
the ratio MB /SB , and an exponential function that depends on `˜B , which plays the role of
a weight. Under the new measure PB , `˜B has zero mean and variance equal to b2 (shown
below in Lemma 6) and it dominates the other term log MB and, hence, the probability of
exceeding zero will happen with much higher probability. Next, we characterize the limiting
ratio and the other factors precisely, by the localization theorem.
4.1.3 Step Three: Establish properties of local field and global term
In (17), our target probability has been decomposed into terms that only depend on (i) the
local field {`s −`B }, 2 ≤ s ≤ Bmax , which are the differences between the log-likelihood ratio
with parameter B and with other parameter values s, and (ii) the global term `˜B , which is
the centered and scaled likelihood ratio with parameter B. We need to first establish some
useful properties of the local field and global term when using the localization theorem. We
will eventually show that the local field and the global term are asymptotically independent.
The following property for the global term can be derived from Lemma 4. The result
shows that under the tilted measure, the global term `˜B is zero mean for any B, with
variance diverging with b.
Lemma 5 (Global term for offline B-statistic) The mean and variance of the global
0 − b), for 2 ≤ B ≤ B
term `˜B = b(ZB
max , are given by
EB [`˜B ] = 0,

VarB [`˜B ] = b2 ,

(18)

0 is approximately normal, the local field ` − ` (or equivalently
Since we assume ZB
s
B
0 )) and the global term `˜ (or equivalently b(Z 0 − b)) are also approximately
b(Zs0 − ZB
B
B
normally distributed.

14

Scan B-Statistic for Kernel Change-Point Detection

Lemma 6 (Local field for offline B-statistic) The mean and variance of the local field
{`s − `B }, for |s − B| = 0, 1, 2, . . ., are given by
EB [`s − `B ] = −b2 (1 − rs,B ),

VarB [`s − `B ] = 2b2 (1 − rs,B ),

with rs,B defined in (10). For any s1 and s2 , the covariance between two local field terms is
given by
CovB (`s1 − `B , `s2 − `B ) = b2 (1 + rs1 ,s2 − rs1 ,B − rs2 ,B ) .
0 ), E [Z 0 ] = b, Var [Z 0 ] = 1. Moreover, due
Proof Note that `s − `B = b(Zs0 − ZB
B B
B B
0 , we have the following decomposition E [` − ` ] =
to the normal assumption of ZB
B s
B
0 )] = E [b(r
0
2
1/2 W − Z 0 )] = −b2 (1 − r
EB [b(Zs0 − ZB
),
where
W is
B
s,B ZB + (1 − rs,B )
s,B
B
0 , representing residual of regression.
a zero-mean random variable and independent of ZB
The variance and covariance can be found using similar decompositions.

Remark 7 (Consequence of Lemma 6) From the expression of the covariance in (10),
we have that for s − B > 0,
rs,B = [1 + (s − B)/B]−1/2 [1 + (s − B)/(B − 1))]−1/2 ,
and for s − B < 0,
rs,B = [1 + (s − B)/B]1/2 [1 + (s − B)/(B − 1)]1/2 .
Consequently,
1. When |s − B| → ∞, rs,B → 0. Therefore, when |s − B| → ∞, EB [`s − `B ] converges
to −b2 and VarB [`s − `B ] converges 2b2 .
2. When |s − B| is small, assume s = B + j, j = 0, ±1, ±2, . . .. Perform the Taylor
expansion of rB+j,B around 0, we have
rB+j,B = 1 −

1 2B − 1
|j| + o(|j|).
2 B(B − 1)

(19)

Define
µ = b{(2B − 1)/[B(B − 1)]}1/2 .

(20)

Note that µ depends on the threshold as well as B, the block size parameter. Using
(19), we have
lim EB [`B+j − `B ] = −

|j|→0

µ2
|j|,
2

lim VarB [`B+j − `B ] = µ2 |j|,

|j|→0

lim

|j1 |→0,|j2 |→0

CovB (`B+j1 − `B , `B+j2 − `B ) = µ2 (|j1 | ∧ |j2 |).

15

Li, Xie, Dai and Song

Therefore, when |j| is small (i.e., in the neighborhood of zero), we can approximate
the local field using a two-sided Gaussian random walk with drift µ2 /2 and variance
of the increment being µ2 :
d

`B+j − `B = µ

|j|
X

ϑi − µ2 j/2,

j = ±1, ±2, . . .

(21)

i=1

where ϑi are i.i.d. standard normal random variables.
4.1.4 Step Four: Approximation using localization theorem
The remaining work is to compute the conditional expectations EB [· · · ; (· · · ) ≥ 0] for each
B in (17). In the following, we drop the subscript B in EB for simplicity. We assume
b → ∞, Bmax → ∞, and b2 /Bmax is held to a fixed positive constant. Introduce an abstract
index κ and let κ = b2 ; this choice is because κ1/2 is the multiplicative factor that balances
the rate of convergence of the global term under the transformed measure. Typically, κ is
equal to the variance of the global term, which is equal to b2 here.
Using a powerful localization theorem (see Theorem 3.1 in (Siegmund et al., 2010) or
Theorem 5.2 in (Yakir, 2013)), we can obtain the limit for each term in the summand of
(17), rewritten as (by changing the index to κ)



Mκ −(`˜κ +log Mκ ) ˜
E
e
; `κ + log Mκ ≥ 0 ,
Sκ

(22)
1

when κ → ∞. Basically, the localization theorem states that (22) scaled by κ 2 converges
to a limit under mild conditions when κ → ∞.
The statement of the theorem involves a local σ-algebra denoted as Fbκ :
Fbκ = σ{`s − `B : |s − B| ≤ g(κ)},

(23)

where a function g(κ) specifies the size of the local region. The choice of g(κ) is critical and
it guarantees subsequent convergence. Following the analysis for scan statistics in (Yakir,
2013), we choose g(κ) = cb−2 for some large constant c. This local σ-field is asymptotically
independent of `˜k , and it carries all information needed to construct the local field.
cκ and Sbκ as the maximization and summation restricted to a smaller subset of
Define M
parameter values {s : |s − B| ≤ g(κ)}, and they are measurable with respect to Fbκ . Note
cκ and Sbκ serve as approximations to Mκ and Sκ . In the limit, the local random field
that M
cκ /Sbκ ] converges to a limit that can
converges to a Gaussian random field, and the ratio E[M
be determined with the parameters of the Gaussian random field.
The localization theorem (Theorem 5.1 in Siegmund et al. (2010) and Sec. 3.4 in Yakir
(2013)) consists of five conditions as follows.
Theorem 8 (Localization theorem) Given  > 0, if for all large κ, all following conditions hold
cκ ≤ Sbκ < ∞ hold in probability one.
I. Both 0 < Mκ ≤ Sκ < ∞ and 0 < M
16

Scan B-Statistic for Kernel Change-Point Detection

cκ | > } ∪ {|Sbκ /Sκ − 1| > }. For some 0 < δ that does
II. Denote Ac = {| log Mκ − log M
not depend on :
h
i
cκ ∈ x + (0, δ]} ∩ {|m̂| ≤ g(κ)} ≤ κ−1/2 ,
max P Ac ∩ {`˜κ + log M
|x|≤3g(κ)

cκ , g(κ)} − log(1 − ).
where m̂κ = min{log M
cκ /Sbκ ] converges to a finite and positive limit denoted by E[M/S].
III. E[M
IV. There exist µκ ∈ R and σκ ∈ R+ such that for every 0 < 0 , δ, for any event E ∈ Fbκ
and for all large enough κ
sup
|x|≤κ1/2

δ µ
κ1/2 P(`˜κ ∈ x + (0, δ], E) − φ
P(E) ≤ 0 .
σ
σ

cκ | > κ1/2 ) and P(log Mκ − log M
cκ < −) are all
V. P(| log Mκ | > κ1/2 ), P(| log M
−1/2
o(κ
).
Then
lim κ

κ→∞

1/2


µ
Mκ −[`˜κ +log Mκ ] ˜
e
; `κ + log Mκ ≥ 0 = σ −1 φ
E[M/S],
E
Sκ
σ


(24)

where φ(·) is the density of the standard normal distribution.
Intuitively, the localization theorem says the following. To find the desired limit of
(22) as κ → ∞, one first approximates Mκ and Sκ by their localized versions, which are
obtained by restricting the maximization and summation in a neighborhood of parameter
values. Then one can show that the localized ratio Mκ /Sκ is asymptotic independent of
the global term `˜κ as κ → ∞. The asymptotic analysis is then done on the local field and
the global term separately. The expected value of the localized ratio E[Mκ /Sκ ] converges
to a constant independent of κ, and the limiting conditional distribution of `˜κ can be found
using the local central limit theorem. Thus, one can calculate the remaining conditional
expectation involving `˜κ .
4.1.5 Step Last: Checking conditions
Let us now verify the validity of the conditions in our setting. First, Condition I is met since
for Gaussian random variables, Mκ > 0, Sκ > 0 with probability 1, and the maximization
of a collection of non-negative numbers is smaller or equal to the summation. Similar
cκ > 0 and Sbκ > 0 when the maximization and
arguments hold for their counterparts M
summation are over a smaller set.
cκ and Sbκ are good approximations
Condition II describes that the localized versions M
of Mκ and Sκ when κ is sufficiently large, for properly defined Fbκ . In Section 3.4.4 of
Yakir (2013), the corresponding Condition II has been rigorously checked, assuming a local
region defined in the same form of our local region and assuming Gaussian random field.
Thus, checking Condition II for our case will follow the same steps, using the properties
established in Section 4.1.3. We omit the details here.
17

Li, Xie, Dai and Song

Condition III is checked by applying distributional approximations to the localized version of Mκ /Sκ . We can show that the expectation of the ratio E[M̂κ /Ŝκ ] converges to a
finite and positive limit denoted by E[M/S], which does not depend on κ. Since the increment `B+j − `B has negative mean as shown in Lemma 6, the values of Mκ and Sκ will be
determined by values j close to 0, so is the ratio Mκ /Sκ . This implies, a relatively small
local region centered on B will be sufficient. From Remark 7, the local field when the index
is close to the shifted measure parameter B can be approximated as a two-sided Gaussian
random walk with drift −µ2 /2 and variance µ2 (with µ defined in (20)), which is denoted
as W (µ2 j) below. Therefore, we have that with high probability,
2

"

max|j|≤cb−2 eW (µ j)
E[M̂κ /Ŝκ ] = E P
W (µ2 j)
|j|≤cb−2 e

#

When c → ∞, it approaches to a limit known as the Mill’s ratio:
#
2
max|j| eW (µ j)
E[M/S] = E P W (µ2 j) ,
|j| e
"

with maximization and summation extending to the entire collection of negative and positive
integers. The Mill’s ratio is related to the Laplace transform of the overshoot of the maxima
of Gaussian random field over a threshold b, and an expression has been obtained based
on nonlinear renewal theory
1985) and Chapter 2.2 of the book (Yakir,
P (see, (Siegmund,
1/2 µ/2)). An easier numerical evaluation is given by
Φ(−j
2013)): E[M/S] = exp(−2 ∞
j=1
E[M/S] ≈ (µ2 /2)v(µ) for a special function ν(µ) defined in (9).
Condition IV can be checked via a local multivariate central limit theorem that is local
in one component and non-local in others (Theorem 5.3 in Yakir (2013)). The theorem
says the following: assuming ξi are independent, identically distributed random vector of
dimension d+1. Assume the mean of each vector is zero, and variance of the first component
converges to a finite σ, the covariance matrix of the last d components converges a finite
matrix Σ, and the correlation between these components and the first oneP
converges to
zero (hence, the overall covariance matrix is block-diagonal). Define Sγ = γi=1 ξi,1 and
a d dimensional vector with element hγ,j = γ −1/2 ξi,j , for 1 ≤ j ≤ d. Then under mild
conditions,
l−u
lim γ 1/2 P(Sγ ∈ [l, u], hγ ∈ A) =
P(h ∈ A)
(25)
γ→∞
(2π)1/2 σ
for any interval [l, u] and an arbitrary set A.
Our setting matches exactly to the above distribution when we set the global term as the
first component and local field as the remaining components. Using the properties in Section
4.1.3, we have shown the finite mean and variance (covariance) of the global and local field
terms. We only need to show the global term and the local fields are independent of each
other asymptotically. It suffices to prove that the conditional covariance of {`B+j − `B }
given `˜B converges to the unconditional covariance, and the conditional means converges
to the unconditional one. With a slight abuse of notation, r1 = rB+j1 ,B and r2 = rB+j2 ,B .
0 (which is proportional
Using the linear regression decomposition, when conditioning on ZB
18

Scan B-Statistic for Kernel Change-Point Detection

to `˜B ), the two local field terms are independent of each other:
0
0
0
0
0
Cov(b(ZB+j
− ZB
), b(ZB+j
− ZB
)|ZB
)
1
2
0
0
0
0
0
= Cov(b(r1 ZB
+ (1 − r12 )1/2 W1 − ZB
), b(r2 ZB
+ (1 − r22 )1/2 W2 − ZB
)|ZB
) = 0.

where W1 and W2 are two mutually independent zero-mean random variables that represent
0 ). On the other hand, using the
the regression residuals (they are also independent of ZB
same decomposition, we can show that without conditioning, the covariance is given by
0
0
0
0
Cov(b(ZB+j
− ZB
), b(ZB+j
− ZB
)) = b2 (1 − r1 )(1 − r2 ).
1
2

Hence, when b → ∞, due to the property of local field in equation (19), for |j1 | ≤ cb−2 , |j2 | ≤
cb−2 , the unconditioned covariance converges to zero, which is equal to the conditioned
0
0 } conditioning
covariance. Similarly, we can show that the conditional means of {ZB+j
− ZB
0
on ZB converges to the unconditional ones.
Now we invoke the local central limit theorem. Since the density of the global term
`˜B is approximately normal, we can calculate a desired form of the probability. From
(18), the variance of the global term increases with b. The density of `˜B can be uniformly
approximated by 1/(2πb2 )1/2 within a small region around the origin |x| ≤ 3(4/ + 1 + ) log b
(Yakir, 2013). Such an approximation also holds for `˜B − x given any value x that is not too
large. Furthermore, notice that log M̂κ is very close to 0 and therefore is negligible; this is
because e`s −`B should attain its maximal value when |s − B| close to 0 as analyzed before.
Let µκ = EB [`˜κ /b] = 0 and σκ2 = VarB [`˜κ /b] = 1. When κ = b2 → ∞, using local central
limit theorem (25), we have that
 


µκ
δ
1/2
˜
φ
.
(26)
κ P `κ ∈ x − log M̂κ + (0, δ] →
σκ
σκ
Condition V is checked as follows. Note that the terms inside the maximization are likelihood ratios with unit expectation since EB [exp(`s −`B )] = 1. This is because from Lemma
6, under the transformed measure, `s −`B is normal random variable with mean −b2 (1−rs,B )
and variance 2b2 (1 − rs,B ). Hence, its moment generating function EB [exp(θ(`s − `B ))] = 1
for any θ. Thus, exp(`s − `B ) is a martingale and by a standard martingale inequality,
P(log Mκ > κ1/2 ) ≤ exp(−κ1/2 ). Then using a similar argument as in Siegmund et al.
cκ is an approximation to Mκ .
(2010), one can show the other two inequalities, since M
Finally, since all conditions are met, we can now apply the localization theorem for
b → ∞ and put things together to obtain


MB −[`˜B +log MB ] ˜
µ2
1
EB
e
; `B + log MB ≥ 0 =
ν(µ) √
(1 + o(1)).
(27)
SB
2
2πb2
Substitute (27) back to the likelihood ratio identity (17), we arrive at the approximation in
Theorem 2.

5. Theoretical approximation of ARL for online B-statistic
In the online setting, two commonly used performance metrics are (see, e.g., Xie and Siegmund (2013)): (1) the average run length (ARL), which is the expected time before incorrectly announcing a change of distribution when none has occurred; (2) the expected
19

Li, Xie, Dai and Song

detection delay (EDD), which is the expected time to fire an alarm when a change occurs
immediately at τ = 0. The EDD considers the worst case and provides an upper bound on
the expected delay to detect a change-point when the change occurs later in the sequence
of observations.
We present an accurate approximation to the ARL of online B-statistics. The approximation is quite useful in setting the threshold. As a result, given a target ARL, one can
determine the corresponding threshold value b from the analytical approximation, avoiding
the more expensive numerical simulations. Our main result in the following
Theorem 9 (ARL in online B-statistic) Let B0 ≥ 2. When b → ∞, the average run
length (ARL) of the stopping time T defined in (6) is given by
2

eb /2
E[T ] =
·
b

(

(2B0 − 1)
√
·ν
2πB0 (B0 − 1)

s
b

2(2B0 − 1)
B0 (B0 − 1)

!)−1
· [1 + o(1)] .

(28)

We verify the accuracy of the approximation numerically, by comparing the thresholds
obtained by Theorem 9 with those obtained from Monte Carlo simulation. Consider several
cases of null distributions: the standard normal N (0, 1), exponential distribution with mean
1, a Erdős-Rényi random graph with 10 nodes and probability of 0.2 of forming random
edges, as well as Laplace distribution with zero mean and unit variance. The simulation
results are obtained from 5000 direct Monte Carlo trials. As shown in Fig. 4, the thresholds
predicted by Theorem 9 are quite accurate. Fig. 4 also demonstrated that theory is quite
accurate for various block sizes (especially for larger B0 ). However, we also note that theory
tends to underestimate the thresholds. This is especially pronounced for small B0 , e.g.,
B0 = 50. The accuracy of the theoretical results can be improved by skewness correction,
shown by black lines in Fig. 4, which are discussed later in Section 7.
2
Theorem 9 shows that ARL is O(eb ) and, hence, b is O((log ARL)1/2 ). Note that
EDD is typically on the order of b/∆ due to Wald’s identity (Siegmund, 1985), where ∆
is the Kullback-Leibler (KL) divergence between the null and alternative distributions (a
constant). Hence, given a desired ARL (typically on the order of 5000 or 10000), the error in
the estimated threshold will only be translated linearly to EDD. This is a blessing, since it
means typically a reasonably accurate b will cause little performance loss in EDD. Similarly,
2
Theorem 2 shows that SL is O(e−b ) and a similar argument can be made for the offline
case.
5.1 Proof sketch of Theorem 9
The method for approximating the ARL is related to what we used to analyze the tail
probability of the offline B-statistic. First, we need to establish the following lemma.
Lemma 10 (Asymptotic null distribution of T ) Under the null, when b → ∞, the
stopping time T defined in (6) is uniformly integrable and asymptotically exponentially distributed, i.e.,
|P{T ≥ m} − exp(−λ0 m)| → 0,
in the range where mλ0 is bounded away from 0.
20

Scan B-Statistic for Kernel Change-Point Detection

8
7

7
6

5

b

b

6

8

Theo. (Gauss. Appr.)
Theo. (Skew. Corr. for Gaussian(0,I))
Simu. for Gaussian(0,I)
Simu. for Exp(1)
Simu. for Random Graph (Node=10, p=0.2)
Simu. for Laplace(0,1)

5

4

4

3

3

2
0

0.2

0.4

0.6
ARL(104)

0.8

2
0

1

Theo. (Gauss. Appr.)
Theo. (Skew. Corr. for Gaussian(0,I))
Simu. for Gaussian(0,I)
Simu. for Exp(1)
Simu. for Random Graph (Node=10, p=0.2)
Simu. for Laplace(0,1)

(a): B0 = 50

0.2

0.4

0.6
ARL(104)

0.8

1

(b): B0 = 200

Figure 4: For a range of target ARL values, comparison of thresholds determined from
simulation, from Theorem 9, and from theory with the skewness correction (Section 7)
under various null distributions. Shaded areas represent standard deviations for skewnesscorrected thresholds.

The proof is provided in Appendix C. It is adapted from arguments in Siegmund and
Venkatraman (1995); Siegmund and Yakir (2008); Yakir (2009). The main idea is to show
that a random variable defined as the number of boundary cross events for detection statistic
over disjoint intervals, converges to Poisson random variable in the total variation norm.
This is resulted from the Poisson limit theorem (Theorem 1 in Arratia et al. (1989)) for
dependent samples.
Uing Lemma 10, we know for large m, P{T ≥ m} is a approximately 1 − exp(−λ0 m) ≈
mλ0 , and E{T } is equal to λ−1
0 asymptotically when b → ∞. So the remaining question is to
0
find the probability and the corresponding λ0 . Consider P{T ≤ m} = P{max2≤t≤m ZB
>
0 ,t
1 2

b}. Suppose m > B0 and log b  m  b−1 e 2 b . We will adopt a similar strategy to approximate this probability using the change-of-measure technique. Note that the covariance
structure for online and offline B-statistic are different, so there will be different drift parameters when we invoke the localization theorem. Using exponential tilting, we introduce
a likelihood ratio
0
ζt = bZB
− b2 /2.
0 ,t
Again using the likelihood ratio identity decomposition, we obtain

P

max

2≤t≤m

0
ZB
0 ,t


 0

m
X
Mt −[ζ̃t +log Mt ]
−b2 /2
0
e
; ζ̃t + log Mt ≥ 0 ,
>b =e
Et
St0
t=2

where
Mt0 = max eζs −ζt ,
2≤s≤m

St0 =

X

eζs −ζt ,

2≤s≤m

21

0
and ζ̃t = b(ZB
− b).
0 ,t

(29)

Li, Xie, Dai and Song

Hence, one can again apply the localization theorem to find the approximation when b → ∞,
and the only differences are in the definition of global and local field terms.
Lemma 11 (Local field of online B-statistic) The mean, variance, and covariance of
the local field {ζs − ζt } are given by
Et [ζs − ζt ] = −b2 (1 − rs,t ),

Vart [ζs − ζt ] = 2b2 (1 − rs,t ),

Covt (ζs1 − ζt , ζs2 − ζt ) = b2 (1 + rs1 ,s2 − rs1 ,t − rs2 ,t ) ,
where
ρs,t =

0
0
Cov(ZB
, ZB
)
0 ,s
0 ,t

=

(B0 −|t−s|)∨0
2

B0
2


.

(30)

The proof can be found in Appendix B.3. Note that when |t − s| is close to 0, Et [ζs − ζt ] is
close to 0. With an increasing |t − s|, Et [ζs − ζt ] decreases until |t − s| > B0 (when there
are no overlapping test data in the sliding block), then Et [ζs − ζt ] becomes −b2 . The values
of Mκ and Sκ will be determined by the values of |j| close to 0.
Now, again, we will use an argument based on Taylor expansion to find the drift term
of the local field. When |s − t| is close to 0, we can approximate {ζs − ζt } as a two-sided
random walk. Using Taylor expansion we have
ρt+j,t = 1 −

2B0 − 1
|j| + o(|j|).
B0 (B0 − 1)

(31)

Let λ = b[2(2B0 − 1)]/[B0 (B0 − 1)]1/2 . Hence, we can show that the mean, variance, and
covariance of the local field are approximately
lim Et [ζt+j − ζt ] = −

|j|→0

λ2
|j|,
2

lim Vart [ζt+j − ζt ] = λ2 |j|,

|j|→0

lim

|j1 |→0,|j2 |→0

Covt (ζt+j1 − ζt , ζt+j2 − ζt ) = λ2 (|j1 | ∧ |j2 |).

As a result, by invoking the localization theorem through a similar set of steps, we obtain
s
!
1 2
2(2B0 − 1)
be− 2 b (2B0 − 1)
·ν b
(1 + o(1)),
(32)
P {T ≤ m} = m · √
B0 (B0 − 1)
2π B0 (B0 − 1)
Matching this to above, we know λ0 is the factor that multiplies m and this leads to the
desired result.

6. Study of power and expected detection delay
In this section, we study the detection power and the expected detection delay of the offline
and online B-statistics, respectively, by comparing with classic methods.
22

Scan B-Statistic for Kernel Change-Point Detection

6.1 Offline change-point detection: Comparison with parametric statistics
We compare the offline B-statistic with two commonly used parametric test statistics:
the Hotelling’s T 2 and the generalized likelihood ratio (GLR) statistics. Assume samples
{x1 , x2 , . . . , xn }.
Hotelling’s T 2 statistic. For a hypothetical change-point location k, the statistic measures
the quadratic distance of the sample means in two segments [1, k] and [k + 1, t] in the
quadratic distance defined by the sample covariance matrix:
T 2 (k) =
where, x̄k =

Pk

i=1 xi /k,

x̄∗k =

b = (n − 2)−1
Σ

k(n − k)
b −1 (x̄k − x̄∗ ),
(x̄k − x̄∗k )T Σ
k
n

Pn

− k) and the pooled covariance estimator
!
k
n
X
X
T
∗
∗ T
(xi − x̄i )(xi − x̄i ) +
(xi − x̄i )(xi − x̄i )
.
i=k+1 xi /(n

i=1

The Hotelling’s

T2

i=k+1

test detects a change whenever max1≤k≤n max T 2 (k) exceeds a threshold.

The generalized likelihood ratio (GLR) statistic can be derived by assuming the null
and the alternative distributions are two multivariate normal distributions when both the
mean and the covariance matrix change, and the parameters all unknown. For a hypothetical change-point location k, the GLR statistic is given by
b n | − klog|Σ
b k | − (n − k)log|Σ
b ∗ |,
`(k) = nlog|Σ
k

P
k
∗
∗ T
T , and Σ
b ∗ = (n − k)−1 n
b k = k −1
where Σ
i=k+1 (xi − x̄i )(xi − x̄i ) .
i=1 (xi − x̄i )(xi − x̄i )
k
The GLR statistic detects a change whenever max1≤k≤n `(k) exceeds a threshold.
For our examples, we set n = Bmax = 200 for the Hotelling and the B-statistic, respectively. Let the change-point occurs at τ = 100, and choose the SL α = 0.05 (the probability
of incorrectly detecting a change under the null). Thresholds for the offline B-statistic are
obtained from Theorem 2 and for the other two methods are obtained from simulations.
Consider the following cases:
Case 1 (mean shift): observe a sequence of observations in R20 , whose distribution shifts
from N (0, I20 ) to N (0.1e, I20 );
Case 2 (mean shift with larger magnitude): observe a sequence of observations in R20 ,
whose distribution shifts from N (0, I20 ) to N (0.2e, I20 );
Case 3 (mean and local covariance change): observe a sequence of observations in R20 ,
whose distribution shifts from N (e, I20 ) to N (0.2e, Σ), where [Σ]11 = 2 and [Σ]ii = 1, i =
2, . . . , 20;
Case 4 (Gaussian to Laplace): observe a sequence of one-dimensional observations,
whose distribution shifts from N (0, 1) to Laplace distribution with zero mean and unit
variance. Note that the mean and the variance remain the same after the change.
We estimate the power for each case using 100 Monte Carlo trials. Table 3 shows that
the B-statistic achieves higher power than the Hotelling’s T 2 statistic and the GLR statistic
in all cases. The GLR statistic performs poorly, since when k small or closer to n, it needs to
estimate the pre-change and post-change sample covariance matrix using a limited number
of samples.
P

23

Li, Xie, Dai and Song

Table 3: Comparison of detection power for offline change-point detection. Thresholds for
all methods are calibrated so that the significance level is α = 0.05.

B-statistic
Hotelling’s T 2
GLR

Case 1
0.71
0.18
0.03

Case 2
1.00
0.88
0.05

Case 3
1.00
0.87
0.12

Case 4
0.44
0.03
0.04

6.2 Online change-point detection: Comparison with Hotelling’s T 2 statistic
Now consider online B-statistic with a fixed block-size B0 = 20. We compare the online
B-statistic with a Shewhart chart based on Hotelling’s T 2 statistic1 . At each time t, we
form a Hotelling’s T 2 statistic using the past B0 samples in [t − B0 + 1, t],
b −1 (x̄t − µ̂0 ),
T 2 (t) = B0 (x̄t − µ̂)T Σ
0
P
b 0 are estimated from reference data. The
where x̄t = ( ti=t−B0 +1 xi )/B0 , and µ̂0 and Σ
procedure detects a change-point whenever T 2 (t) exceeds a threshold for the first time.
The threshold for online B-statistic is obtained from Theorem 9, and from simulations for
the Hotelling’s T 2 statistic. To simulate EDD, let the change occur at the first point of the
testing data. Consider the following cases:
Case 1 (mean shift): distribution shifts from N (0, I20 ) to N (0.31, I20 );
Case 2 (covariance change): distribution shifts from N (0, I20 ) to N (0, Σ), where [Σ]ii =
2, where i = 1, 2, . . . , 5 and [Σ]ii = 1, i = 6, . . . , 20;
Case 3 (covariance change): distribution shifts from N (0, I20 ) to N (0, 2I20 );
Case 4 (Gaussian to Gaussian mixture): distribution shifts from N (0, I20 ) to mixture
Gaussian 0.3N (0, I20 ) + 0.7N (0, 0.1I20 );
Case 5 (Gaussian to Laplace)2 : distribution shifts from N (0, 1) to Laplace distribution
with zero mean and unit variance.
We evaluate the EDD for each case using 500 Monte Carlo trials. The results are summarized in Table 4. Note that in detecting changes in either Gaussian mean or covariance,
the online B-statistic performs competitively with Hotelling’s T 2 , which is tailored to the
Gaussian distribution. In the more challenging scenarios such as Case 4 and Case 5, the
Hotelling’s T 2 completely fails to detect the change-point whereas the online B-statistic can
detect the change fairly quickly.

7. Skewness correction
0 are normally disWe have shown that approximations to SL and ARL, assuming that ZB
0 does not converge to normal distribution
tributed, are reasonably accurate. However, ZB
even when B is large (see Appendix G) and it has a nonvanishing skewness, as illustrated

1. Here we do not compare the online B-statistic with the GLR statistic, since in our experiments Hotelling’s
T 2 consistently outperforms GLR when dimension is large.
2. For this difficult case, we report the EDD comparisons based on the selected 500 sequences where Bstatistic can successfully detect the changes which are defined as crossing the threshold within 50 steps
once the change occurs. Hotelling’s T 2 fails to detect the changes for all sequences.

24

Scan B-Statistic for Kernel Change-Point Detection

Table 4: Comparison of EDD in online change-point detection. The parameter is B0 = 20
and thresholds for all methods are calibrated so that ARL = 5000. Dashed line − means
that the procedure fails to detection the change, i.e., EDD is longer than 50.
Case 1
4.20
2.47

B-statistic
Hotelling’s T 2

Case 2
9.10
25.46

Case 3
1.00
1.27

Case 4
23.38
−

Case 5
23.03
−

by the following numerical example. Form 10000 instances of ZB computed using samples
from N (0, I20 ). Figs. 5(a)-(b) show the empirical distributions of ZB when N = 5, and
B = 2 or B = 200, respectively. Also plotted are the Gaussian probability density functions
with mean equal to the sample mean, and the variance predicted by Lemma 1. Note that
the empirical distributions of ZB match with Gaussian distributions to a certain extent but
the skewness becomes larger for larger B. Figs. 5(c)-(d) show the corresponding Q-Q plots.
400

400

B=2, N=5
300

300

200

200

100

100

0
−0.4

−0.2

0

0.2

B=200, N=5

0
−3

0.4

−2

−1

0

1

2

3
−3

x 10

(a): B = 2, N = 5, empirical distribution

(b): B = 200, N = 5, empirical distribution
−3

Quantiles of Input Sample

Quantiles of Input Sample

QQ Plot of Sample Data versus Standard Normal

0.4
0.2
0
−0.2
−0.4
−4

−2
0
2
Standard Normal Quantiles

4

(c): B = 2, N = 5, Q-Q plot

4

x 10

2
0
−2
−4
−4

−2
0
2
Standard Normal Quantiles

4

(d): B = 200, N = 5, Q-Q plot

Figure 5: Empirical distributions of ZB when B = 2 and B = 200, respectively. Note that
although Gaussian distribution seems to be a reasonable fit to the statistic, the skewness
becomes larger for larger B.

Incorporating the skewness of ZB , one can improve the accuracy of the approximations
for SL in Theorem 2 and for ARL in Theorem 9. Note that the log-moment-generating
function ψ(θ) defined in (11) corresponds to the so-called cumulant generating function
25

Li, Xie, Dai and Song

(McCullagh and Kolassa, 2009) and it has an expansion for θ close to zero:
ψ(θ) = κ1 θ +

κ2 2 κ3 3
θ + θ + o(θ3 ).
2
3!

0 ] = 0, the cumulants take values κ = E[Z 0 ] = 0, κ = Var[Z 0 ] = 1,
Since E[ZB
1
2
B
B
0 )3 ] − 3E[(Z 0 )2 ]E[Z 0 ] + 2(E[Z 0 ])3 = E[(Z 0 )3 ]. Recall that when deriving apκ3 = E[(ZB
B
B
B
B
proximations using change-of-measurement, we choose parameter θ such that the moment
0 is a standard normal, ψ(θ) = θ 2 /2, and hence θ = b.
generating function ψ̇(θ) = b. If ZB
Now with skewness correction, we approximate ψ(θ) as θ2 /2 + κ3 θ3 /6 when solving for θ.
Hence, we solve for
0 3 2
ψ̇(θ) ≈ θ + E[(ZB
) ]θ /2 = b,

and denote the solution to be θB (note that this time the solution may depend on B).
Moreover, with skewness correction, we will change the leading exponent term in (8) and
2
0
0
(32) from e−b /2 to be eψ(θB )−θB b . From numerical experiments, we find that the skewness
correction is especially useful when SL is small (e.g., α = 0.01) for the offline case or when
block size B0 is small (see Table 1, 2 and Fig. 4). The skewness correction can be important,
since it appears in the exponent of the expressions.
The remaining task is to estimate the skewness of B-statistic. Since ZB is zero-mean,
0 is related to the variance and third moment of Z via
the skewness of ZB
B
0 3
κ3 = E[(ZB
) ] = Var[ZB ]−3/2 E[Z3B ].

We already know how to estimate the variance of ZB from Lemma 1. The following lemma
3 ] in terms of the moments of the kernel h defined in
express the third-order moment E[ZB
(2):
Lemma 12 (Third-order moment of ZB )


1 
8(B − 2)
3
E h(x, x0 , y, y 0 )h(x0 , x00 , y 0 , y 00 )h(x00 , x, y 00 , y)
E[ZB ] = 2
2
2
B (B − 1)
N

3(N − 1) 
+
E h(x, x0 , y, y 0 )h(x0 , x00 , y 0 , y 00 )h(x000 , x0000 , y 00 , y)
2
N


(N − 1)(N − 2) 
0
0
00 000 0 00
0000 00000 00
+
E h(x, x , y, y )h(x , x , y , y )h(x , x , y , y)
N2


4
1 
+ 2
E h(x, x0 , y, y 0 )3
2
2
B (B − 1)
N


3(N − 1)
0
0 2
00 000
0
+
E
h(x,
x
,
y,
y
)
h(x
,
x
,
y,
y
)
N2


(N − 1)(N − 2) 
0
0
00 000
0
0000 00000
0
+
E
h(x,
x
,
y,
y
)h(x
,
x
,
y,
y
)h(x
,
x
,
y,
y
)
.
N2

(33)

Lemma 12 enables us to estimate the skewness efficiently, by reducing it to evaluating simpler terms in (33) that only requires estimating the statistic of the kernel function h(·, ·, ·, ·)
with tuples of samples. For instance, to evaluate E [h(x, x0 , y, y 0 )h(x0 , x00 , y 0 , y 00 )h(x00 , x, y 00 , y)],
we may each time draw a tuple of six samples without replacement from the reference data,
26

Scan B-Statistic for Kernel Change-Point Detection

treat them as x, x0 , x00 , y, y 0 and y 00 , compute the value of the function, repeat this many
times, and take average.
0 does not converge to Gaussian, the difference between its moment
Finally, although ZB
generating functions and that the standard normal distribution can be bounded, as we show
below. By applying a result on Page 220 of Yakir (2013), we have
0

E[eθZB ] − (1 +

|θ|3
θ2
0 3
0 2
) ≤ min{
E[|ZB
| ], θ2 E[|ZB
| ]}.
2
6

0 , we have
and if considering the skewness κ3 of ZB
0

E[eθZB ] − (1 +

θ4
θ2 θ3 κ3
0 4 1
0 3
+
) ≤ min{ E[|ZB
| ], |θ|3 E[|ZB
| ]}.
2
6
24
3

8. Real-data
We test the performance of scanning B-statistics for change-point detection on real data.
Our datasets include: (1) CENSREC-1-C: a real-world speech data set in the Speech Resource Consortium (SRC) corpora provided by National Institute of Informatics (NII)3 ; (2)
Human Activity Sensing Consortium (HASC) challenge 2011 data4 . We compare our modified B-statistic with a baseline algorithm, the relative density-ratio (RDR) estimate (Song
et al., 2013). One limitation of the RDR algorithm, however, is that it is not suitable for
high-dimensional data because estimating density ratio in the high-dimensional setting is
an ill-posed problem. To achieve reasonable performance for the RDR algorithm, we adjust
the bandwidth and the regularization parameter at each time step and, hence, the RDR
algorithm is computationally more expensive than using B-statistics. We adopt the Area
Under Curve (AUC) (Song et al., 2013) as a performance metric: the larger the AUC, the
better.
Our modified B-statistics have very competitive performance compared with the baseline RDR algorithm on the real data. Here we only report the main results and leave the
details in Appendix E. For speech data, our goal is to detect online the onset of a speech
signal emergent from the background. The backgrounds are taken from real acoustic signals,
such as noise in highway, airport and subway stations. The overall AUC for the B-statistic
is 0.8014 and for the baseline algorithm is 0.7578. For human activity detection data,
our goal is to detect a transition from one activity to another as quickly as possible. Each
instance consists six possible of human activity data collected by portable three-axis accelerometers. The overall AUC for the B-statistic is 0.8871 and for the baseline algorithm
is 0.7161.

9. Discussions
There are a few possible directions to extend our scan B-statistics. (1) Thus far, we have
assumed that data are i.i.d. from a null distribution P and when the change happens, data
are i.i.d. from an alternative distribution Q. Under these assumptions, we have developed
3. Available from http://research.nii.ac.jp/src/en/CENSREC-1-C.html
4. Available from http://hasc.jp/hc2011

27

Li, Xie, Dai and Song

the offline and online change-point detection algorithms based on the kernel two-sample
test statistic MMD. One may relax the temporal independence assumption and extend Bstatistics for dependent data by incorporating ideas from (Chwialkowski and Gretton, 2014).
(2) We have demonstrated how the number of blocks and block size affect the performance
of B-statistics. One can also explore how kernel bandwidth, as well as the dimensionality of
the data, would affect the performance. An empirical observation is that the performance
of MMD statistic degrades with the increasing dimensions in data. Some recent results for
kernel based test can be found in (Ramdas et al., 2015). We may adopt the idea of (Ramdas
et al., 2015) to extend our modified B-statistics for detecting a change in the dependence.
(3) For an exceedingly high dimensional data set with large Gram matrix, one can perform
random subsampling to reduce complexity similar to (Xie et al., 2015).

Acknowledgement
This research was supported in part by NSF CMMI-1538746, NSF CCF-1442635, NSF
CAREER CCF-1650913, a grant from Atlanta Police Foundation, gift donation from Adobe
Research to Yao Xie; NSF/NIH BIGDATA 1R01GM108341, ONR N00014-15-1-2340, NSF
IIS-1218749, NSF IIS-1639792, NSF CAREER IIS-1350983, grant from Intel and NVIDIA
to Le Song.

References
S. Arlot, A. Celisse, and Z. Harchaoui. Kernel change-point detection. arXiv preprint
arXiv:1202.3878, 6, 2012.
R. Arratia, L. Goldstein, and L. Gordon. Two moments suffice for Poisson approximations:
the Chen-Stein method. The Annals of Probability, pages 9–25, 1989.
M. Bibinger, M. Jirak, M. Vetter, et al. Nonparametric change-point analysis of volatility.
The Annals of Statistics, 45(4):1542–1578, 2017.
E. Brodsky and B. Darkhovsky. Nonparametric methods in change point problems, volume
243. Springer Science & Business Media, 2013.
K. Chwialkowski and A. Gretton. A kernel independence test for random processes. In
International Conference on Machine Learning (ICML), 2014.
M. Csörgő and L. Horváth. Invariance principles for changepoint problems. Journal of
Multivariate Analysis, 27(1):151–168, 1988.
M. Csörgö and L. Horváth. Limit theorems in change-point analysis, volume 18. John Wiley
& Sons Inc, 1997.
H. Dehling, R. Fried, I. Garcia, and M. Wendler. Change-point detection under dependence
based on two-sample u-statistics. In Asymptotic Laws and Methods in Stochastics, pages
195–220. Springer, 2015.
A. Dembo and O. Zeitouni. Large deviations techniques and applications. Springer, 2nd
edition, 2009.
28

Scan B-Statistic for Kernel Change-Point Detection

F. Desobry, M. Davy, and C. Doncarli. An online kernel change detection algorithm. IEEE
Transactions on Signal Processing, 53(8):2961–2974, 2005.
F. Enikeeva and Z. Harchaoui. High-dimensional change-point detection with sparse alternatives. arXiv:1312.1900, 2014.
G. Fasano and A. Franceschini. A multidimensional version of the kolmogorov–smirnov
test. Monthly Notices of the Royal Astronomical Society, 225(1):155–170, 1987.
L. Gordon and M. Pollak. An efficient sequential nonparametric scheme for detecting a
change of distribution. Annals of Statistics, 22(2):763–804, 1994.
A. Gretton, K.M. Borgwardt, M.J. Rasch, B. Schölkopf, and A. Smola. A kernel two-sample
test. The Journal of Machine Learning Research, 13(1):723–773, 2012.
Z. Harchaoui and O. Cappé. Retrospective mutiple change-point estimation with kernels.
In IEEE Workshop on Statistical Signal Processing (SSP), pages 768–772. IEEE, 2007.
Z. Harchaoui, F. Bach, and E. Moulines. Kernel change-point analysis. In Advances in
Neural Information Processing Systems (NIPS), 2008.
Z. Harchaoui, F. Bach, O. Cappe, and E. Moulines. Kernel-based methods for hypothesis
testing. IEEE Signal Processing Magazine, pages 87–97, 2013.
F. J . Massey Jr. The kolmogorov-smirnov test for goodness of fit. Journal of the American
statistical Association, 46(253):68–78, 1951.
D. Kifer, S.B. David, and J. Gehrke. Detecting change in data streams. In Proceedings of
the 30th International Conference on Very Large Data Bases-Volume 30, pages 180–191.
VLDB Endowment, 2004.
H.W. Lilliefors. On the kolmogorov-smirnov test for normality with mean and variance
unknown. Journal of the American statistical Association, 62(318):399–402, 1967.
D. Maragoni-Simonsen and Y. Xie. Sequential changepoint approach for online community
detection. IEEE Signal Processing Letters, 22(8):1035–1039, 2015.
D. Matteson and N. James. A nonparametric approach for multiple change point analysis
of multivariate data. Journal of the American Statistical Association, 109(505):334–345,
2014.
P. McCullagh and J. Kolassa. Cummulants. Scholarpedia, 4(3):4699, 2009.
D. Picard. Testing and estimating change-points in time series. Advances in applied probability, 17(04):841–867, 1985.
A. Ramdas, S.J. Reddi, B. Póczos, A. Singh, and L. Wasserman. On the decreasing power of
kernel and distance based nonparametric hypothesis tests in high dimensions. In TwentyNinth AAAI Conference on Artificial Intelligence (AAAI), 2015.
29

Li, Xie, Dai and Song

Z. E. Ross and Y. Ben-Zion. Automatic picking of direct P , S seismic phases and fault
zone head waves. Geophysical Journal International, 2014.
B. Schölkopf and A. Smola. Learning with kernels: support vector machines, regularization,
optimization, and beyond. MIT press, 2001.
B. Schölkopf and A. Smola. Learning with kernels. MIT Press, Cambridge, MA, 2002.
B. Schölkopf, K. Tsuda, and J.-P. Vert. Kernel methods in computational biology. MIT
Press, Cambridge, MA, 2004.
R. Serfling. Approximation theorems of mathematical statistics, volume 162. John Wiley &
Sons, 2009.
R. J. Serfling. U-Statistics. Approximation theorems of mathematical statistics. John Wiley
& Sons, 1980.
D. Siegmund. Sequential analysis: tests and confidence intervals. Springer, 1985.
D. Siegmund and E. S. Venkatraman. Using the generalized likelihood ratio statistic for
sequential detection of a change-point. Annals of Statistics, (23):255–271, 1995.
D. Siegmund and B. Yakir. Detecting the emergence of a signal in a noisy image. Statistics
and Its Inference, (1):3–12, 2008.
D. Siegmund, B. Yakir, and N. Zhang. Tail approximations for maxima of random fields by
likelihood ratio transformations. Sequential Analysis, 29(3):245–262, 2010.
L. Song, Y. Makoto, C. Nigel, and S. Masashi. Change-point detection in time-series data
by direct density-ratio estimation. Neural Networks, 43:72–83, 2013.
A. Tartakovsky, I. Nikiforov, and M. Basseville. Sequential analysis: Hypothesis testing and
changepoint detection. CRC Press, 2014.
T. Wan, K. Qu, Q. C. Zhang, R. A. Flynn, O. Manor, Z. Ouyang, J. Zhang, R. C. Spitale,
M.P. Snyder, E. Segal, et al. Landscape and variation of rna secondary structure across
the human transcriptome. Nature, 505(7485):706, 2014.
T. Wang, J. J. Wei, D.M. Sabatini, and E.S. Lander. Genetic screens in human cells using
the crispr-cas9 system. Science, 343(6166):80–84, 2014.
R. Wilcox. Kolmogorov–smirnov test. Encyclopedia of biostatistics, 2005.
B. Xie, Y. Liang, and L. Song. Scale up nonlinear component analysis with doubly stochastic
gradients. In Advances in neural information processing systems (NIPS), 2015.
Y. Xie and D. Siegmund. Sequential multi-sensor change-point detection. Annals of Statistics, 41(2):670–692, 2013.
B. Yakir. Multi-channel change-point detection statistic with applications in dna copynumber variation and sequential monitoring. In Proceedings of Second International
Workshop in Sequential Methodologies, pages 15–17, 2009.
30

Scan B-Statistic for Kernel Change-Point Detection

B. Yakir. Extremes in random fields: A theory and its applications. Wiley, 2013.
W. Zaremba, A. Gretton, and M. Blaschko. B-test: low variance kernel two-sample test.
In Advances in Neural Information Processing Systems (NIPS), 2013.
C. Zou, G. Yin, L. Feng, Z. Wang, et al. Nonparametric maximum likelihood approach to
multiple change-point problems. The Annals of Statistics, 42(3):970–1002, 2014a.
S. Zou, Y. Liang, H. V. Poor, and X. Shi. Nonparametric detection of anomalous data via
kernel mean embedding. arXiv:1405.2294, 2014b.

Appendix A. Recursive implementation of online B-statistic
The online B-statistic can be computed recursively via a simple update scheme. By its
construction, when time elapses from t to (t + 1), a new sample is added into the postchange block, and the oldest sample is moved to the reference pool. Each reference block is
updated similarly by adding one sample randomly drawn from the pool of reference data,
and the oldest sample is purged. Hence, only a limited number of entries in the Gram matrix
that are due to the new sample need to be updated. The update scheme is illustrated in Fig.
6 and explained in more details therein. The recursive scheme reduces the computational
cost of the online B-statistic to be linear in time. Similarly, the offline B-statistic can
also be computed recursively by utilizing the fact that ZB for B ∈ {2, . . . , Bmax } shares
many common terms. The recursive scheme reduces the computational cost of the offline
2 ).
B-statistic to be O(N Bmax

Appendix B. Variance and covariance calculation
(B)

Below, Xi,j , where i = 1, . . . , N , and j = 2, . . . , Bmax , denotes the j-th sample in the i-th
(B)

(B)

block Xi , and Yj denotes the j-th sample in Y (B) . The superscript B denotes the block
size. We start with proving Lemma 13 and Lemma 14, which are useful in proving Lemma
1.
Lemma 13 (Variance of MMD, under the null.) Under the null hypothesis,
h

2

Var MMD

(B)
(Xi , Y (B) )

i

 −1
B
=
E[h2 (x, x0 , y, y 0 )],
2

i = 1, . . . , N.

(34)

Proof For notational simplicity, below we drop the superscript B. Furthermore, we use x,
d
d
d
d
x0 , y and y 0 to denote generic samples, i.e., Xi,l = x, Xi,j = x0 , Yl = y, Yj = y 0 and they
d

are mutually independent of each other. Here the notation = means two random variables
have the same distribution. Below, we follow the same convention. For any i = 1, 2, . . . , n,
31

Li, Xie, Dai and Song

Time%t+1:%

Time%t:%
Xi(B0 ,t )

Xi(B0 ,t+1)

Y (B0 ,t )

Y (B0 ,t+1)

B

0
1
∑ k(Xi,(Bj 0 ,t ), Xi,l(B0 ,t ) )
B0 (B0 −1) j,l, j≠l

Xi(B0 ,t )

B

0
1
∑ k(Xi,(Bj 0 ,t ),Yl(B0 ,t ) )
B0 (B0 −1) j,l, j≠l

Xi(B0 ,t+1)

B

Y (B0 ,t )

0
1
∑ k(Yj(B0 ,t ),Yl(B0 ,t ) )
B0 (B0 −1) j,l, j≠l

Y (B0 ,t+1)

Figure 6: Recursive update scheme to compute the online B-statistics. The online Bstatistic is formed with N background blocks and one testing block and, hence, we keep
track of N Gram matrices. For illustration purposes, we partition the Gram matrix into
four windows (in red, black and blue, as shown on the left panel). At time t, to obtain
(B ,t)
MMD2 (Xi 0 , Y (B0 ,t) ), we compute the shaded elements and take an average within each
window. The diagonal entries in each window are removed to obtain an unbiased estimate.
(B ,t)
At time t + 1, we update Xi 0 and Y (B0 ,t) with the new data point and purge the oldest
data point, and update the Gram matrix by moving the colored window as shown on the
right panel, computing the elements within the new windows, and taking an average. Note
that we only need to compute the right-most column and the bottom row.

by definition of U-statistic, we have


 −1 X


B
Var MMD2 (Xi , Y ) = Var 
h(Xi,l , Xi,j , Yl , Yj )
2
l<j
 −2   



B
B
2 B−2
=
Var Ex,y [h(x, x0 , y, y 0 )]
2
2
1
2−1
  




B
2 B−2
+
Var h(x, x0 , y, y 0 ) .
2
2
2−2

(35)

Under null distribution, Ex,y [h(x, x0 , y, y 0 )] = 0. Thus, Var [Exi y [h(x, x0 , y, y 0 )]] = 0, and


Var h(x, x0 , y, y 0 ) = E[h2 (x, x0 , y, y 0 )] − E[h(x, x0 , y, y 0 )]2 = E[h2 (x, x0 , y, y 0 )].
Substitute these results into (35), we obtain the desired result (34).

32

Scan B-Statistic for Kernel Change-Point Detection

Lemma 14 (Covariance of MMD, under the null, same block size.) For s 6= 0, under null hypothesis
h
i
(B)
(B)
Cov MMD2 (Xi , Y (B) ), MMD2 (Xi+s , Y (B) )
 −1


B
=
Cov h(x, x0 , y, y 0 ), h(x00 , x000 , y, y 0 ) .
2

Proof For notational simplicity, we drop the superscript B. For i = 1, 2, . . . , N , and
s = (1 − i), (2 − i), . . . , (N − i), s 6= 0,


Cov MMD2 (Xi , Y ), MMD2 (Xi+s , Y )


 −1 X
 −1 X
B
B
= Cov 
h(Xi,l , Xi,j , Yl , Yj ),
h(Xi+s,p , Xi+s,q , Yp , Yq )
2
2
p<q
l<j
 −2   



B
B
2 B−2
=
Cov h(x, x0 , y, y 0 ), h(x00 , x000 , y, y 00 )
2
2
1
2−1

 −2   


2 B−2
B
B
Cov h(x, x0 , y, y 0 ), h(x00 , x000 , y, y 0 ) .
+
2−2
2
2
2
Under null distribution,


Cov h(x, x0 , y, y 0 ), h(x00 , x000 , y, y 00 )
Z
= h(x, x0 , y, y 0 )h(x00 , x000 , y, y 00 )dP(x, x0 , x00 , x000 , y, y 0 , y 00 )




Z Z
Z Z




0
0
0 0 
00 000
00
00 00 
000

= 
 h(x, x , y, y )dP(x , y ) dP(x) ·  h(x , x , y, y )dP(x , y ) dP(x ) = 0.
|
{z
}
|
{z
}
=0

=0

Above, with a slight abuse of notation, we use dP(·) to denote the probability measure of
appropriate arguments. Finally, we have:



Cov MMD2 (Xi , Y ), MMD2 (Xi+s , Y ) =

 −1


B
Cov h(x, x0 , y, y 0 ), h(x00 , x000 , y, y 0 ) .
2

33

Li, Xie, Dai and Song

B.1 Variance of B-statistics.
Proof [Proof for Lemma 1] For notational simplicity, we drop the superscript B. Using
results in Lemma 13 and Lemma 14, we have
#
N
1 X
2
Var[ZB ] = Var
MMD (Xi , Y )
N
i=1


X


1
Cov MMD2 (Xi , Y ; B), MMD2 (Xj , Y ) 
= 2 N Var[MMD2 (Xi , Y )] +
N
i6=j

 −1 


N −1
B
1
2
0
0
0
0
00 000
0
E[h (x, x , y, y )] +
Cov h(x, x , y, y ), h(x , x , y, y ) .
=
N
N
2
"

Next, we introduce Lemma 15 and Lemma 16, which are useful in proving Lemma 3.
Lemma 15 (Covariance of MMD, different block sizes, same block index.) For blocks
with the same index i but with distinct block sizes, under the null hypothesis we have
h
i B ∨ (B + v)−1
(B)
(B+v)
2
(B)
2
(B+v)
Cov MMD (Xi , Y ), MMD (Xi
,Y
) =
E[h2 (x, x0 , y, y 0 )].
2
(36)

Proof Note that
h
i
(B)
(B+v)
Cov MMD2 (Xi , Y (B) ), MMD2 (Xi
, Y (B+v) )


 −1 X

−1 B+v
B
X
B
B+v
= Cov 
h(Xi,l , Xi,j , Yl , Yj ),
h(Xi,p , Xi,q , Yp , Yq )
2
2
p<q
l<j


 −1 
−1
B
B+v
X
X
B
B+v
=
Cov 
h(Xi,l , Xi,j , Yl , Yj ),
h(Xi,p , Xi,q , Yp , Yq )
2
2
p<q
l<j
 −1 
−1 

B
B+v
B ∧ (B + v)
=
Var[h(x, x0 , y, y 0 )]
2
2
2


B ∨ (B + v) −1
=
E[h2 (x, x0 , y, y 0 )],
2
where the second last equality is due to a similar argument as before to drop block indices
as they are i.i.d. under the null.

34

Scan B-Statistic for Kernel Change-Point Detection

Lemma 16 (Covariance of MMD, different block sizes, different block indices.)
Under the null we have
h
i B ∨ (B + v)−1
(B)
(B+v)
2
(B)
2
(B+v)
Cov MMD (Xi , Y ), MMD (Xi+s , Y
) =
·
2


Cov h(x, x0 , y, y 0 ), h(x00 , x000 , y, y 0 ) .
Proof Note that

=

=

=
=

h
i
(B)
(B+v)
2
(B)
2
(B+v)
Cov MMD (Xi , Y ), MMD (Xi+s , Y
)


 −1 X

−1 B+v
B
X
B
B+v
(B)
(B)
(B)
(B)
(B+v)
(B+v)
Cov 
h(Xi,l , Xi,j , Yl , Yj ),
h(Xi+s,p , Xi+s,q , Yp(B+v) , Yq(B+v) )
2
2
p<q
l<j


 −1 
−1
B+v
B
X
X
B
B+v
(B)
(B)
(B)
(B)
(B+v)
(B+v)
h(Xi,l , Xi,j , Yl , Yj ),
h(Xi+s,p , Xi+s,q , Yp(B+v) , Yq(B+v) )
Cov 
2
2
p<q
l<j
 −1 
−1 



B
B+v
B ∧ (B + v)
Cov h(x, x0 , y, y 0 ), h(x00 , x000 , y, y 0 )
2
2
2

−1


B ∨ (B + v)
Cov h(x, x0 , y, y 0 ), h(x00 , x000 , y, y 0 ) ,
2

where the second last equality is due to a similar argument as before to drop block indices
as they are i.i.d. under the null.

B.2 Covariance of offline B-statistics.
Proof [Proof of Lemma 3] For the offline case, we have that the correlation
1
1
p
rB,B+v = p
Cov [ZB , ZB+v ] ,
Var[ZB ] Var[ZB+v ]
where

N
n
X
X
1
1
(B)
(B+v)
Cov (ZB , ZB+v ) = Cov 
MMD2 (Xi , Y (B) ),
MMD2 (Xj
, Y (B+v) )
N
N
i=1
j=1
h
i
1
(B)
(B+v)
= Cov MMD2 (Xi , Y (B) ), MMD2 (Xi
, Y (B+v) )
N
h
i
1 X
(B)
(B+v)
+ 2
Cov MMD2 (Xi , Y (B) ), MMD2 (Xj
, Y (B+v) ) .
N


i6=j

35

Li, Xie, Dai and Song

Using results from Lemma 15 and Lemma 16, we have:



1 B ∨ (B + v) −1
E[h2 (x, x0 , y, y 0 )]
Cov (ZB , ZB+v ) =
N
2




N − 1 B ∨ (B + v) −1
+
Cov h(x, x0 , y, y 0 ), h(x00 , x000 , y, y 0 )
N
2

−1 
B ∨ (B + v)
1
=
E[h2 (x, x0 , y, y 0 )]
N
2



N −1
0
0
00 000
0
+
Cov h(x, x , y, y ), h(x , x , y, y ) .
N

Finally, plugging in the expressions for Var[ZB ] and Var[ZB+v ], we have (10) for the offline
case.
(𝐵0 ,𝑡)

𝑀𝑀𝐷2 (𝑋𝑖
𝑋𝑖

, Y (B0,t) )

𝐵0 ,𝑡

𝑌

(𝐵0 ,𝑡+𝑠)

𝑀𝑀𝐷2 (𝑋𝑖
𝑋𝑖

𝐵0 ,𝑡+𝑠

(𝐵0 ,𝑡)

𝑀𝑀𝐷2 (𝑋𝑖

𝐵0 ,𝑡

𝑋𝑖

, Y (B0,t+s) )
𝑌

𝐵0 ,𝑡+𝑠

𝐵0 ,𝑡

(𝐵0 ,𝑡)

𝑀𝑀𝐷2 (𝑋𝑗
𝑋𝑗

, Y (B0,t) )
𝑌

𝐵0 ,𝑡

, Y (B0,t) )

𝐵0 ,𝑡+𝑠

𝑌

𝐵0 ,𝑡+𝑠

B0-s
B0-s

B0-s

B0-s

(a)

(b)

Figure 7: (a): Illustration of how MMD2 (Xi(B0 ,t) , Y (B0 ,t) ) and MMD2 (Xi(B0 ,t+s) , Y (B0 ,t+s) ) are constructed in the online change-point detection, where the shaded areas represent the overlapping data.
(B ,t)
(B ,t+s)
Illustration of how MMD2 (Xi 0 , Y (B0 ,t) ) and MMD2 (Xj 0
, Y (B0 ,t+s) ), j 6= i are constructed
in the online change-point detection, where the shaded areas represent the overlapping data. (b):
(B ,t)
(B ,t+s)
Illustration of how MMD2 (Xi 0 , Y (B0 ,t) ) and MMD2 (Xj 0
, Y (B0 ,t+s) ), j 6= i are constructed
in the online change-point detection, where the shaded areas represent the overlapping data.

B.3 Covariance of online B-statistics


0
0
Similarly, for the online case we need to analyze r0 = Cov ZB
,
Z
B0 ,t+s . Without loss of
0 ,t
generality, assume s > 0. We adopt the same strategy as the above for a fixed block size
36

Scan B-Statistic for Kernel Change-Point Detection

B0 to obtain


(B ,t)
(B ,t+s)
Cov MMD2 (Xi 0 , Y (B0 ,t) ), MMD2 (Xi 0
, Y (B0 ,t+s) )


 −1 X
 −1 X
B0
B0
B0
B0
(t+s)
(t+s)
(t)
(t)
(t)
(t)
h(Xi,p , Xi,q , Yp(t+s) , Yq(t+s) )
h(Xi,l , Xi,j , Yl , Yj ),
= Cov 
2
2
p<q
l<j

 −2 
(B0 − s) ∨ 0
B0
Var[h(x, x0 , y, y 0 )].
(37)
=
2
2
(B ,t)

(B ,t+s)

Figure 7 demonstrates how MMD2 (Xi 0 , Y (B0 ,t) ) and MMD2 (Xi 0
, Y (B0 ,t+s) ) are
constructed. The shaded areas represent the overlapping data.
Similarly, we have


(B ,t)
(B ,t+s)
Cov MMD2 (Xi 0 , Y (B0 ,t) ), MMD2 (Xj 0
, Y (B0 ,t+s) )
"  B
#
 −1 X
B0
0
B0 −1 X
B0
(t)
(t)
(t)
(t)
(t+s)
(t+s)
= Cov
h(Xj,p , Xj,q , Yp(t+s) , Yq(t+s) )
h(Xi,l , Xi,k , Yl , Yk ),
2
2
p<q
l<k
 −2 

B0
(B0 − s) ∨ 0
=
Cov(h(x, x0 , y, y 0 ), h(x00 , x000 , y, y 0 )),
(38)
2
2
(B ,t)

(B ,t+s)

Figure 7 demonstrates how MMD2 (Xi 0 , Y (B0 ,t) ) and MMD2 (Xj 0
, Y (B0 ,t+s) ), j 6= i
are constructed. The shaded areas represent the overlapping data. Thus,
Cov (ZB0 ,t , ZB0 ,k+s )


N
N
X
X
1
1
(B ,t)
(B ,t+s)
= Cov 
MMD2 (Xi 0 , Y (B0 ,t) ),
MMD2 (Xj 0
, Y (B0 ,t+s) )
N
N
i=1
j=1
h
 −2 
1
(B0 − s) ∨ 0
B0
=
Var(h(x, x0 , y, y 0 ))
2
2
N
i
N −1
+
Cov(h(x, x0 , y, y 0 ), h(x00 , x000 , y, y 0 )) .
N
Finally, plugging in the expressions for Var[ZB0 ,t ] and Var[ZB0 ,t+s ], we have (30) for the
online case.

Appendix C. Proof of Lemma 10
Proof First, we show that the stopping time T is asymptotically exponentially distributed.
The analysis of the distribution of the stopping time is based on Poisson approximation.
0
Define an indicator of the event 1j such that the event 1{max(j−1)m≤t≤jm ZB
> b}.
0 ,t
Consider the time interval [0, x]. Note that the stopping time is not activated in the interval
[0, x], if and only if, all the relevant indicators are zero. For simplicity, we assume x is
c = Px/m 1j . Hence, {W
c = 0} = {Tb > x}.
divisible by m. Define the random variable W
j=1
37

Li, Xie, Dai and Song

Thus, to final the tail probability of the stopping time, we need to show that the sum of
the indicator functions converge to a Poisson distribution.
For online B-statistics, the standard Poisson limit cannot be directly applied, since the
events {1j }, j = 1, . . . , x/m, are not independent, and we need the generalized Poisson
limit theorem (Arratia et al., 1989), which allows for dependence between the variables.
The setup for the theorem is as follows. Let I be an arbitrary index set, andPfor α ∈ I, let
Xα be a Bernoulli random variable with pα = P(Xα = 1) > 0. Let W = α∈I Xα . For
each α ∈ I, suppose we choose Bα ⊂ I with α ∈ Bα . Think of Bα as a “neighborhood
of dependence” for each α, such P
that XP
α is independent or
Pnearly
P independent of all of
the
X
for
β
∈
/
B
.
Define
p
=
p
p
,
p
=
α
1
2
β
α∈I
β∈Bα α β
αI
α6=β∈Bα E(Xα Xβ ), p3 =
P
α∈I E|E(Xα − pα |σ(Xβ : β ∈ I − Bα ))|, where σ(·) represents the σ-field generated by
the corresponding random field. Loosely speaking, p1 measures the neighborhood size, p2
measures the expected number of neighbors of a given occurrence and p3 measures the
dependence between an event and the number of occurrences outside its neighborhood.
Then, we have the following
Theorem 17 (Poisson approximation, Theorem 1 in (Arratia et al., 1989)) Let W
be the number of occurrences of dependent events, and let Z be a Poisson random variable
with EZ = EW = λ > 0. Then the total variation distance between the distributions of W
and Z is bounded by
sup |Eh(W ) − Eh(Z)| ≤ p1 + p2 + p3 .
khk=1

where h :

Z+

→ R, khk = supk≥0 |h(k)|.

The theorem is a consequence of the powerful Chen-Stein method.
Invoking the above theorem in our online B-statistics setting, we can bound the total
variation distance between the random variable, defined as the number of boundary cross
events for the statistic over disjoint intervals, and a Poisson random variable with the
same rate. In our setting, let I = {1, 2, . . . , x/m} and N (j) = {j − 1, j, j + 1} where
j = 2, . . . (x/m − 1) (with obvious modifications for j = 1 and j = x/m). Then we can
specify:
X X
p1 =
P{1j = 1}P{1i = 1} = 2(x/m − 2)P{11 = 1}2 + 2P{11 = 1},
(39)
j∈I i∈N (j)\{j}

p2 =

X

X

P{1j = 1, 1i = 1} = 2(x/m − 1)P{11 = 1, 12 = 1},

(40)

j∈I i∈N (j)\{j}

p3 =

X

E {|E{1j |σ{1i : i 6∈ N (j)}} − E{1j }|} .

(41)

j∈I

We will show that p1 , p2 , and p3 converge to 0 as b → ∞. For p1 , the last summand
in (39) is associated with the two edge elements. It follows that p1 is asymptotically to
(2C + 2)P{11 = 1}, which will converge to zero as b → ∞ since P{11 = 1} converges to
1 2
zero when m is sub-exponential, i.e., log b  m  b−1 e 2 b . Next, let us examine p2 in (40).
Redefine parameter sub-region
S1 = [0, m − B0 /2],

S2 = [m − B0 /2, m + B0 /2],
38

S3 = [m + B0 /2, 2m],

Scan B-Statistic for Kernel Change-Point Detection

0
and denote Yi , i = 1, 2, 3 as {Yi = 1} = {maxt∈Si ZB
> b}, which are the indicator func0 ,t
tions of crossings of the threshold in the approximate sub-regions. Notice that the indicator
functions Y1 and Y3 are independent of each other and they share the same distribution. We
use the fact that unless the crossing occurs in a shared sub-region, it must simultaneously
occur in two disjoint sub-regions in order to have double crossing. As a consequence, we
obtain the upper bound 11 · 12 ≤ Y2 + Y1 · Y3 , and

P{11 = 1, 12 = 1} ≤ P{Y2 = 1} + P{Y1 = 1}2 ≤ P{Y2 = 1} + P{11 = 1}2 .
1 2

The probability P{Y2 = 1} is proportional to B0 ·be− 2 b . Consequently, p2 is asymptotically
1 2
bounded by 2C(B0 /m + P{11 = 1}). Hence, p2 converges to zero if log b  m  b−1 e 2 b
whenever b → ∞. For p3 in (41), 1j and 1i are computed over non-overlapping observations
and are therefore independent. Thus, the term p3 vanishes.
Next prove that the collection of stopping times {Tb } indexed by b is uniformly integrable.
Again consider the sequence of indicators {1j }, j = 2k and k = 1, 2, . . . . Define the random
variable τ that identifies the index of the first indicator in the sequence that obtains the
value one: τ = inf{k : 12k = 1}. Note that τ has a geometric distribution. Moreover, since
Tb ≤ 2mτ we obtain that
P{Tb > x} ≤ P{τ > x/(2m)} = (1 − P(12 = 1))bx/(2m)c .
The conclusion then follows from that 1/m · P(12 = 1) converges to 0.

Appendix D. Dependence of EDD on block-size B0
We investigate the dependence of EDD of the B-statistic for online change-point detection,
on the pre-defined block size B0 . This example may also shed some light on how to choose
B0 in practice. The rationale is that, on the one-hand, the expected detection delay (EDD)
is greater than B0 , since we need at least B0 samples to compute one B-statistic. Hence,
a large B0 will artificially impose a longer EDD. On the other hand, when block size is
too small, we may not pool enough post-change samples and the statistical power of the
B-statistic is weak, which will also result in a large EDD. Hence, there should be an optimal
choice for B0 that minimizes the EDD. This is validated by our numerical example. Consider
the distribution shifts from N (0, I20 ) to N (µ, I20 ), with µ being element-wise equal to a
non-zero constant. In Figure 8(a), the shift size of the mean is 0.2, where the minimum
EDD is achieved by B0 = 28. Figure 8(b) shows the optimal block sizes for a range of values
of the mean shift size (from 0.1 to 1.2).
We fix B0 = 20 and let x ∈ R20 . Consider a simple case where signal changes from a zeromean Gaussian N (0, I) to a non-zero mean Gaussian N (µ, I20 ), and the post-change mean
vector µ is element-wise equal to a constant. Note when we implement this experiment, we
let the mean shift change occur at the first point of the testing data. Under this scenario,
the EDD would be at least B0 , since we need at least B0 testing samples to compute one
M -statistic. In this way, we intentional get the worst case for EDD, which serve as an upper
bound for the actually EDD(s) in more general sequential settings.
39

Li, Xie, Dai and Song

8

100

7

80

Optimal Block Size

log(EDD)

In Figure 8(a), we fix the Gaussian mean shift µ to be 0.2. Figure 8(a) shows the
log (EDD) versus B0 and demonstrates the tradeoff in choosing the block size for the M statistic: when block size is too small, the statistical power of the M -statistic is weak and,
hence, the EDD is large; on the other hand, when block size is too large, although statistical
power is good, the EDD is also large because we have to pool a large number of samples for
the post-change block which incurs a larger delay. In view of the above, there is an optimal
choice for the block size, and here B0 = 28 yields the smallest (EDD), which is defined as
the optimal block size for a mean shift equals 0.2. Similarly, we evaluate the optimal block
sizes for a range of values of the mean shift (from 0.1 to 1.2), as illustrated in Figure 8(b).
In doing this, we hope to provide a guidance in how to set B0 in practice.

6

5
Optimal B0=28
4

3
0

60

40

20

20

40
60
Block Size

80

0

100

0.2

0.4

0.6
0.8
Mean Shift Size

(a)

1

1.2

(b)

Figure 8: Online B-statistic for a case where the distribution shifts from N (0, I20 ) to
N (µ, I20 ), with µ being element-wise equal to a non-zero constant. (a) log(EDD) versus
block size and the optimal block size B0 = 28 corresponds to a minimum EDD; (b) optimal
block sizes versus the size of the mean shift.

Appendix E. More details for real-data experiments
E.1 CENSREC-1-C Speech dataset
CENSREC-1-C is a real-world speech dataset in the Speech Resource Consortium (SRC)
corpora provided by National Institute of Informatics (NII)5 . This dataset contains two
categories of data: (1) Simulated data. The simulated speech data are constructed by
concatenating several utterances spoken by one speaker. Each concatenated sequence is
then added with 7 different levels of noise from 8 different environments. So there are totally
56 different noise. Each noise setting contains 104 sequences from 52 males and 52 females
speakers. (2) Recording data. The recording data is from two real-noisy environments
(in university restaurant and in the vicinity of highway), and with two Signal Noise Ratio
(SNR) settings (lower and higher). Ten subjects were employed for recording, and each one
has four speech sequence data.
Experiment Settings. We will compare our algorithm with the baseline algorithm from
(Song et al., 2013). (Song et al., 2013) only utilized 10 sequences from “STREET SNR HIGH”
5. Available from http://research.nii.ac.jp/src/en/CENSREC-1-C.html

40

Scan B-Statistic for Kernel Change-Point Detection

Table 5: AUC results in CENSREC-1-C speech dataset. Simulated data are from 8
noise categories, and with two different noise levels (clean(C) and SNR 20db (S)); Recording data are from RESTAURANT SNR HIGH (RH), RESTAURANT SNR LOW (RL),
STREET SNR HIGH (SH) and STREET SNR LOW (SL).
(a) Recording data

Ours
Baseline

RH
0.7800
0.7503

RL
0.7282
0.6835

SH
0.6507
0.4329

SL
0.6865
0.6432

(b) Simulate clean data

Ours
Baseline

C1
0.9413
0.9138

C2
0.9446
0.9262

C3
0.9236
0.8691

C4
0.9251
0.9128

C5
0.9413
0.9138

C6
0.9446
0.9216

C7
0.9236
0.8691

C8
0.9251
0.9128

S7
0.6796
0.6487

S8
0.7145
0.6541

(c) Simulated data with SNR=20db

Ours
Baseline

S1
0.7048
0.7083

S2
0.7160
0.6681

S3
0.7126
0.6490

S4
0.7129
0.7119

S5
0.7094
0.6994

S6
0.7633
0.6815

setting in recording data. Here we will use all the settings in recording data, the SNR level
20db and clean signals from simulated data. See Figure 9 for some examples of the testing
data, as well as the statistics computed by our algorithm. For each sequence, we decompose
it into several segments. Each segment consists of two types of signals (noise vs speech).
Given the reference data from noise, we want to detect the point where the signal changes
from noise to speech.
Evaluation Metrics. We use Area Under Curve (AUC) to evaluate the computed statistics, like in (Song et al., 2013). Specifically, for each test sequence that consists of two signal
distributions, we will mark the points as change-points whose statistics exceed the given
threshold. If the distance between detected point and true change-point is within the size
of detection window, then we consider it as True Alarm (True Positive). Otherwise it is a
False Alarm (False Positive).
We use 10% of the sequences to tune the parameters of both algorithms, and use the rest
90% for reporting AUC. The kernel bandwidth is tuned in {0.1dmed , 0.5dmed , dmed , 2dmed , 5dmed },
where dmed is the median of pairwise distances of reference data. Block size is fixed to 50,
and the number of blocks is simply tuned in {10, 20, 30}.
Results. Table 5 shows the AUC of two algorithms on different background settings. Our
algorithm surpasses the baseline on most cases. Both algorithms are performing quite well
on the simulated clean data, since the difference between speech signals and background
is more significant than the noisy ones. The averaged AUC of our algorithm on all these
settings is .8014, compared to .7578 achieved by baseline algorithm. See the ROC curves
in Figure 10 for a better comparison.
41

Li, Xie, Dai and Song

Original Signal

2000

Original Signal

2000

Original Signal

5000

0
0

0

-2000
-4000

400

500

600

700

800

900

Online Statistics

10

-2000

200

400

5

5 b=3.56

0

600

800

-5000

1000

Online Statistics

400

500

600

700

800

900

800

900

800

900

800

900

Online Statistics

10
b=3.56

b=3.56
0

0

-5

-5

-10

400

500

600

700

800

900

(a) RESTAURANT SNR HIGH
2

#10

4

Original Signal

-2

600

800

-10

1000

#10

4

Original Signal

400

500

600

700

800

900

Online Statistics

-1

0

500

600

700

(c) Clean-1
Original Signal

0

400

600

800

1000

Online Statistics

5

b=3.56

400

5000

0

5
0

400

(b) STREET SNR HIGH
1

0

200

-5000

400

b=3.56

500

600

700

Online Statistics

5
b=3.56
0

-5
-10

-5
400

500

600

700

800

900

(d) RESTAURANT SNR LOW

-10

400

600

800

1000

(e) STREET SNR LOW

-5

400

500

600

700

(f) SNR 20db-1

Figure 9: Examples of speech dataset. The red vertical bar shown in the upper part of each
figure is the ground truth of change-point; The green vertical bar shown in the lower part
is the change-point detected by our algorithm (the point where the statistic exceeds the
threshold). We also plot the threshold as a red dash horizontal line in each figure. Once
the statistics touch the threshold, we will stop the detection.
E.2 HASC human activity dataset
This data is from Human Activity Sensing Consortium (HASC) challenge 20116 . Each
data consists of human activity information collected by portable three-axis accelerometers.
Following the setting in (Song et al., 2013), we use the `2 -norm of 3-dimensional data (i.e.,
the magnitude of acceleration) as the signals.
We use the ‘RealWorldData’ from HASC Challenge 2011, which consists of 6 kinds of
human activities:
walk/jog, stairUp/stairDown, elevatorUp/elevatorDown,
escalatorUp/escalatorDown, movingWalkway, stay.
We make pairs of signal sequences from different activity categories, and remove the sequences which are too short. We finally get 381 sequences. We tune the parameters using
the same way as in CENSREC-1-C experiment. The AUC of our algorithm is .8871, compared to .7161 achieved by baseline algorithm, which greatly improved the performance.
Examples of the signals are shown in Figure 11. Some sequences are easy to find the
change-point, like Figure 11a, and 11d. Some pairs of the signals are hard to distinguish visually, like Figure 11b and 11c. The examples show that our algorithm can tell
6. http://hasc.jp/hc2011

42

Scan B-Statistic for Kernel Change-Point Detection

0

True Positive Rate

True Positive Rate

0.4

0.2

0

0.5

0.8

0.2

0

True Positive Rate

True Positive Rate

0.4

0.2

0.5

0.2

1

0

simulate/N1_SNR20

True Positive Rate

True Positive Rate

0.4

0.2

0.5

simulate/N5_SNR20

True Positive Rate

True Positive Rate

0.4

0.2

0

0.5

False Positive Rate

(q) S5

1

0.2

0.5

0

ours
baseline

0.8

0.6

0.4

0.2

0

1

0

1

False Positive Rate

0.4

0.2

0.5

1

False Positive Rate

(r) S6

(s) S7

Figure 10: AUC comparison on speech dataset
43

simulate/N8_SNR20

1

ours
baseline

0

1

(p) S4

0.6

0

0.5

False Positive Rate

simulate/N7_SNR20

0.8

1

simulate/N4_SNR20

1

0.2

0.5

0.5

(l) C8

0.4

1

0.4

0

0.2

(o) S3

0.6

0

0.4

False Positive Rate

ours
baseline

0.8

0.6

False Positive Rate

0.6

0

1

ours
baseline

0.8

0

1

ours
baseline

0.8

0

1

simulate/N6_SNR20

1

ours
baseline

0.5

simulate/N3_SNR20

(n) S2

0.6

0

0

0.5

simulate/clean8

1

0.2

False Positive Rate

(m) S1

0.8

0.5

0

(h) C4

0.4

1

0.2

0

0.2

(k) C7

0.4

0

1

0.4

False Positive Rate

0.6

False Positive Rate

1

1

ours
baseline

0.8

0.6

False Positive Rate

0.6

0

simulate/N2_SNR20

1

ours
baseline

0

0.8

1

ours
baseline

0.8

0

1

ours
baseline

(j) C6

0.6

0

0.5

0.5

simulate/clean7

False Positive Rate

(i) C5

0.8

0.2

1

0.4

False Positive Rate

1

0.4

0

0.5

simulate/clean4

1

(g) C3

0.6

0

0

False Positive Rate

ours
baseline

0.8

0.2

(d) SL

0.6

0

1

simulate/clean6

1

ours
baseline

0

0.8

(f) C2

simulate/clean5

0.4

False Positive Rate

ours
baseline

False Positive Rate

0.6

0

0.5

0.6

0

1

simulate/clean3

1

0.4

0

1

0.5

ours
baseline

0.8

(c) SH

0.6

(e) C1

0.8

0

recording/STREET_SNR_LOW

False Positive Rate

ours
baseline

False Positive Rate

1

1

simulate/clean2

1

ours
baseline

0.6

0

0.2

(b) RL

simulate/clean1

0.8

0.4

False Positive Rate

(a) RH
1

0.5

0.6

True Positive Rate

1

False Positive Rate

0.8

0

True Positive Rate

0.5

0

True Positive Rate

0

0.2

True Positive Rate

0

0.4

1

ours
baseline

True Positive Rate

0.2

0.6

recording/STREET_SNR_HIGH

True Positive Rate

0.4

0.8

True Positive Rate

0.6

1

ours
baseline

True Positive Rate

0.8

1

True Positive Rate

True Positive Rate

ours
baseline

True Positive Rate

recording/RESTAURANT_SNR_LOW

1

True Positive Rate

recording/RESTAURANT_SNR_HIGH

ours
baseline

0.8

0.6

0.4

0.2

0

0

0.5

False Positive Rate

(t) S8

1

Li, Xie, Dai and Song

the change-point from walk to stairUp/stairDown, or from stairUp/stairDown to escalatorUp/escalatorDown. There are some cases when our algorithm raises false alarm. See
Figure 11h. It find a change-point during the activity ‘elevatorUp/elevatorDown’. It is reasonable, since this type of action contains the phase from acceleration to uniform motion,
and the phase from uniform motion to acceleration.
Original Signal

2
1

2

0

0

500
600
700
Online Statistics

20
0

800

0

500

600

700

800

500
600
700
Online Statistics

800

-20

0

500

600

700

800

1

5000

-20

800

0

500

800

600

700

Original Signal

1.1

2

500
600
700
Online Statistics

800

(c) A1 vs A2

Original Signal

4

500
600
700
Online Statistics
b=3.38

(b) A1 vs A4

Original Signal

0

0
20

b=3.38

(a) A1 vs A6
2

Original Signal

4
2

20

b=3.38

-20

Original Signal

4

1

400

10

500
600
Online Statistics

700

0.9
500

500
600
700
Online Statistics

800

500

800

b=3.38
0 b=3.38
-5000

0 b=3.38

0

500

600

700

800

-10

400

(d) A6 vs A2

700

2000

500
600
Online Statistics

700

0

600

700

Original Signal

4

1

400

-500

(f) A4 vs A3

Original Signal

2

2
0

600

(e) A2 vs A4

Original Signal

4

500

2

200

5000

400
600
800
Online Statistics

1000

0

400

10

500
600
Online Statistics

700

b=3.38
0 b=3.38
-2000

400

0 b=3.38

500

600

(g) A6 vs A4

700

-5000

200

0

400

600

800

(h) A3 vs A1

1000

-10

400

500

600

700

(i) A2 vs A6

Figure 11: Examples of HASC dataset. The markers in this figure are the same as in
Figure 9.

Appendix F. Skewness of B-statistics
In the following, Lemma 18, Lemma 19, and Lemma 20 are used to derive the final expression
for the skewness of the statistic:

44

Scan B-Statistic for Kernel Change-Point Detection

Lemma 18 Under null hypothesis,
h
3 i
E MMD2 (Xi , Y )
=





8(B − 2)
4
E h(x, x0 , y, y 0 )h(x0 , x00 , y 0 , y 00 )h(x00 , x, y 00 , y) + 2
E h(x, x0 , y, y 0 )3 .
2
2
2
B (B − 1)
B (B − 1)

Proof
h
E

MMD2 (Xi , Y )

3 i


!3 
 −3
X
B
=
E
h(Xi,a , Xi,b , Ya , Yb ) 
2
a<b
 −3 X
B
=
Ck E [hab hcd hef ] ,
2
k

where for simplicity we write hab = h(Xi,a , Xi,b , Ya , Yb ) and define Ck the corresponding
number of combination under specific structure. Most of the terms in E [hab hcd hef ] vanish under the null. By enumerating all the combinations, only two terms are nonzero:
E [hab hbc hca ] and E [hab hab hab ]. Then,
 −3  
 −3  
h
3 i
B
B
B
B
E MMD2 (Xi , Y )
=
2(B − 2)E [hab hbc hca ] +
E [hab hab hab ]
2
2
2
2
8(B − 2)
E [h(Xi,a , Xi,b , Ya , Yb )h(Xi,b , Xi,c , Yb , Yc )h(Xi,c , Xi,a , Yc , Ya )]
= 2
B (B − 1)2


4
+ 2
E h(Xi,a , Xi,b , Ya , Yb )3
2
B (B − 1)




8(B − 2)
4
0
0
0 00 0 00
00
00
0
0 3
= 2
E
h(x,
x
,
y,
y
)h(x
,
x
,
y
,
y
)h(x
,
x,
y
,
y)
+
E
h(x,
x
,
y,
y
)
.
B (B − 1)2
B 2 (B − 1)2

Lemma 19 Under null hypothesis,
h
i
2
E MMD2 (Xi , Y ) MMD2 (Xj , Y )

i6=j



8(B − 2)
E h(x, x0 , y, y 0 )h(x0 , x00 , y 0 , y 00 )h(x000 , x0000 , y 00 , y)
= 2
2
B (B − 1)


4
+ 2
E h(x, x0 , y, y 0 )2 h(x00 , x000 , y, y 0 ) .
2
B (B − 1)
Proof
h

i
2
MMD2 (Xi , Y ) MMD2 (Xj , Y )
i6=j

!2
!
 −3
X
X
B
=
E
h(Xi,a , Xi,b , Ya , Yb )
h(Xj,a , Xj,b , Ya , Yb ) 
2
a<b
a<b
 −3 X
B
Ck E [hi,ab hi,cd hj,ef ] ,
=
2
E

k

45

Li, Xie, Dai and Song

where for simplicity we write hi,ab = h(Xi,a , Xi,b , Ya , Yb ) and define Ck the corresponding
number of combination under specific structure. Similarly, most of the terms in E [hi,ab hi,cd hj,ef ]
vanish under the null. By enumerating all the combinations, only two terms are nonzero:
E [hi,ab hi,bc hj,ca ] and E [hi,ab hi,ab hj,ab ]. Then,
h
i
2
E MMD2 (Xi , Y ) MMD2 (Xj , Y )
i6=j
 −3  
 −3  
B
B
B
B
=
2(B − 2)E [hi,ab hi,bc hj,ca ] +
E [hi,ab hi,ab hj,ab ]
2
2
2
2
8(B − 2)
E [h(Xi,a , Xi,b , Ya , Yb )h(Xi,b , Xi,c , Yb , Yc )h(Xj,c , Xj,a , Yc , Ya )]
= 2
B (B − 1)2


4
+ 2
E h(Xi,a , Xi,b , Ya , Yb )2 h(Xj,a , Xj,b , Ya , Yb )
2
B (B − 1)


8(B − 2)
0
0
0 00 0 00
000 0000 00
E
h(x,
x
,
y,
y
)h(x
,
x
,
y
,
y
)h(x
,
x
,
y
,
y)
= 2
B (B − 1)2


4
0
0 2
00 000
0
+ 2
E
h(x,
x
,
y,
y
)
h(x
,
x
,
y,
y
)
.
B (B − 1)2

Lemma 20 Under null hypothesis,


E MMD2 (Xi , Y )MMD2 (Xj , Y )MMD2 (Xr , Y ) i6=j6=r
=



8(B − 2)
E h(x, x0 , y, y 0 )h(x00 , x000 , y 0 , y 00 )h(x0000 , x00000 , y 00 , y)
2
2
B (B − 1)


4
0
0
00 000
0
0000 00000
0
+ 2
E
h(x,
x
,
y,
y
)h(x
,
x
,
y,
y
)h(x
,
x
,
y,
y
)
.
B (B − 1)2

Proof Note that


E MMD2 (Xi , Y )MMD2 (Xj , Y )MMD2 (Xr , Y ) i6=j6=r


!
!
 −3
X
X
X
B
=
E
h(Xi,a , Xi,b , Ya , Yb )
h(Xj,c , Xj,d , Yc , Yd ) 
h(Xr,e , Xr,f , Ye , Yf )
2
a<b
c<d
e<f
 −3 X
B
=
Ck E [hi,ab hj,cd hr,ef ] .
2
k

Similarly, most of the terms in E [hi,ab hj,cd hr,ef ] vanish under the null. By enumerating
all the combinations, only two terms are nonzero: E [hi,ab hj,bc hr,ca ] and E [hi,ab hj,ab hr,ab ].
46

Scan B-Statistic for Kernel Change-Point Detection

Then,


E MMD2 (Xi , Y )MMD2 (Xj , Y )MMD2 (Xr , Y ) i6=j6=r
 −3  
 −3  
B
B
B
B
=
2(B − 2)E [hi,ab hj,bc hr,ca ] +
E [hi,ab hj,ab hr,ab ]
2
2
2
2
8(B − 2)
= 2
E [h(Xi,a , Xi,b , Ya , Yb )h(Xj,b , Xj,c , Yb , Yc )h(Xr,c , Xr,a , Yc , Ya )]
B (B − 1)2
4
E [h(Xi,a , Xi,b , Ya , Yb )h(Xj,a , Xj,b , Ya , Yb )h(Xr,a , Xr,b , Ya , Yb )]
+ 2
B (B − 1)2


8(B − 2)
= 2
E h(x, x0 , y, y 0 )h(x00 , x000 , y 0 , y 00 )h(x0000 , x00000 , y 00 , y)
2
B (B − 1)


4
+ 2
E h(x, x0 , y, y 0 )h(x00 , x000 , y, y 0 )h(x0000 , x00000 , y, y 0 ) .
2
B (B − 1)

Proof We can write the raw third-order moment as

!3 
N
X
1
3
MMD2 (Xi , Y ) 
E[ZB
] =E 
N
i=1


! N
!
N
N
X
X
1  X
= 3E
MMD2 (Xi , Y ) 
MMD2 (Xj , Y )
MMD2 (Xr , Y ) 
N
r=1
i=1
j=1
  
 h
h
i
i

2
1
N −1
1 3 N
3
= 3 N E MMD2 (Xi , Y )
E MMD2 (Xi , Y ) MMD2 (Xj , Y )
+ 3
N
N 2
1
1
i6=j
 




1 N
N −1 N −2
+ 3
E MMD2 (Xi , Y )MMD2 (Xj , Y )MMD2 (Xr , Y ) i6=j6=r
N
1
1
1
h
i
i

2
1
3(N − 1) h
3
2
2
= 2 E MMD2 (Xi , Y )
+
E
MMD
(X
,
Y
)
MMD
(X
,
Y
)
i
j
N
N2
i6=j


(N − 1)(N − 2)
+
E MMD2 (Xi , Y )MMD2 (Xj , Y )MMD2 (Xr , Y ) i6=j6=r
2
N






1
8(B − 2)
4
0
0
0 00 0 00
00
00
0
0 3
= 2
E
h(x,
x
,
y,
y
)h(x
,
x
,
y
,
y
)h(x
,
x,
y
,
y)
+
E
h(x,
x
,
y,
y
)
N
B 2 (B − 1)2
B 2 (B − 1)2



3(N − 1)
8(B − 2)
E h(x, x0 , y, y 0 )h(x0 , x00 , y 0 , y 00 )h(x000 , x0000 , y 00 , y)
+
2
2
2
N
B (B − 1)



4
0
0 2
00 000
0
+ 2
E h(x, x , y, y ) h(x , x , y, y )
B (B − 1)2



(N − 1)(N − 2)
8(B − 2)
E h(x, x0 , y, y 0 )h(x00 , x000 , y 0 , y 00 )h(x0000 , x00000 , y 00 , y)
+
2
2
2
N
B (B − 1)



4
0
0
00 000
0
0000 00000
0
+ 2
E h(x, x , y, y )h(x , x , y, y )h(x , x , y, y )
B (B − 1)2
47

Li, Xie, Dai and Song



8(B − 2)
1 
= 2
E h(x, x0 , y, y 0 )h(x0 , x00 , y 0 , y 00 )h(x00 , x, y 00 , y)
2
2
B (B − 1)
N

3(N − 1) 
0
0
0 00 0 00
000 0000 00
E
h(x,
x
,
y,
y
)h(x
,
x
,
y
,
y
)h(x
,
x
,
y
,
y)
+
N2


(N − 1)(N − 2) 
0
0
00 000 0 00
0000 00000 00
+
E
h(x,
x
,
y,
y
)h(x
,
x
,
y
,
y
)h(x
,
x
,
y
,
y)
N2


1 
4
+ 2
E h(x, x0 , y, y 0 )3
2
2
B (B − 1)
N

3(N − 1) 
0
0 2
00 000
0
+
E
h(x,
x
,
y,
y
)
h(x
,
x
,
y,
y
)
N2


(N − 1)(N − 2) 
0
0
00 000
0
0000 00000
0
E h(x, x , y, y )h(x , x , y, y )h(x , x , y, y ) .
+
N2

Appendix G. Argument that ZB does not converge to Gaussian
Note that the third-order moment of ZB scales as O(B −3 ) (due to (33)), but when dividing
by its variance which scales as O(B −2 ), the skewness becomes a constant with respect to
B. Furthermore, examining the Taylor expansion of moment generating function at θ = 0,
we have
θ2
θ3
0
0
0
0 2
0 3 θZB
E[eθZB ] = 1 + E[ZB
]θ +
E[(ZB
) ] + E[(ZB
) e
] + o(θ3 ).
| {z }
2 | {z } 6
0

1

Recall that the moment generating function of a standard normal Z is given by E[eθZ ] =
1 + θ2 /2 + o(θ3 ). The difference between the two moment generating functions is given by
0

E[eθZB ] − E[eθZ ] =

|θ|3
|θ|3
0
0 3 θ 0 ZB
0 3
|E[(ZB
) e
]| + o(θ3 ) >
c|E[(ZB
) ]| + o(θ3 ),
6
6
0

(42)

0

where the inequality is due to the fact that eθ ZB > 0 and we may assume it is larger than
an absolute constant c. Note that the first term on the right hand side of (42) is given by
(cθ3 /6)Var[ZB ]−3/2 |E[ZB 3 ]|, which is clearly bounded away from zero. Hence,
0

E[eθZB ] − (1 +

θ2
|θ|3
) >
γ + o(θ3 )
2
6

for some constant γ > 0. This shows that the difference between the moment generating
0 and a standard normal is always non-zero and, hence, Z 0 does not converge
functions of ZB
B
to a standard normal in any sense.

48

