# Nonparametric Detection of Multiple Location-Scale Change Points via Wild
Binary Segmentation

Gordon J. Ross
University of Edinburgh
<EMAIL_ADDRESS>

###### Abstract

While parametric multiple change point detection has been widely studied, less
attention has been given to the nonparametric task of detecting multiple
change points in a sequence of observations when their distribution is
unknown. Most existing work on this topic is either based on penalized cost
functions which can suffer from false positive detections, or on binary
segmentation which can fail to detect certain configurations of change points.
We introduce a new approach to change point detection which adapts the
recently proposed Wild Binary Segmentation (WBS) procedure to a nonparametric
setting. Our approach is based on the use of rank based test statistics which
are especially powerful at detecting changes in location and/or scale. We show
via simulation that the resulting nonparametric WBS procedure has favorable
performance compared to existing methods, particularly when it comes to
detecting changes in scale. We apply our procedure to study a problem in
stylometry involving change points in an author’s writing style, and provide a
full implementation of our algorithm in an associated R package.

##  1 Introduction

The multiple change point problem involves the simultaneous estimation of the
number and location of $K\geq 0$ change points
$\boldsymbol{\tau}=\\{\tau_{1},\ldots,\tau_{K}\\}$ that partition a sequence
$\textbf{y}=(y_{1},\ldots,y_{n})$ into $K+1$ segments such that:

$S_{1}=\\{y_{1},\ldots,y_{\tau_{1}}\\},\quad
S_{2}=\\{y_{\tau_{1}+1},\ldots,y_{\tau_{2}}\\},\quad\ldots,\quad
S_{K+1}=\\{y_{\tau_{K}+1},\ldots,y_{n}\\}.$

We will make the commonly used assumption that the $y_{i}^{\prime}s$ are
univariate observations which are independent and identically distributed
between each pair of change points so that $y_{i}\sim_{i.i.d.}F_{j}$ if
$y_{i}\in S_{j}$ for some set of unknown continuous distributions
$F_{1},\ldots,F_{K+1}$. The literature on the single change point problem
where $K\in\\{0,1\\}$ is vast, with an overview of traditional methods
provided by Basseville and Nikiforov (1993). For the multiple change point
problem where $K$ can be greater than $1$, most existing literature has
assumed that the distributional form of each $F_{j}(\cdot;\theta_{j})$ is
known, with the change points consisting of shifts in a finite dimensional
parameter $\theta$. A recent overview of parametric techniques for multiple
change point detection is provided by Niu et al. (2016) and (Truong et al.,
2020). Most existing methods fall into one of four categories:

  1. 1.

Sequential approaches such as the CUSUM (Page, 1954; Moustakides, 1986; Lai,
1995) or change point model framework (Hawkins et al., 2003) which process the
sequence one observation at a time, and which are only concerned with
estimating the location of the most recent change point.

  2. 2.

Penalized cost function approaches which seek to find the global minimum of
some appropriately chosen cost function $C(\mathbf{y},\boldsymbol{\tau})$
(such as the negative log-likelihood) with a penalty for the number of change
points (Hawkins, 2001; Killick et al., 2012; Jackson et al., 2005; Lavielle,
2005).

  3. 3.

Stepwise approaches such as binary segmentation which transform the multiple
change point problem into the task of recursively testing for the existence of
a single change point (Bai, 1997; Cho and Fryzlewicz, 2015; Olshen et al.,
2014; Vostrikov, 1981; Inclan and Tiao, 1994; Fryzlewicz, 2014). The key
advantage of such approaches is that they allow for a strict control on the
false positive rate, i.e. the probability of incorrectly concluding that at
least one change point exists in the sequence when in fact none are present.

  4. 4.

Various Bayesian models such as (Barry and Hartigan, 1993; Green, 1995;
Fearnhead, 2006; Chib, 1998). Although interesting, the Bayesian paradigm
takes a quite different approach to the multiple change point problem, and we
will not discuss it further.

In many real world situations, it is not reasonable to expect the
distributional form of each $F_{j}$ to be known. This has led to the
development of nonparametric change point detection algorithms which do not
require information about the sequence distributions. Much of the classic
nonparametric literature has focused on the single change point setting
(Carlstein, 1988; Pettitt, 1979; Brodsky and Darkhovsky, 1993; Chen, 2019)
however the multiple change point setting is also of considerable interest.
Extensions of the sequential approach to the nonparametric setting have been
considered by authors such as (Hawkins and Deng, 2010; Ross et al., 2011;
Padilla et al., 2019), however such approaches are inefficient in a non-
sequential setting since inference for a change point at time $t$ is based
only on the observations $y_{1},\ldots,y_{t-1}$ rather than on the whole
sequence. The penalized cost function approach is limited in the nonparametric
setting since there is no likelihood to act as a cost function. However, Zou
et al. (2014) has recently shown encouraging results using a cost function
based on a nonparametric likelihood, and follow-up work by (Haynes et al.,
2017) has rendered their approach computationally tractable even for long
sequences. While this direction is promising, a limitation is the lack of
control on the false positive probability which as we show in Section 3.1 can
be very substantial.

The stepwise binary segmentation approach directly extends to the
nonparametric setting, since the test for a single change point at each stage
can be carried out using a nonparametric single change point test such as that
based on maximized Mann-Whitney statistic as used in Pettitt (1979). Such an
approach was taken by Matteson and James (2014) in a multivariate context.
This allows for a strict control on the false positive rate. However, binary
segmentation can perform poorly in multiple change point settings due to its
greedy nature which only ever searches for the best fitting single change
point at each step. In some situations, which we discuss further in Section
2.2, multiple change points will mask each other which may lead to a failure
to detect any of them.

Recently in an influential paper, (Fryzlewicz, 2014) proposed a modified
version of binary segmentation which solves some of these problems. Their
approach – named Wild Binary Segementation (WBS) – is based on maximizing a
test statistic over randomly sampled subsequences of $y_{1},\ldots,y_{n}$. By
considering subsequences, WBS is able to detect change points which only
affect local parts of the sequence which would potentially be masked when
using standard binary segmentation. Both their original paper and the
subsequent literature shows a high degree of promise for WBS.

In (Fryzlewicz, 2014), only the parametric setting of the change point problem
is considered, with a focus on detecting change points in the mean of a
sequence of Gaussian observations. In this work, we develop a nonparametric
version of their WBS procedure which is able to detect location-scale (i.e.
mean and/or variance) changes in sequences with an arbitrary and potentially
unknown distribution, while maintaining strict control of the false positive
rate. Our approach is based on nonparametric rank statistics which have a null
distribution which is independent of the form of the $F_{j}$ distributions.
This independence allows critical values of the test statistics to be computed
in a distribution-free manner. We compute these values using large scale Monte
Carlo simulation and the resulting procedure is implemented in an R package
npwbs available at https://cran.r-project.org/web/packages/npwbs/index.html

The remainder of this paper proceeds as follows: In Section 2.1 we discuss the
single change point problem where $K\in\\{0,1\\}$. Section 2.2 generalizes
this to the multiple change point problem via binary segmentation, and
explains the problems that this procedure can suffer from. In Section 2.3 we
introduce our nonparametric multiple change point detection procedure using
Wild Binary Segmentation, with further implementation details discussed in
Sections 2.4 and 2.5. Section 3 compares its performance to several existing
nonparametric change point methods. This includes the method proposed by
(Padilla et al., 2021) who also developed a nonparametric WBS scheme that
differs from ours in both the choice of test statistic and the thresholding
procedure. We show in this section that their method tends to suffer from poor
performance on most datasets. Finally in Section 4 we conclude with a real
data example taken from the field of stylometry where the task is to detect
changes in the writing style of an author.

##  2 Methodology

###  2.1 Detecting a Single Change Point

We first consider the task of testing for the existence of a single change
point $\tau$ at a specific known location $\tau=k$ in a sequence
$y_{1},\ldots,y_{n}$ of observations. This can be phrased as a two sample
hypothesis test (Basseville and Nikiforov, 1993):

$H_{0}:y_{1},\ldots,y_{n}\sim_{i.i.d}F_{0}$ (1)
$H_{1}:y_{1},\ldots,y_{k}\sim_{i.i.d}F_{1},\quad
y_{k+1},\ldots,y_{n}\sim_{i.i.d}F_{2},\quad F_{1}\neq F_{2}$

Let $T^{k}_{1,n}$ denote a two sample test statistic that has been chosen with
regards to the type of change that we wish to detect (for example, a shift in
location or scale). To make the notation clear, we write $T^{k}_{p,q}$ to
denote a two sample test statistic for the samples
$\mathbf{y}_{1}=(y_{p},y_{p+1},\ldots,y_{k})$ and
$\mathbf{y}_{2}=(y_{k+1},\ldots,y_{q})$. In the parametric case where the
functional forms of $F_{0},F_{1},F_{2}$ are all known with only their
parameters unknown, a statistic based on the generalised likelihood ratio test
is a common choice (Basseville and Nikiforov, 1993). Now suppose that $k$ is
unknown and requires estimation. In this case, the null hypothesis is the same
as above but the alternative hypothesis becomes:

$H_{1}:\exists k:y_{1},\ldots,y_{k}\sim_{i.i.d}F_{1},\quad
y_{k+1},\ldots,y_{n}\sim_{i.i.d}F_{2},\quad F_{1}\neq F_{2}$

The test statistic can be chosen as the maximal value of $T^{k}_{1,n}$ over
all possible $k$, i.e.:

$T_{1,n}=\max_{1\leq k<n}|T^{k}_{1,n}|$

The null hypothesis of no change is rejected if $T_{1,n}>\gamma$ for some
appropriately chosen threshold $\gamma(n)$, where we write . In the case where
the null distribution of $T_{1,n}$ is known and where $\gamma$ is chosen to be
equal to its $1-\alpha^{th}$ percentile then this leads to a control of the
false positive rate at level $\alpha$, i.e. if the null hypothesis of no
change is true, then the probability of incorrectly detecting a change is less
than $\alpha$.

If the null hypothesis is rejected then an estimate of $\tau$ is given by
$\hat{\tau}=\tilde{k}$ where $\tilde{k}$ is the value of $k$ for which
$|T^{k}_{1,n}|$ is maximal. There is a vast literature studying the
theoretical properties of this change detection procedure in a parametric
context, including derivations of the (usually asymptotic) null distribution
of $T_{1,n}$ under various choices of the test statistic $T^{k}_{1,n}$, along
with the construction of confidence bands for the estimate of $\tau$ .

The extension of this procedure to a nonparametric setting is
straightforwards: simply choose $T^{k}_{1,n}$ to be a test statistic with a
null distribution that does not depend on the underlying distribution of
$y_{i}$. One way to achieve this is to take $T^{k}_{1,n}$ to be a statistic
associated with a nonparametric test based on sample ranks (Pettitt, 1979;
Brodsky and Darkhovsky, 1993). In this case, the null distribution of
$T_{1,n}$ will not depend on the distribution of the $y_{i}$ observations and
the false positive rate can hence be controlled without requiring
distributional information about $y_{i}$.

###  2.2 Binary Segementation for Detecting Multiple Change Points

Next suppose that there are an unknown number $K\geq 0$ of change-points
$\boldsymbol{\tau}=\\{\tau_{1},\ldots,\tau_{K}\\}$ and the task is to estimate
both $K$ and $\boldsymbol{\tau}$. The above procedure could be directly
extended by replacing the maximal test statistic $T_{1,n}$ with one that is
maximized over every possible configuration of change points for every value
$K\in\\{1,2,\ldots,n-1\\}$. While finding the configuration which maximizes
$T_{1,n}$ might appear challenging, it can often by achieved in a
computationally efficient manner through the use of dynamic programming
(Killick et al., 2012). However, determining the null distribution of the
resulting $T_{1,n}$ is a substantially more difficult task. This has led to
the development of step-wise procedures which avoid this direct maximization
and instead recursively apply hypothesis tests based on the single change
point alternative. Specifically, consider again the single change point
hypothesis test:

${}_{0}:y_{1},\ldots,y_{n}\sim_{i.i.d}F_{0}$ (2) $H_{1}:\exists
k:y_{1},\ldots,y_{k}\sim_{i.i.d}F_{1},\quad
y_{k+1},\ldots,y_{n}\sim_{i.i.d}F_{2},\quad F_{1}\neq F_{2}$

Let $T_{1,n}=\max_{1\leq k<n}|T^{k}_{1,n}|$ for an appropriate $T^{k}_{1,n}$
be defined as above. Again, the null hypothesis of no change is rejected if
$T_{1,n}>\gamma(n)$ where we choose $\gamma$ as a function of the sample size
for reasons that will be discussed below. Suppose the null hypothesis is
rejected and let $\tilde{k}$ be the estimated change point location. The
binary segmentation procedure consists of performing a further hypothesis test
to check for a change in the observations $y_{1},\ldots,y_{\tilde{k}}$ to the
left of the change point. This is carried out in the same manner as was done
for the original sequence, i.e. by computing $T_{1,\tilde{k}}=\max_{1\leq
k<\tilde{k}}|T^{k}_{1,\tilde{k}}|$. If this
$T_{1,\tilde{k}}>\gamma(\tilde{k})$ then we deduce that there is a second
change point in the observations $y_{1},\ldots,y_{\tilde{k}}$. This procedure
is then applied recursively, with the segment $y_{1},\ldots,y_{\tilde{k}}$
split into two pieces around this second point, and both segments tested for
further change points and so on, until no more change points are found. The
same procedure is applied to the observations
$y_{\tilde{k}+1},\ldots,y_{\tilde{n}}$ to the right of the original change
point, with a further change point in this segment detected if
$T_{\tilde{k}+1,n}=\max_{\tilde{k}+1\leq
k<n}|T^{k}_{\tilde{k}+1,n}|>\gamma(n-\tilde{k})$. The final output of this
procedure is an estimate $\hat{K}$ of the number of change points and their
locations
$\hat{\boldsymbol{\tau}}=(\hat{\tau}_{1},\ldots,\hat{\tau}_{\hat{K}})$.

The binary segmentation procedure was first introduced by Vostrikov (1981) and
has been widely used since (Inclan and Tiao, 1994; Fryzlewicz, 2014). It is
sometimes claimed that the main advantage of binary segmentation is
computational since it avoids searching through every possible change point
configuration (Niu et al., 2016). While this is true, another major advantage
is that it allows the false positive rate to be strictly controlled in the
case where the sequence does not contain any change points which is not
typically the case with multiple change point methods which are based on
optimizing a penalized cost function.

Figure 1: Example of a sequence with 3 change points (dotted lines) at
locations $\boldsymbol{\tau}=\\{100,110,120\\}$ where binary segmentation will
struggle to detect any changes due to the masking effect.

Despite this, binary segmentation suffers from a major limitation. At each
stage in the recursive segmentation process, a decision is taken whether the
segment under consideration contains exactly 0 or 1 change points. This may
result in a failure to detect change in the case where two or more change
points mask each other. Consider the sequence shown in Figure 1. In this case,
there are three change points in the sequence mean that are easily visible to
the eye. However there is no obvious location for a single change point that
would split this sequence into two segments with substantially different
means. As such, performing the hypothesis test in Equation 2 based on a single
change point alternative may fail to reject the null hypothesis, in which case
the estimate of the number of change points will be $\hat{K}=0$. To solve this
problem, Fryzlewicz (2014) introduced Wild Binary Segmentation which replaces
the above test statistic with one based on maximizing over random length
subsequences of $y_{1},\ldots,y_{n}$. We will now discuss this procedure along
with its nonparametric extension.

###  2.3 Nonparametric Wild Binary Segmentation

Our goal is the nonparametric detection of multiple location-scale change
points in a univariate sequence of continuous valued observations. Following
the above discussion, we use a hypothesis testing framework with a
nonparametric two-sample test statistic that is sensitive to changes in both
location and scale. A natural choice would be an omnibus statistic such as the
Kolmogorov-Smirnov or Cramer-von-Mises which can detect arbitrary
distributional changes. However these statistics can be underpowered for
detecting location-scale changes compared to other statistics which are more
specialized towards these alternatives . Instead, we will consider a Lepage-
like statistic which is based on a linear combination of standardized Mann-
Whitney and Mood test statistics (Lepage, 1971; Ross et al., 2011).

Let $\boldsymbol{y}_{1}=(y_{1},\ldots,y_{k})$ and
$\boldsymbol{y}_{2}=(y_{k+1},\ldots,y_{n})$ denote two samples of independent
observations and consider a test for distributional equality where the null
hypothesis is that both samples have identical distributions. For each
observation $y_{i}$, let $r_{i}$ denote its rank in the combined sample, i.e.
$r_{i}=\sum_{j=1}^{n}I(y_{i}\geq y_{j})$ so that smallest observation has rank
$1$, the second smallest has rank $2$ and so on.. The Mann-Whitney test (Mann
and Whitney, 1947) is a two sample rank test for the alternative hypothesis
that the distribution of the two samples differs in location. It is defined
as:

$U^{k}_{1,n}=\sum_{i=1}^{k}r_{i}-\frac{k(k+1)}{2}$

Since this statistic depends on the observations only through their ranks, it
is easy to show that its distribution under the null hypothesis is independent
of the distribution of the observations. Specifically, its mean and variance
under the null are (Mann and Whitney, 1947):

$E[U^{k}_{1,n}]=\frac{k(n-k)}{2},\quad
Var[U^{k}_{1,n}]=\frac{k(n-k)(n+1)}{12}$

While the Mann-Whitney test will reject the null hypothesis if
$\boldsymbol{y}_{1}$ and $\boldsymbol{y}_{2}$ come from distributions with
different locations, it will generally not reject it if the distributions
differ only in scale. For this purpose an alternative test statistic should be
used. Our choice is the Mood statistic (Mood, 1954), defined as:

$M^{k}_{1,n}=\sum_{i=1}^{k}\left(r_{i}-\frac{n+1}{2}\right)^{2}$

Like the Mann-Whitney test, the distribution of the Mood statistic under the
null hypothesis is independent of the distribution of the observations. Its
mean and variance under the null hypothesis are Mood (1954):

$E[M^{k}_{1,n}]=\frac{k(n^{2}-1)}{12},\quad
Var[M^{k}_{1,n}]=\frac{k(n-k)(n+1)(n^{2}-4)}{180}$

In cases where the distribution of the two samples can differ in either
location or scale, a test statistic which is sensitive to both possibilities
can be defined by combining the Mann-Whitney and Mood statistics together.
Specifically, we use a Lepage-type statistic similar to Lepage (1971) that is
based on the sum of squared standardized versions of these statistics, i.e:

$T^{k}_{1,n}=\left(\frac{U^{k}_{1,n}-E[U^{k}_{1,n}]}{\sqrt{Var[U^{k}_{1,n}]}}\right)^{2}+\left(\frac{M^{k}_{1,n}-E[M^{k}_{1,n}]}{\sqrt{Var[M^{k}_{1,n}]}}\right)^{2}$

Note that the original Lepage statistic defined in Lepage (1971) use an
Ansari-Bradley test statistic rather than the Mood test statistic for
detecting in scale. Our decision to use the Mood is based on some empirical
results suggesting it has favorable power (Duran, 1976) although we did not
find a substantial difference between the two.

Now suppose that we wish to test the hypothesis in Equation 2 for a single
change point alternative in the sequence $y_{1},\ldots,y_{n}$. Due to the
concerns discussed in the above section, the previous maximized test statistic
$T_{1,n}=\max_{k}|T^{k}_{1,n}|$ may fail to detect a change in the case where
multiple change points mask each other. Our WBS procedure avoids this problem
by randomly selecting subsequences of $y_{1},\ldots,y_{n}$ and computing the
maximized single change point Lepage statistic on each of these subsequences.
The statistic used to test the hypothesis in Equation 2 is then the maximum of
these subsequence statistics. The intuition behind this procedure can be seen
in Figure 1. The failure of standard binary segmentation in this situation
comes from the fact that there is no single location which would divide the
sequence into two segments with substantially different means. However, there
are subsequences within which a change point will be very apparent, for
example the subsequence $y_{95},\ldots,y_{105}$ which will contain the first
change point located at $\tau_{1}=100$. As such, the value of the Lepage
statistic would be substantially higher when evaluated only on this
subsequence rather than on the sequence as a whole.

More formally, let $F^{M}_{n}$ be a set of $M$ randomly chosen intervals where
the start and end points of each interval have been sampled uniformly from the
set $\\{1,2,\ldots,n\\}$. For each such interval $[s_{m},e_{m}]\in F^{M}_{n}$,
let $T^{k}_{s_{m},e_{m}}$ denote the value of the Lepage statistic for the two
sample hypothesis test where $\mathbf{y}_{1}=(y_{s_{m}},\ldots,y_{k})$ and
$\mathbf{y}_{2}=(y_{k+1},\ldots,y_{e_{m}})$ and let
$T_{s_{m},e_{m}}=\max_{s_{m}\leq k<e_{m}}|T^{k}_{s_{m},e_{m}}|$ be the maximal
value of the test statistic on this subsequence. The test statistic used in
the WBS procedure is then:

$T^{WBS}_{1,n}=\max_{m:[s_{m},e_{m}]\in F^{M}_{n}}T_{s_{m},e_{m}}$ (3)

The null hypothesis of no change is rejected if $T^{WBS}_{1,n}>\gamma(n)$
where $\gamma$ is again a function of the number of observations and will be
discussed in the next section. If rejection occurs then let
$(s_{0},e_{0})=\arg\max_{m}T_{s_{m},e_{m}}$ be the interval for which the
Lepage statistic was maximal. The estimate of the change point is then
$\hat{\tau}=\tilde{k}$ where $\tilde{k}=\arg\max_{k}T^{k}_{s_{0},e_{0}}$. The
WBS procedure then finds additional change points in a recursive manner, using
the same procedure as binary segmentation above. First, the observations
$y_{1},\ldots,y_{\tilde{k}}$ are tested for a change point by maximizing the
Lepage statistic over $M$ randomly chosen subsequences of
$y_{1},\ldots,y_{\tilde{k}}$ and comparing it to a threshold
$\gamma(\tilde{k})$. The observations $y_{\tilde{k}+1},\ldots,y_{n}$ to the
right of the $\tilde{k}$ are tested for further change points in a similar
way. If any further change points are found, the sequence is again split into
two segments and both are tested for additional change points, until no
further change points are found. A full description of the procedure is given
in Algorithm 1.

###  2.4 Threshold Determination

We next turn our attention to choosing appropriate values of the thresholds
$\gamma(n)$ used in the WBS-Lepage procedure. In order to control the false
positive rate, we will choose $\gamma(n)$ to be the $1-\alpha^{th}$ percentile
of the null distribution of $T^{WBS}_{1,n}$ in which case the false positive
rate will be controlled at level $\alpha$. For example, if we choose
$\gamma(n)$ to be equal to the $95^{th}$ percentile then the probability of
incorrectly detecting a change when none exists will be less than $0.05$. Note
that this is a different approach to that taken by (Fryzlewicz, 2014) who
instead used an information criteria based method to determine the number of
change points, which does not give any direct control over the false positive
rate.

Due to the nonparametric nature of the Lepage test statistics, the null
distribution of $T^{WBS}_{1,n}$ does not depend on the distribution of the
$y_{i}$ observations. As such, Monte Carlo simulation can be used to find the
required percentiles of the null distribution. While this requires a
substantial amount of computation, it only needs to be done a single time and
the resulting percentiles will hold for any length $n$ sequence of independent
continuous observations regardless of their distribution.

Figure 2: Values of the threshold $\gamma(n)$ required to give a false
positive probability of 0.05 (red line) and 0.01 (black line) for WBS using
the Lepage statistic.

We implemented this computation as follows: for each
$n\in\\{10,11,12,\ldots,99,100\\}$ we simulated $1,000$ sequences of length
$n$ from a $N(0,1)$ distribution and computed $T^{WBS}_{1,n}$ for each
sequence. For each $n$, the resulting $T^{WBS}_{1,n}$ statistics were ranked
from smallest to largest, with the $50^{th}$ largest being an estimate of the
$95^{th}$ percentile of $T^{WBS}_{1,n}$. Since the distribution of
$T^{WBS}_{1,n}$ does not depend on the distribution of the data, this will
also be a valid estimate of its $95^{th}$ percentile for any other choice of
sequence distribution, not just $N(0,1)$.

The same procedure was then carried our for each
$n\in\\{105,110,\ldots,995,1000\\}$, for each
$n\in\\{1100,1200,\ldots,5000\\}$, then for each
$n\in\\{6000,7000,\ldots,10000\\}$. Linear interpolation based on these values
was then used to produce thresholds for every $n\leq 10,000$. Figure 2 shows a
plot of the resulting thresholds $\gamma(n)$ for significance level
$\alpha=0.05$. We also carried out a similar procedure to find thresholds
corresponding to a significance level of $0.01$, which may be useful in cases
where there is even greater desire to avoid false positives.

These thresholds provide the values of $\gamma(n)$ which are used to test for
an initial change point in $y_{1},\ldots,y_{n}$. Suppose a change point is
found a location $\tilde{k}$. When the test for the next change point is
applied to observations $y_{1},\ldots,y_{\tilde{k}}$ we use the threshold that
corresponds to the length of this segment, i.e. $\gamma(\tilde{k})$. This
procedure is then repeated at each stage of the segmentation, with
$\gamma(\cdot)$ always chosen to be the value corresponding to the length of
the segment being tested for a change point. The full procedure for our
nonparametric WBS algorithm is provided in Algorithm 1. This algorithm is
implemented in the npwbs R package Ross (2021) which also contains the
precomputed values of $\gamma(n)$ above for both $\alpha=0.05$ and
$\alpha=0.01$.

The following two points should be noted:

  1. 1.

Since our procedure only detects a change point if the observed value
$T^{WBS}_{p,q}$ is greater than its $95^{th}$ (or $99^{th}$) null percentile,
it may not be possible to detect a change point if the segment under
consideration is too short. This is because our use of rank-based test
statistics means that the distribution of $T^{WBS}_{p,q}$ is highly discrete
when $q-p$ is small. For example, suppose we attempt to detect changes in a
sequence of length $3$, so the observations are $y_{1},y_{2},y_{3}$.
Regardless of the observed values of the $y_{i}^{\prime}s$, there are only $6$
possible sequences that can be obtained once the observations are converted to
ranks: (1,2,3), (1,3,2), (2,1,3), (2,3,1), (3,1,2), (3,2,1). As such, it is
clear that the number of unique values of $\max_{k}|T^{k}_{1,3}|$ must be less
than 6 (it will not be exactly equal to 6, since several of these rank
sequences will produce identical values of the test statistic). As such, it
will not be possible to choose a threshold $\gamma$ such that the false
positive probability is controlled at level $\alpha=0.05$ and where it is also
possible for $\max_{k}|T^{k}_{1,3}|$ to be greater than this value. In general
for the Lepage test statistic, it can be shown that the minimum number of
observations required for both these conditions to hold is $10$. As such when
we compute thresholds and subsample sequences, we restrict the minimum length
of the sequences considered to be 10. Note that this does not mean that we
require change points to be at separated by at least 10 observations in order
to be detectable, since a detected change point in such a sequence of $10$
observations will segment it into two sequences of smaller length.

  2. 2.

An alternative to this pre-computing procedure would be to adaptively generate
each $\gamma(\cdot)$ in real-time during the segmentation procedure using a
permutation test, as was done in a CUSUM context by Olshen et al. (2014). A
potential advantage of this approach is that it allows the choice of each
$\gamma$ to be sensitive to the previous split points, and would also (e.g.)
allow the threshold level $\alpha$ to be changed adaptively, perhaps gradually
lowering it after each split in accordance with a Bonferonni-like multiple
comparisons procedure. While this has some appeal, it does require a
significant amount of computational resources which may not be practical. As
such, we use the pre-generation approach when it comes to assessing
performance in Section 3.

Algorithm 1 WBS-Lepage

1:Choose a value for $\alpha$ on which the $\gamma$ thresholds are based

2:Initialization: Call WBS-LP$(\mathbf{y},1,n)$ to find all change points in
$\mathbf{y}=(y_{1},\ldots,y_{n})$

3:

4:procedure WBS-LP($\mathbf{y},p,q$)

5: For each $y_{i}$ where $p\leq i\leq q$, replace $y_{i}$ with its rank
$r_{i}=\sum_{j=p}^{q}I(y_{i}\geq y_{j})$

6: Compute $T^{WBS}_{p,q}$ from Equation 3.

7: Let $\gamma(q-p)$ be the precomputed threshold corresponding to $\alpha$.

8: if $T^{WBS}_{p,q}>\gamma(q-p)$ then

9: Let $(s_{0},e_{0})=\arg\max_{m}T_{s_{m},e_{m}}$

10: Let $\tilde{k}=\arg\max_{k}T^{k}_{s_{0},e_{0}}$.

11: Return $\bigcup\left[\\{\text{WBS-
LP}(\mathbf{y},p,\tilde{k})\\},\\{\tilde{k}\\},\\{\text{WBS-
LP}(\mathbf{y},\tilde{k}+1,q)\\}\right]$

12: else

13: Return $\\{\varnothing\\}$

14:Prune the resulting change points as described in Section 2.5.

###  2.5 Pruning

A potential issue when using WBS is that, like binary segmentation, the
location of each change point is estimated in a greedy manner and is not
adjusted based on additional changes points that are found at later stages in
the segmentation procedure. This can sometimes result in an overestimation of
the number of change points if a single change point ends up getting detected
twice. Suppose that when WBS is used on a sequence, the first change point is
estimated to occur at location $\tilde{k}_{1}$. The sequence is then split
into two segments, and WBS is applied recursively to each. Suppose this
results in two further change points being found at locations $\tilde{k}_{0}$
and $\tilde{k}_{2}$ such that $\tilde{k}_{0}<\tilde{k}_{1}<\tilde{k}_{2}$. It
may be the case that given $\tilde{k}_{0}$ and $\tilde{k}_{2}$, the change
point at $\tilde{k}_{1}$ is not needed, i.e. that the observations between
$\tilde{k}_{0}$ and $\tilde{k}_{2}$ are identically distributed.

This issue was noticed as far back as Inclan and Tiao (1994) but hasn’t always
been taken into account in the subsequent literature. We will use a procedure
similar to that discussed in Inclan and Tiao (1994) to prune the change points
found by WBS via a post-processing step. Suppose that the WBS-Lepage algorithm
finds $\hat{K}$ change points that have been estimated to occur at locations
$\hat{\tau}_{1}<\ldots<\hat{\tau}_{K}$. For notational convenience we write
$\hat{\tau}_{0}=0$ and $\hat{\tau}_{K+1}=n$. For each $k\in 1,\ldots,K$ we
will make a decision whether $\hat{\tau}_{k}$ should be pruned. This is done
by performing the following hypothesis test, which tests whether the
observations between $\hat{\tau}_{k-1}$ and $\hat{\tau}_{k+1}$ are identically
distributed:

$H_{0}:y_{\hat{\tau_{k-1}+1}},\ldots,y_{\hat{\tau_{k+1}}}\sim_{i.i.d}F_{0}$
$H_{1}:\exists k:y_{\hat{\tau_{k-1}+1}},\ldots,y_{k}\sim_{i.i.d}F_{1},\quad
y_{k+1},\ldots,y_{\hat{\tau_{k+1}}}\sim_{i.i.d}F_{2},\quad F_{1}\neq F_{2}$

This is carried out using the same WBS procedure described above, i.e. we
compute the maximized Lepage statistic over $M$ randomly chosen subsequences
of $y_{\hat{\tau_{k-1}+1}},\ldots,y_{\hat{\tau_{k+1}}}$ to form
$T^{WBS}_{\hat{\tau}_{k-1}+1,\hat{\tau}_{k+1}}$ and compare this to the
threshold $\gamma(\hat{\tau}_{k+1}-\hat{\tau}_{k-1})$. If the null hypothesis
is not rejected, then we delete $\hat{\tau}_{k}$ from the list of change
points. Note that since this post-processing step only involves the removal of
change points rather than the creation of new ones, it cannot increase the
false positive rate of tests used in the WBS procedure. The improvement in
performance that comes from this pruning step will be studied in Section 3.3.

##  3 Experiments

We now investigate the performance of our proposed WBS-Lepage method for
detecting location-scale changes. We will compare its performance to three
other recently proposed methods for nonparametric change detection. The first
is the divisive partitioning approach from Matteson and James (2014) which is
implemented in the ecp R package (James and Matteson, 2015). This uses a type
of binary segmentation to perform change detection in multivariate sequences,
but is also applicable in our univariate setting. The second method is that of
(Zou et al., 2014) which estimates change points using a penalized cost
function approach. Specifically, they utilize a nonparametric maximum
likelihood estimate, penalized by a BIC-type penalty. While their approach is
computationally demanding, an extension using the PELT (Pruned Linear Exact
Time (Killick et al., 2012)) algorithm was introduced by Haynes et al. (2017)
which has been implemented in the changepoint.np R package (Haynes et al.,
2021) and we use this version. The third method we compare to is the PYWR
method from (Padilla et al., 2021) which also uses a WBS-type approach
although with different test statistics and thresholding than ours. We use the
R implementation which they provide alongside their paper. For our own WBS-
Lepage procedure, we choose $\gamma(\cdot)$ such that the false positive rate
is controlled at the level of $\alpha=0.05$, and utilize the pruning step
described in Section 2.5.

###  3.1 False Positives

In many situations it will not be known in advance whether the sequence being
studied contains any change points at all. In this case, it is desirable for a
change point detection algorithm to avoid false positive detections, and to
correctly estimate that no change points have occurred.

Due to the way the thresholds were constructed in Section 2.4, our WBS
algorithm is guaranteed to control the false positive rate at the desired
level (such as $\alpha=0.05$). The ECP algorithm is also based around
hypothesis testing and hence offers a similar guarantee. However approaches
such as PELT which are based only on optimizing a penalized cost function
offer no such guarantees. This is not necessarily a flaw, since the guaranteed
control of false positives reflects a strongly frequentist attitude to
statistical inference which is not universally accepted. Nonetheless, a change
detection algorithm which is highly likely to find that change points have
occurred when in fact they have not, will have limited practical utility in
some situations.

To explore this further, we carry out the following simulation: we
independently simulated $50,000$ sequences from each from the following three
distributions: 1) Normal(0,1), 2) Student-t(3), and 3) Lognormal(1,1/2). The
latter two are examples of a heavy-tailed and a skewed distribution
respectively. Each of these $50,000$ sequences consists of $100$ independent
observations, and does not contain any change points. The WBS, ECP and PELT
algorithms were then applied to each sequence. Since there are no change
points in any of the sequences, any detected change points will be false
positives. For both WBS and ECP we set the false positiver rate to be
$\alpha=0.05$. For PELT, we tried a variety of cost functions implemented in
the changepoint.np package, specifically the Bayesian Information Criterion
(BIC), the Modified BIC (MBIC) and the Strengthened Information Criterion
(SIC). For PYWR we used the default parameter settings in the authors code.

Table 1 shows the proportion of sequences in which a change point was detected
. For the WBS-Lepage and ECP procedures, the false positive rate is correctly
bounded by $\alpha=0.05$ regardless of the data distribution, i.e. there is
less than a $5\%$ probability of incorrectly flagging that at least one change
point has occurred, when the sequence actually contains no change points. This
reflects their nonparametric nature and allows the false positive rate to be
controlled regardless of the sequence distribution. While the PYWR algorithm
does not provide any guaranteed control of the false positive rate, it does in
practice also incur only a low number of false positive detections. However
the PELT approach suffers from a severely high rate of false positives, with
the BIC and SIC penalties incorrectly flagging that has a change has occurred
in $99\%$ of sequences, regardless of the data distribution. The MBIC penalty
does slightly better and only flags in around $40\%$ of cases. This makes the
PELT approach questionable in situations where it reasonable to expect that
the sequence may not contain any change points.

Model | Normal | Student-t(3) | Lognormal(1,1/2)
---|---|---|---
WBS-LP | 0.047 | 0.048 | 0.047
ECP | 0.048 | 0.049 | 0.050
PYWR | 0.028 | 0.028 | 0.027
PELT-MBIC | 0.396 | 0.394 | 0.392
PELT-BIC | 0.990 | 0.990 | 0.990
PELT-SIC | 0.990 | 0.990 | 0.990
Table 1: Probability of each change detection method producing a false
positive detection in a length $100$ sequence of independent observations with
no change points, simulated from various distributions.

###  3.2 Change Detection Performance

We now explore how our method performs on a number of simulated data sets.
Specifically, we consider the following four data models which have previously
been studied in the change point literature:

(a) fms: Normal

(b) fms: Student-t(3)

(c) fms: Lognormal(1,1/2)

(d) mix: Normal

(e) mix: Student-t(3)

(f) mix: Lognormal(1,1/2)

(g) interval: Normal

(h) interval: Student-t(3)

(i) interval: Lognormal(1,1/2)

(j) dhk: Normal

(k) dhk: Student-t(3)

(l) dhk: Lognormal(1,1/2)

(m) kfe: Normal

(n) kfe: Student-t(3)

(o) kfe: Lognormal(1,1/2)

Figure 3: Plots showing a typical sequence simulated from each of the 15 data
models. Each row respectively shows the fms, mix, interval, dhk, and kfe data.

  * •

fms data, previously studied by Fryzlewicz (2014); Frick et al. (2014). This
data consists of 497 observations with $K=6$ change points at locations
$\boldsymbol{\tau}=\\{39,226,243,300,309,333\\}$. The sequence means in each
segment are respectively $0.18,0.08,1.07,-0.53,0.16,-0.69,-0.16$ and the
standard deviation $\sigma=0.3$ is constant.

  * •

mix data, previously studied by Fryzlewicz (2014). This consists of 560
observations with $K=13$ change points at locations
$\boldsymbol{\tau}=\\{11,21,41,61,91,121,161,201,251,301,361,421,491\\}$. The
sequence means in each segment are respectively
$7,-7,6,-6,5,-5,4,-4,3,-3,2,-2,1,-1$ and the standard deviation $\sigma=4$ is
constant.

  * •

dhk data, previously studied by Davies et al. (2012). This consists of 1000
observations with $K=9$ change points at locations
$\boldsymbol{\tau}=\\{100,200,300,400,500,600,700,800,900\\}$. The sequence
standard deviations in each segment are respectively
$2.5,1,2.5,1,2.5,1,2.5,1,2.5,1$ and the sequence mean $\mu=0$ is constant.

  * •

kfe data, previously studied by Killick et al. (2012). This consists of 1000
observations with $K=5$ change points. Unlike the other data models, these
change points are not fixed, and are randomly generated from a uniform
distribution on $[30,970]$ with the constraint that there must be at least 30
observations in each segment. The mean of the sequence $\mu=0$ is constant,
and the sequence standard deviations are randomly generated from a
Lognormal($0,\log(10/2)$) distribution. Note that in the original version of
this data studied by Killick et al. (2012) there were 10 change points; we
have reduced the number to 5 since the nonparametric version of the
segmentation task is substantially more challenging than the parametric
version which they consider.

The first two data modes represent shifts in location, while the latter two
are shifts in scale. We also consider a fifth data model which has not been
previously studied and is intended to represent a short-duration location
change inside a longer sequence: which is precisely the situation in which
classic binary segmentation algorithms struggle since there is no obvious
location to place a single change point:

  * •

interval data. This data consists of 1000 observations with $K=2$ change
points at locations $\boldsymbol{\tau}=\\{490,510\\}$. The sequence means in
each segment are respectively $0,2,0$ and the standard deviation $\sigma=1$ is
constant.

For each of the five data models, we will consider the case where the error
distribution is 1) Normal, 2) Student-t with 3 degrees of freedom, and 3)
Lognormal(1,1/2), giving 15 data models in total. The latter two distributions
represent heavy-tailed and skewed distributions respectively. For each of
these 15 types of data, $100$ independent sequences were simulated. Figure 3
shows a typical sequence in each of these 15 data sets.

To measure performance, we use two metrics. The first assesses whether the
change detection algorithm can correctly recover the true number of change
points. The second assesses whether the algorithm can correctly identify their
locations. For a particular data model, let $K$ denote the true number of
change points and $\hat{K}_{i}$ denote the estimated number of change points
in the $i^{th}$ simulated data set where $i\in\\{1,2,\ldots,100\\}$.
Similarly, let $\boldsymbol{\tau}=\\{\tau_{1},\ldots,\tau_{K}\\}$ denote the
true location of the change points, and
$\boldsymbol{\hat{\tau}}_{i}=\\{\hat{\tau}_{i,1},\ldots,\hat{\tau}_{i,\hat{K}}\\}$
be their estimated locations in the $i^{th}$ data set. The two metrics are:

  1. 1.

$\frac{1}{100}\sum_{i=1}^{100}|K-\hat{K}_{i}|$, i.e. the average absolute
error in the estimated number of change points

  2. 2.

$\frac{1}{100K}\sum_{i=1}^{100}\sum_{k=1}^{K}I(\tau_{k}\in\boldsymbol{\hat{\tau}}_{i})$,
i.e. the proportion of true change points which are estimated to be in the
correct location.

The second metric only counts a change point as being correctly detected if
its location is predicted with exact accuracy. Since this is fairly
challenging for sequences where the changes are not of a large magnitude, we
also consider a third metric which counts a change point as being correctly
detected if it is estimated to have occurred within 2 observations (inclusive)
on either side, i.e.:

  1. 3.

$\frac{1}{100K}\sum_{i=1}^{100}\sum_{k=1}^{K}I(\exists\hat{\tau}_{i,j}\in\boldsymbol{\hat{\tau}}_{i}\text{
such that }|\tau_{k}-\hat{\tau}_{i,j}|\leq 2)$,

Data | Method | Normal | Student-t(3) | LNorm(1,1/2)
---|---|---|---|---
fms | WBS-LP | 0.32 | 0.37 | 0.19
ECP | 0.09 | 0.11 | 0.06
PELT-MBIC | 0.12 | 0.45 | 0.22
| PYWR | 1.92 | 1.31 | 1.43
mix | WBS-LP | 1.44 | 1.02 | 0.77
ECP | 0.95 | 0.25 | 0.55
PELT-MBIC | 1.17 | 0.86 | 0.74
| PYWR | 5.60 | 4.29 | 4.72
interval | WBS-LP | 0.21 | 0.31 | 0.15
ECP | 0.10 | 0.18 | 0.13
PELT-MBIC | 1.30 | 1.88 | 1.96
| PYWR | 2.21 | 2.36 | 2.11
dhk | WBS-LP | 0.45 | 0.52 | 0.34
ECP | 7.10 | 7.84 | 0.38
PELT-MBIC | 1.78 | 1.85 | 1.71
| PYWR | 8.55 | 8.80 | 4.41
kfe | WBS-LP | 0.97 | 1.34 | 1.08
ECP | 2.64 | 2.75 | 2.78
PELT-MBIC | 1.35 | 1.57 | 1.52
| PYWR | 3.18 | 3.03 | 3.03
Table 2: Average absolute discrepancy $\frac{1}{100}\sum_{i=1}^{100}|K-K_{i}|$
for each change point method computed over $100$ realizations of each of the
15 data models. Lower numbers indicate better performance. The best method for
each data model is highlighted in bold text.

We compare performance of nonparametric WBS to the ecp, changepoint.np (PELT)
and PYWR methods discussed above. We use the MBIC version of PELT due to its
lower false positive rate discussed above. Tables 2, 3 and 4 show the results
of all these methods evaluated over the simulated data. We can draw the
following conclusions:

  1. 1.

When it comes to the identifying the correct number of change points, the WBS
method performs the best when it comes to scale changes, while ECP is superior
for location changes. However the extent to which ECP struggles with scale
changes is notable, with it often failing to detect even a single change point
on the dfk data despite it containing 9 change points.

  2. 2.

The PYWR method performs relatively poorly on essentially every data set under
all performance metrics.

  3. 3.

The PELT method generally performs well for location changes although it
somewhat struggles on the interval data, and consistently overestimates the
number of change points. This is due to its tendency to incur false positive
when faced with a long segment of observations which does not contain any
change points.

  4. 4.

When it comes to identifying the correct location of the change points, the
ECP method performs very poorly. Its performance improves drastically when
performance is judged on estimating the change point to within two
observations of its true location. This suggests that it not a preferred
approach when absolute accuracy is necessary. The performance of WBS and PELT
is very similar when it comes to identifying the location of change points
that correspond to location changes, with WBS being superior for scale
changes.

In summary, these results show a strong performance by the WBS algorithm which
tends to be the best at detecting scale shifts, while also performing
reasonably well at detecting location shifts. While it is not the best method
in every situation, it did not perform extremely badly in any situation
either. The very poor performance of ECP in some of the scale shift settings
means that ECP cannot be recommended in these settings, and PYWR genearlly
does not perform well. However it should be noted that an advantage of ECP is
its ability to also detect changes in multivariate data, for which our WBS
approach would not be appropriate. The difference between WBS and PELT is less
drastic, with both giving similar performances for detecting location changes,
with WBS being better for scale changes. However PELT’s tendency to flag for
false positive detections in situations where there are no change points is a
potential handicap in situations where there may not be any changes at all, or
where the number of change points is small compared to the length of the
sequence, as in the interval data.

###  3.3 The Impact of Pruning

In Section 2.5 we mentioned that Wild Binary Segmentation can still suffer
from overestimating the number of change points due to its greedy nature, and
recommended a pruning strategy to combat this. The above results utilized this
pruning technique. It is potentially interesting to compare the extent to
which pruning actually improves performance. To this end, Table 5 compares the
performance of the pruned and unpruned version of the algorithm, according to
the mean absolute discrepancy performance metric. It can be seen that pruning
results in consistent and noticeable improvements in performance across most
of the data sets. As such, its use can generally be recommended in practice
compared to the unpruned version.

##  4 Real Data

We conclude with a real data example drawn from the field of stylometry, where
statistical methods are used to answer literary authorship questions.
Specifically we consider the case of Sir Terry Pratchett, a celebrated British
author who wrote the Discworld series of fantasy novels, which consists of
$41$ books written over a $32$ year period from 1983 to 2015. During 2007,
Pratchett was diagnosed with Alzheimers disease (BBC News, 2015), although he
continued to write for several years after this diagnosis. In (Ross, 2020) the
authors developed a multivariate parametric change point model to analyse
whether there was a detectable change in Pratchett’s writing style following
this diagnosis. We reanalyze their dataset in a nonparametric manner.

In stylometry, it is common to represent texts as vectors containing the
counts of frequently occurring words. Specifically, let $w_{1},\ldots,w_{200}$
denote the $200$ most commonly occurring words across the whole corpus of
Discworld books. Each of the $41$ books can then be written as a vector
$\mathbf{y}_{i}=(c_{i,1},\ldots,c_{i,200})$ where $c_{i,j}$ denotes the number
of times that word $w_{j}$ appeared in book $i$, with the books arranged in
temporal order so that the first Discworld novel corresponds to vector
$\mathbf{y}_{1}$ and the last corresponds to vector $\mathbf{y}_{41}$. A
$201^{st}$ element $c_{i,201}$ is also appended to each vector to count the
number of words in the book which were not one of the $200$ common words, so
that the elements of each vector sum to the total number of words in the
corresponding book. For reference, the list of the 200 most common words in
the corpus is shown in Figure 4, reproduced from (Ross, 2020).

Given this data representation, (Ross, 2020) developed a change point mdoel
using the Dirichlet-Multonomial distribution to test for a change point in the
sequence of books $\mathbf{y}_{1},\ldots,\mathbf{y}_{41}$, to investigate
whether Alzheimers may cause a detectable change in writing style that would
show up in the usage pattern of these common words. This hypothesis is based
on a substantial volume of stylometric evidence which has found that literary
writing style can be effectively characterized through the use of commonly
occurring words (Argamon and Levitan, 2005) . The change point approach taken
in (Ross, 2020) is multivariate and explicitly treats each book as a $201$
element vector. Since our WBS method is univariate, we will use a lower
dimensional representation of this dataset. Specifically we project the 43x201
matrix of word counts down onto the principal components of its associated
correlation matrix. We found that the first principal component did not reveal
much useful information. However the second principal is shown in Figure 5 and
displays evidence of structural change towards both the beginning and end of
the sequence.

the, and, to, a, of, I, in, he, was, it, you, that, his, said, her, she, had,
with, as, not, but, at, for, on, is, be, have, him, they, all, were, what, me,
there, one, my, this, if, no, from, so, would, up, out, by, been, them, or,
do, could, when, we, an, are, who, like, know, about, your, did, mr, then,
will, very, its, into, their, now, more, man, well, down, which, some, see,
back, time, think, just, than, dont, little, can, only, here, any, way, how,
over, thought, other, good, say, never, too, looked, much, before, come, two,
go, again, old, even, has, made, might, where, head, right, eyes, got, after,
mrs, still, yes, hand, off, something, face, should, away, through, must,
people, sir, get, though, miss, look, long, us, came, going, went, am,
himself, make, why, men, own, big, around, im, those, take, lord, seemed,
first, tell, being, always, another, quite, woman, upon, want, things,
nothing, last, door, these, didnt, such, oh, knew, once, took, great, really,
put, thing, day, young, told, voice, our, let, most, enough, thats, because,
every, room, turned, may, left, without, saw, many, course, anything, looking,
ever, asked, heard, yet, night, find, done

Figure 4: List of the 200 most common words in the Discword corpus of 41
novels.

’

(a)

(b)

Figure 5: The left hand plot shows the $41$ Discworld books projected onto the
second principal component of the corpus. The right hand plot shows the same
data with the two detected change-points superimposed as vertical lines.

We deployed our WBS-Lepage method on this univariate sequence of $41$
observations. Two change points were detected at locations $i=8$ and $i=35$
and are shown on Figure 5. The first change point most likely corresponds to a
maturation of writing style from Pratchett’s early books into his main body of
writing. The second change point immediately follows the release of his novel
‘Wintersmith’ in the year 2006. The next Discworld novel to be written after
this change point was titled ‘Making Money’ and released in 2007. Since
Pratchett’s Alzheimer’s diagnosis was publicly announced in the year 2007,
this provides some evidence that the change in writing style may indeed be
connected to the diagnosis. These results are consistent with the findings of
(Ross, 2020) who detected a similar change point using their Dirichlet-
Multnomial approach.

##  5 Conclusion

The task of detecting multiple change points in a nonparametric setting has
been under-explored, however recent developments have produced methods such as
nonparametric maximum likelihood (Zou et al., 2014) and binary segmentation
(Matteson and James, 2014) which are well-suited to this purpose. We have
presented an alternative approach that combines wild binary segmentation with
rank based hypothesis testing that competes extremely well with these existing
methods on a variety of simulated data sets. Compared to the nonparametric
maximum likelihood approach, our method is substantially better for detecting
changes in scale, while being almost as effective as detecting changes in
location. It also allows for the control of the false positive rate, which can
be important in situations where the sequence does not contain any change
points, or when the number of change points is small compared to the number of
observations. In this situation, we find that the nonparametric maximum
likelihood approach tends to overestimate the number of change points, which
may not be desirable. A drawback of our method is the amount of computational
time required to simulate from the null distribution of the test statistic in
order to generate the $\gamma(n)$ thresholds. As such, we have provided an R
implementation of our method in the npwbs package which contains our algorithm
using these precomputed thresholds.

## References

  * Argamon and Levitan [2005] S. Argamon and S. Levitan.  Measuring the Usefulness of Function Words for Authorship Attribution.  _Proceedings of the 2005 ACH/ALLC Conference,_ , 2005.
  * Bai [1997] J. Bai.  Estimating Multiple Breaks One at a Time.  _Econometric Theory_ , 13(3):315–352, 1997.
  * Barry and Hartigan [1993] D. Barry and J. Hartigan.  A Bayesian Analysis for Change Point Problems.  _Journal of the American Statistical Association_ , 88:309–319, 1993.
  * Basseville and Nikiforov [1993] M. Basseville and I. V. Nikiforov.  _Detection of Abrupt Change Theory and Application_.  Prentice Hall, 1993.
  * BBC News [2015] BBC News.  _Obituary: Sir Terry Pratchett [Online]. Available from http://www.bbc.co.uk/news/entertainment-arts-25401679_.  2015\.
  * Brodsky and Darkhovsky [1993] B. E. Brodsky and B. S. Darkhovsky.  _Nonparametric Methods in Change-Point Problems_.  Springer Netherlands, 1993.
  * Carlstein [1988] E. Carlstein.  Nonparametric Change-Point Estimation.  _Annals of Statistics_ , 16(1):188–197, 1988\.
  * Chen [2019] H. Chen.  Sequential change-point detection based on nearest neighbors.  _The Annals of Statistics_ , 47(3):1381–1407, 2019.
  * Chib [1998] S. Chib.  Estimation and comparison of multiple change-point models.  _Journal of Econometrics_ , 86(2):221–241, 1998\.  ISSN 0304-4076.
  * Cho and Fryzlewicz [2015] H. Cho and P. Fryzlewicz.  Multiple-change-point detection for high dimensional time series via sparsified binary segmentation.  _Journal of the Royal Statistical Society. Series B (Statistical Methodology)_ , 77(2):475–507, 2015.
  * Davies et al. [2012] L. Davies, C. Höhenrieder, and W. Krämer.  Recursive computation of piecewise constant volatilities.  _Computational Statistics & Data Analysis_, 56(11):3623–3631, 2012.
  * Duran [1976] B. Duran.  Survey of Nonparametric Tests for Scale.  _Communications in Statistics - Theory and Methods_ , 5(14):1287–1312, 1976.
  * Fearnhead [2006] P. Fearnhead.  Exact and efficient Bayesian inference for multiple changepoint problems.  _Statistics and Computing_ , 16:203–213, 2006.
  * Frick et al. [2014] K. Frick, A. Munk, and H. Sieling.  Multiscale change point inference.  _Journal of the Royal Statistical Society Series B_ , 76(3):495–580, 2014.
  * Fryzlewicz [2014] P. Fryzlewicz.  Wild binary segmentation for multiple change-point detection.  _The Annals of Statistics_ , 42(6):2243–2281, Dec. 2014.
  * Green [1995] P. Green.  Reversible Jump Markov Chain Monte Carlo Computation and Bayesian Model Determination.  _Biometrika_ , 84(2):711–732, 1995.
  * Hawkins [2001] D. M. Hawkins.  Fitting multiple change-point models to data.  _Computational Statistics & Data Analysis_, 37:323–341, 2001.
  * Hawkins and Deng [2010] D. M. Hawkins and Q. Deng.  A Nonparametric Change-Point Control Chart.  _Journal of Quality Technology_ , 42(2):165–173, 2010.
  * Hawkins et al. [2003] D. M. Hawkins, P. H. Qiu, and C. W. Kang.  The Changepoint Model for Statistical Process Control.  _Journal of Quality Technology_ , 35(4):355–366, 2003.
  * Haynes et al. [2017] K. Haynes, P. Fearnhead, and I. A. Eckley.  A computationally efficient nonparametric approach for changepoint detection.  _Statistics and Computing_ , 27(5):1293–1305, 2017.
  * Haynes et al. [2021] K. Haynes, R. Killick, P. Fearnhead, I. Eckley, and D. Grose.  changepoint.np: Methods for Nonparametric Changepoint Detection, 2021.  URL https://CRAN.R-project.org/package=changepoint.np.
  * Inclan and Tiao [1994] C. Inclan and G. C. Tiao.  Use of Cumulative Sums of Squares for Retrospective Detection of Changes of Variance.  _Journal of the American Statistical Association_ , 89(427):913–923, 1994.
  * Jackson et al. [2005] B. Jackson, J. Scargle, D. Barnes, S. Arabhi, A. Alt, P. Gioumousis, E. Gwin, P. Sangtrakulcharoen, L. Tan, and T. T. Tsai.  An algorithm for optimal partitioning of data on an interval.  _IEEE Signal Processing Letters_ , 12(2):105–108, 2005.
  * James and Matteson [2015] N. A. James and D. S. Matteson.  ecp: An R Package for Nonparametric Multiple Change Point Analysis of Multivariate Data.  _Journal of Statistical Software_ , 62(1):1–25, 2015.
  * Killick et al. [2012] R. Killick, P. Fearnhead, and I. A. Eckley.  Optimal Detection of Changepoints With a Linear Computational Cost.  _Journal of the American Statistical Association_ , 107(500):1590–1598, 2012.
  * Lai [1995] T. L. Lai.  Sequential Changepoint Detection in Quality Control and Dynamical Systems.  _Journal of the Royal Statistical Society. Series B - Methodological_ , 57(4):613–658, 1995.
  * Lavielle [2005] M. Lavielle.  Using penalized contrasts for the change-point problem.  _Signal Processing_ , 85(8):1501–1510, 2005.
  * Lepage [1971] Y. Lepage.  Combination of Wilcoxian’s and Ansari-Bradley’s Statistics.  _Biometrika_ , 58(1):213–217, 1971.
  * Mann and Whitney [1947] H. B. Mann and D. R. Whitney.  On a Test of Whether One of Two Random Variables is Stochastically Larger than the Other.  _Annals of Mathematical Statistics_ , 18(1):50–60, 1947.
  * Matteson and James [2014] D. S. Matteson and N. A. James.  A Nonparametric Approach for Multiple Change Point Analysis of Multivariate Data.  _Journal of the American Statistical Association_ , 109(505):334–345, 2014.
  * Mood [1954] A. Mood.  On the Asymptotic Efficiency of Certain Nonparametric Two-Sample Tests.  _Annals of Mathematical Statistics_ , 25:514–533, 1954\.
  * Moustakides [1986] G. V. Moustakides.  Optimal Stopping Times for Detecting Changes in Distributions.  _Annals of Statistics_ , 14(4):1379–1387, 1986\.
  * Niu et al. [2016] Y. S. Niu, N. Hao, and H. Zhang.  Multiple Change-Point Detection: A Selective Overview.  _Statistical Science_ , 31(4):611–623, 2016.
  * Olshen et al. [2014] A. B. Olshen, E. S. Venkatraman, R. Lucito, and M. Wigler.  Circular binary segmentation for the analysis of array-based DNA copy number data.  _Biostatistics_ , 5(4):557–572, 2014.
  * Padilla et al. [2019] O. H. M. Padilla, A. Athey, A. Reinhart, and J. G. Scott.  Sequential Nonparametric Tests for a Change in Distribution: An Application to Detecting Radiological Anomalies.  _Journal of the American Statistical Association_ , 114(526):514–528, 2019.
  * Padilla et al. [2021] O. H. M. Padilla, Y. Yu, D. Wang, and A. Rinaldo.  Optimal nonparametric change point analysis.  _Electronic Journal of Statistics_ , 15(1):1154–1201, 2021.
  * Page [1954] E. S. Page.  Continuous Inspection Schemes.  _Biometrika_ , 41(1/2):100–115, 1954.
  * Pettitt [1979] A. N. Pettitt.  A Non-Parametric Approach to the Change-Point Problem.  _Journal of the Royal Statistical Society Series C - Applied Statistics_ , 28(2):126–135, 1979.
  * Ross [2020] G. J. Ross.  Tracking the evolution of literary style via Dirichlet–multinomial change point regression.  _Journal of the Royal Statistical Society: Series A (Statistics in Society)_ , 183(1):149–167, 2020.
  * Ross [2021] G. J. Ross.  npwbs: Nonparametric Multiple Change Point Detection Using Wild Binary Segmentation, 2021.  URL https://cran.r-project.org/web/packages/npwbs/index.html.
  * Ross et al. [2011] G. J. Ross, D. K. Tasoulis, and N. M. Adams.  Nonparametric Monitoring of Data Streams for Changes in Location and Scale.  _Technometrics_ , 53(4):379–389, 2011.
  * Truong et al. [2020] C. Truong, L. Oudre, and N. Vayatis.  Selective review of offline change point detection methods.  _Signal Processing_ , 167:107299, 2020.
  * Vostrikov [1981] L. Vostrikov.  Detecting ’Disorder’ in Multidimensional Random Processes.  _Soviet Mathematics Doklady_ , 24:55–59, 1981.
  * Zou et al. [2014] C. Zou, G. Yin, L. Feng, and Z. Wang.  Nonparametric maximum likelihood approach to multiple change-point problems.  _The Annals of Statistics_ , 42(3):970–1002, 2014\.

Data | Method | Normal | Student-t(3) | LNorm(1,1/2)
---|---|---|---|---
fms | WBS.LP | 0.63 | 0.73 | 0.70
ECP | 0.12 | 0.13 | 0.14
PELT.MBIC | 0.64 | 0.74 | 0.76
| PYWR | 0.14 | 0.15 | 0.18
mix | WBS-LP | 0.48 | 0.62 | 0.60
ECP | 0.09 | 0.06 | 0.07
PELT-MBIC | 0.51 | 0.61 | 0.62
| PYWR | 0.00 | 0.00 | 0.00
interval | WBS-LP | 0.57 | 0.75 | 0.70
ECP | 0.11 | 0.10 | 0.08
PELT-MBIC | 0.57 | 0.74 | 0.71
| PYWR | 0.28 | 0.61 | 0.56
dhk | WBS.LP | 0.29 | 0.21 | 0.24
ECP | 0.02 | 0.01 | 0.10
PELT.MBIC | 0.28 | 0.15 | 0.25
| PYWR | 0.01 | 0.00 | 0.12
kfe | WBS.LP | 0.31 | 0.30 | 0.34
ECP | 0.04 | 0.05 | 0.04
PELT.MBIC | 0.31 | 0.24 | 0.31
| PYWR | 0.04 | 0.04 | 0.06
Table 3: Percentage of change points which were estimated to be in their correct location, for each detector, averaged over $100$ realizations of each of the 15 data models. Higher numbers indicate better performance. The best method for each data model is highlighted in bold text. Data | Method | Normal | Student-t(3) | LNorm(1,1/2)
---|---|---|---|---
fms | WBS-LP | 0.90 | 0.95 | 0.93
ECP | 0.88 | 0.94 | 0.94
PELT-MBIC | 0.91 | 0.94 | 0.96
| PYWR | 0.57 | 0.66 | 0.62
mix | WBS-LP | 0.76 | 0.88 | 0.86
ECP | 0.82 | 0.90 | 0.85
PELT-MBIC | 0.79 | 0.86 | 0.86
| PYWR | 0.42 | 0.56 | 0.52
interval | WBS-LP | 0.95 | 0.99 | 0.96
ECP | 0.93 | 0.98 | 0.97
PELT-MBIC | 0.93 | 0.99 | 0.94
| PYWR | 0.41 | 0.77 | 0.70
dhk | WBS-LP | 0.70 | 0.58 | 0.64
ECP | 0.09 | 0.04 | 0.57
PELT-MBIC | 0.70 | 0.43 | 0.65
| PYWR | 0.02 | 0.00 | 0.23
kfe | WBS-LP | 0.59 | 0.56 | 0.62
ECP | 0.38 | 0.35 | 0.35
PELT-MBIC | 0.57 | 0.49 | 0.57
| PYWR | 0.15 | 0.17 | 0.19
Table 4: Percentage of change points which were estimated to be within 2 observations of their correct location, for each detector, averaged over $100$ realizations of each of the 15 data models. Higher numbers indicate better performance. The best method for each data model is highlighted in bold text. Data | Method | Normal | Student-t(3) | LNorm(1,1/2)
---|---|---|---|---
fms | WBS-Unpruned | 0.48 | 0.52 | 0.41
WBS-Pruned | 0.32 | 0.37 | 0.19
mix | WBS-Unpruned | 1.34 | 1.04 | 1.16
WBS-Pruned | 1.44 | 1.02 | 0.77
interval | WBS-Unpruned | 0.24 | 0.31 | 0.20
WBS-Pruned | 0.21 | 0.71 | 0.15
dnk | WBS-Unpruned | 0.62 | 0.71 | 0.60
WBS-Pruned | 0.45 | 0.52 | 0.34
kfe | WBS-Unpruned | 1.10 | 1.39 | 1.05
WBS-Pruned | 1.07 | 1.38 | 1.05
Table 5: Average discrepancy $\frac{1}{100}\sum_{i=1}^{100}(K-K_{i})$ for the
pruned and unpruned version of the WBS method across each data set. Lower
numbers indicate better performance.