Epidemic changepoint detection
in the presence of nuisance changes

Julius Juodakis, Stephen Marsland
School of Mathematics and Statistics, Victoria University of Wellington, New Zealand
Corresponding author: [email protected]

Abstract

Many time series problems feature epidemic changes – segments where a parameter deviates from a background baseline. The number and location of such changes can be estimated in a principled way by existing detection methods, providing that the background level is stable and known. However, practical data often contains nuisance changes in background level, which interfere with standard estimation techniques. Furthermore, such changes often differ from the target segments only in duration, and appear as false alarms in the detection results.

To solve these issues, we propose a two-level detector that models and separates nuisance and signal changes. As part of this method, we developed a new, efficient approach to simultaneously estimate unknown, but fixed, background level and detect epidemic changes. The analytic and computational properties of the proposed methods are established, including consistency and convergence. We demonstrate via simulations that our two-level detector provides accurate estimation of changepoints under a nuisance process, while other state-of-the-art detectors fail. Using real-world genomic and demographic datasets, we demonstrate that our method can identify and localise target events while separating out seasonal variations and experimental artefacts.

Keywords: Change point detection, Epidemic change points, Piecewise stationary time series, Segmentation, Stochastic gradient methods

1 Introduction

The problem of identifying when the probability distribution of a time series changes – changepoint detection – has been studied since the middle of the 20th century. Early developments, such as Shewhart’s control charts (Shewhart, 1930) and Page’s CUSUM test (Page, 1954), stemmed from operations research. However, as automatic sensing systems and continuous data collection became more common, many new use cases for changepoint detection have arisen, such as seismic events (Li et al., 2016), epidemic outbreaks (Texier et al., 2016), physiological signals (Vaisman et al., 2010), market movements (Bracker and Smith, 1999), and network traffic spikes (Hochenbaum et al., 2017). Stimulated by such practical interest, the growth of corresponding statistical theory has been rapid, as reviewed in e.g. Aminikhanghahi and Cook (2016) and Truong et al. (2020).

Different applications pose different statistical challenges. If a single drastic change may be expected, such as when detecting machine failure, the goal is to find a method with minimal time to detection and a controlled false alarm rate (Lau and Tay, 2019). In many other applications the number of changepoints is unknown a priori; the challenge is then to identify the change locations regardless, preferably in a computationally efficient way. Some problems, such as peak detection in sound (Aminikhanghahi and Cook, 2016) or genetic data (Hocking et al., 2017), feature epidemic segments – changepoints followed by a return to the background level – and incorporating this constraint can improve detection or simplify post-processing of ouputs.

However, current detection methods that do incorporate a background level assume it to be stable throughout the dataset (e.g. Fisch et al., 2018; Zhao and Yau, 2019). This is not realistic in several common applications. In genetic analyses such as measurements of protein binding along DNA there may be large regions where the observed background level is shifted due to structural variation in the genome or technical artefacts (Zhang et al., 2008). Similarly, a standard task in sound processing is to detect speech in the presence of dynamic background chatter, so-called babble noise (Pearce et al., 2000). In various datasets from epidemiology or climatology, such as wave height measurements (Killick et al., 2012), seasonal effects are observed as recurring background changes and will interfere with detection of shorter events. Methods that assume a constant background will be inaccurate in these cases, while ignoring the epidemic structure entirely would cost detection power and complicate the interpretation of outputs.

Our goal is to develop a general method for detecting epidemic changepoints in the presence of nuisance changes in the background. Furthermore, we assume that the only available feature for distinguishing the two types of segments is their duration: this would allow analysing the examples above, which share the property that the nuisance process is slower. The most similar research to ours is that of Lau and Tay (2019) for detecting failure of a machine that can be switched on, and thus undergo an irrelevant change in the background level. However, the setting there concerned a single change with specified background and nuisance distributions; in contrast, we are motivated by the case when multiple changes may be present, with only duration distinguishing their types. Detection then requires two novel developments: 1) rapid estimation of local background level, 2) modelling and distinguishing the two types of potentially overlapping segments.

These developments are presented in this paper as follows: after a background section we provide a new algorithm that simultaneously detects epidemic changepoints and estimates the unknown background level (Section 3). The convergence and consistency of this algorithm are proved in the Supplementary Material. While this algorithm is of its own interest, we use it to build a detector that allows local variations in the background, i.e., nuisance changes, in Section 4, again with proven properties. In Section 5 we test various properties of the algorithms with simulations before showing how the proposed nuisance-robust detector can be applied in two problems: detecting histone modifications in human genome, while ignoring structural variations, and detecting the effects of the COVID-19 pandemic in Spanish mortality data, robustly to seasonal effects. Compared to state-of-the-art methods, the proposed detector produced lower false-alarm rates (or more parsimonious models), while retaining accurate detection of true signal peaks.

2 Background

In the general changepoint detection setup, the data consists of a sequence of observations $x_{0:n}=\{x_{0},x_{1},\dots,x_{n}\}$ , split by $k$ changepoints $0<\tau_{1}<\tau_{2}<\dots\tau_{k}<n$ into $k+1$ segments. Within each segment the observations are drawn from a distribution specific to that segment; these are often assumed to belong to some fixed family $f_{\theta}(x)$ , with different values of parameter $\theta$ for each segment. The most common and historically earliest example is the change in mean of a Gaussian (see Truong et al. (2020)), i.e., for each $t\in[\tau_{i},\tau_{i+1})$ , $x_{t}\sim\mathcal{N}(\mu_{i},\sigma^{2})$ , for known, fixed $\sigma^{2}$ . (We assume $\theta\in\mathbb{R}^{1}$ to keep notation clearer, but multidimensional problems are also common, such as when both the mean and variance of a Gaussian are subject to change.)

While early research focused on optimal hypothesis testing of a single change against no change in a parameter, modern models aim to estimate the position of all changepoints $\{\tau_{i}\}$ in a given time series, and optionally the vector of segment parameters $\bm{\theta}$ . A common approach is to use a penalised cost: choose a segment cost function $C(x_{a:b};\theta)$ , usually based on negative log-likelihood of the observations $x_{a:b}$ given $\theta$ , and a penalty $p(k)$ for the number of changepoints $k$ . The full cost of $x_{\tau_{0}:\tau_{k}}$ (where $\tau_{0}=0$ and $\tau_{k}=n$ ) is then:

F(n;\bm{\tau},\bm{\theta},k)=\sum_{i=0}^{k-1}C(x_{\tau_{i}+1:\tau_{i+1}};\theta_{i})+p(k).

(1)

Changepoint number and positions are estimated by finding $F(n)=\min F(n;\bm{\tau},\bm{\theta},k)$ . Such estimation has been shown to be consistent for a range of different data generation models (Fisch et al., 2018; Zheng et al., 2019).

For a discrete problem such as this, computation of the true minimum is not trivial – a naïve brute force approach would require $\mathcal{O}(2^{n})$ tests. Approaches to reducing this cost fall in two broad classes: 1) simplifying the search algorithm by memoisation and pruning of paths (Jackson et al., 2005; Killick et al., 2012; Rigaill, 2015); 2) using greedy methods to find approximate solutions faster (Fryzlewicz, 2014; Baranowski et al., 2019). In both classes, there are methods that can consistently estimate multiple changepoints in linear time under certain conditions.

The first category is of more relevance here. It is based on the Optimal Partitioning (OP) algorithm (Jackson et al., 2005). Let the data be partitioned into discrete blocks $B_{i}:\bigcup_{i}B_{i}=\{x\}_{n}$ , so $B_{i}\cap B_{j}=\emptyset,\forall i\neq j$ . A function $V$ that maps each set of blocks $P_{j}=\{B_{i}\}$ to a cost is block-additive if:

\forall P_{1},P_{2},V(P_{1}\cup P_{2})=V(P_{1})+V(P_{2}).

(2)

If each segment incurs a fixed penalty $\beta=p(k)/k$ , then the function $F$ defined in (1) is block-additive over segments, and can be defined recursively as:

F(s)=\min_{t}\left(F(t)+C(x_{t+1:s})+\beta\right).

In OP, this cost is calculated for each $s\leq n$ , and thus its minimisation requires $\mathcal{O}(n^{2})$ evaluations of $C$ . Furthermore, when the cost function $C$ is such that for all $a\leq b<c$ :

C(x_{a:c})\geq C(x_{a:b})+C(x_{b+1:c})

(3)

then at each $s$ , it can be quickly determined that some candidate segmentations cannot be “rescued” by further segments, and so they can be pruned from the search space. This approach reduces the complexity to $\mathcal{O}(n)$ in the best case. It gave rise to the family of algorithms called PELT, Pruned Exact Linear Time (Killick et al., 2012; Zhao and Yau, 2019). Note that OP and PELT do not rely on any probabilistic assumptions and find the exact same minimum as the global exponential-complexity search. Since the introduction of PELT, a number of alternative pruning schemes have been developed (Maidstone et al., 2016).

A variation on the basic changepoint detection problem is the identification of epidemic changepoints, when a change in regime appears for a finite time and then the process returns to the background level. The concept of a background distribution $f_{B}$ is introduced, and the segments are defined by pairs of changepoints $s_{i},e_{i}$ , outside of which the data is drawn from the background model. The data model then becomes:

f(x_{t})=\begin{cases}f_{S}(x_{t};\theta_{i})&\text{ if }\exists i:s_{i}\leq t\leq e_{i}\\ f_{B}(x_{t})&\text{ otherwise}.\end{cases}

(4)

An equivalent of the CUSUM test for the epidemic situation was first proposed by Levin and Kline (1985), and a number of methods for multiple changepoint detection in this setting have been proposed (Olshen et al., 2004; Fisch et al., 2018; Zhao and Yau, 2019). These use a cost function that includes a separate term $C^{0}$ for the background points:

F(n;\{(s_{i},e_{i})\}_{k},\bm{\theta},k)=\sum_{i=1}^{k}C(x_{s_{i}:e_{i}};\theta_{i})+C^{0}(\{x_{t}:t\notin\bigcup_{i}[s_{i};e_{i}]\};\theta_{0})+p(k)

(5)

A common choice for the background distribution is some particular “null” case of the segment distribution family, so that $f_{B}(x)=f_{S}(x;\theta_{0})$ and $C^{0}(\cdot)=C(\cdot;\theta_{0})$ . However, while the value of $\theta_{0}$ is known in some settings (such as when presented with two copies of DNA), in other cases it may need to be estimated. Since this parameter is shared across all background points, the cost function is no longer block-additive as in (2), and algorithms such as OP and PELT cannot be directly applied.

One solution is to substitute the unknown parameter with some robust estimate of it, based on the unsegmented series $x_{0:n}$ . The success of the resulting changepoint estimation then relies on this estimate being sufficiently close to the true value, and so the non-background data fraction must be small (Fisch et al., 2018). This is unlikely to be true in our motivating case, when nuisance changes in the background level are possible.

Another option is to define

F(n;\theta_{0})=\min_{k,\{(s_{i},e_{i})\}_{k},\{\theta_{i}\}_{k}}F(n;\{(s_{i},e_{i})\}_{k},\{\theta_{i}\}_{k},\theta_{0},k)

which can be minimised using OP or PELT, and then use gradient descent for the outer minimisation over $\theta_{0}$ (Zhao and Yau, 2019). The main drawback of this approach is an increase in computation time proportional to the number of steps needed for the optimisation.

3 Detection of changepoints with unknown background parameter

To estimate the background parameter and changepoints efficiently, while allowing a large proportion of non-background data, we introduce Algorithm 1. The algorithm makes a pass over the data, during which epidemic segments are detected in a standard way, i.e., by minimising penalized cost using OP (steps 3-12). The estimate of parameter $\theta_{0}$ is iteratively updated at each time point. We will show that this process converges almost surely to the true $\theta_{0}$ for certain models of background and segment levels, or reaches its neighbourhood $(\theta_{0}-\epsilon;\theta_{0}+\epsilon)$ under more general conditions.

Algorithm 1 Detection of changepoints with unknown background parameter

1: Input:

C^{0},C,\beta,x_{0:n}

2: Initialize

F(0)=0,B(0)=\{x_{0}\},\theta_{0}=x_{0}

3: for

t\in 1,\dots,n

4: Calculate

F_{B}=F(t-1)+C^{0}(x_{t};\theta_{0})

, and

F_{S}=\min_{1\leq k\leq l}F(t-k)+C(x_{t-k+1:t})+\beta

5: if

F_{B}<F_{S}

then

6: Assign

B(t)=B(t-1)\cup x_{t}

, and recalculate

\theta_{0}

from

B

F(t)=F_{B}

8: else

9: Assign

B(t)=B(t-k)

(here

k

\operatorname*{argmin}F_{S}

in step 4)

10:

F(t)=F_{S}

11: end if

12: end for

13: Repeat steps 3-12, without updating

\theta_{0}

14: Output:

\theta_{0}

, changepoint positions

The algorithm then makes a second pass over the data (step 13), repeating the segmentation using the final estimate of $\theta_{0}$ . The step is identical to a single application of a known-background changepoint detector, such as the one described in Fisch et al. (2018). Its purpose is to update the changepoint positions that are close to the start of the data and so had been determined using less precise estimates of $\theta_{0}$ . This simplifies the theoretical performance analysis, but an attractive option is to use this algorithm in an online manner, without this step. We evaluate the practical consequences of this omission later in Section 5.

By segmenting and estimating the background simultaneously, this algorithm is computationally more efficient than dedicated optimisation over possible $\theta_{0}$ values as in Zhao and Yau (2019), while allowing a larger fraction of non-background points than using a robust estimator on unsegmented data, as in Fisch et al. (2018).

We now demonstrate that Algorithm 1 will converge almost surely under certain constraints, and is consistent.

3.1 Convergence

Theorem 1.

Consider the problem of an epidemic change in mean, with data $x_{0:n}$ generated as in (4). Assume the distribution $f_{B}(x)$ and marginal distribution of segment points $f_{S}(x)=\int f_{S}(x;\theta)f_{\Theta}(\theta)d\theta$ are symmetric and strongly unimodal, with unknown background mean $\theta_{0}$ , and that data points within each segment are iid. Here, $f_{\Theta}$ is the probability density of parameter $\theta$ corresponding to each data point. Denote by $w_{t}$ the estimate of $\theta_{0}$ obtained by analysing $x_{0:t}$ by Algorithm 1. The sequence $\{w_{t}\}$ converges:

1.

to $\theta_{0}$ almost surely, if $\mathbb{E}f_{S}=\theta_{0}$ .
2.

to a neighbourhood $(\theta_{0}-\epsilon,\theta_{0}+\epsilon)$ almost surely, where $\epsilon\rightarrow 0$ as the number of background points $n$ between successive segments $n\rightarrow\infty$ .

We refer the reader to Appendix A for the proof. It is based on a result by Bottou (1998), who established conditions under which an online minimisation algorithm almost surely converges to the optimum (including stochastic gradient descent as a special case). We show that these conditions are either satisfied directly by the updating process in our algorithm, or longer update cycles can be defined that will satisfy the conditions, in which case the convergence efficiency is determined by $n$ .

3.2 Consistency

The changepoint model can be understood as a function over an interval that is sampled to produce $n$ observed points. As the sampling density increases, more accurate estimation of the number and locations of changepoints is expected; this property is formalised as consistency of the detector. Fisch et al. (2018) showed that detectors based on minimising penalised cost are consistent in the Gaussian data case, and their result can be adapted to prove consistency of Algorithm 1. We will use the strengthened SIC-type penalty $\alpha\log(n)^{1+\delta},\delta>0$ as presented in Fisch et al. (2018), but a similar result can be obtained with $\delta=0$ .

Additionally, the minimum signal-to-noise ratio at which segments can be detected is formalised as this assumption:

\forall i,(e_{i}-s_{i})\Delta_{i}>\log(n)^{1+\delta},

(6)

where $\Delta_{i}$ is a function of the relative strength of the parameter changes $|\mu_{i}-\mu_{0}|$ and $\sigma_{k}/\sigma_{0}+\sigma_{0}/\sigma_{k}$ .

Theorem 2.

Let the data $x_{0:n}$ be generated from an epidemic changepoint model as in (4), with $f_{B}$ and $f_{S}$ Gaussian, and the changing parameter $\theta$ is either its mean or variance (assume the other parameter is known). Further, assume (6) holds for $k$ changepoints. Analyse the data by Algorithm 1 with penalty $\beta=\alpha\log(n)^{1+\delta}$ , $\alpha,\delta>0$ . The estimated number and position of changepoints will be consistent, i.e. $\forall\epsilon>0,n>B$ :

P\left(\hat{k}=k,|\hat{s_{i}}-s_{i}|<\frac{A}{\Delta_{i}}\log(n)^{1+\delta},|\hat{e_{i}}-e_{i}|<\frac{A}{\Delta_{i}}\log(n)^{1+\delta},\forall 1\leq i\leq k\right)\geq 1-Cn^{-\epsilon},

(7)

for some $A,B,C$ that do not increase with $n$ .

The proof is given in Appendix B. We use the connection between our Algorithm 1 and stochastic gradient descent to establish error bounds on the backgorund parameter estimates. These bounds then allow us to apply a previous consistency result (Fisch et al., 2018) to our case.

3.3 Pruning

Much of the improvement in performance of change estimation algorithms comes from pruning of the search space. Standard pruning (Killick et al., 2012) can be applied to Algorithm 1. It continues to find optimal solutions in this case, as is shown in Proposition 1.

To implement pruning, the search set for $F_{S}$ in step 4 would be changed to:

F_{S}=\min_{t^{\prime}\in K}(F(t^{\prime})+C(x_{t^{\prime}+1:t})+\beta),

and step 11a would be introduced to update $K$ as:

\text{11a. }K=K\cap\{s:F(s)+C(x_{s+1:t})<F(t),~{}t+1-l\leq s<t\}\cup\{t\}.\\

Proposition 1.

Assume a cost function $C$ such that (3) applies. If, for some $s>t-l$ :

F(s)+C(x_{s+1:t})\geq F(t),

then $t-k$ will not be a segment start point in the optimal solution, and can be excluded from consideration for all subsequent $t^{\prime}>t$ .

Proof.

For all $t^{\prime}\leq s+l$ , the proof applies as for other pruned algorithms (Killick et al., 2012):

F(t)+C(x_{t+1:t^{\prime}})+\beta\leq F(s)+C(x_{s+1:t})+C(x_{t+1:t^{\prime}})+\beta\leq F(s)+C(x_{s+1:t^{\prime}})+\beta.

For $t^{\prime}>s+l$ , segment $(s,t^{\prime})$ will exceed the length constraint and thus cannot be part of the final segmentation. ∎

4 Detecting changepoints with a nuisance process

4.1 Problem setup

In this section, we consider the changepoint detection problem when there is an interfering nuisance process. We assume that this process, like the signal, consists of segments, and we index the start and end times as $s^{N}_{j},e^{N}_{j}$ . Data within these segments is generated from a nuisance-only distribution $f_{N}$ , or from some distribution $f_{NS}$ if a signal occurs at the same time. In total, four states are possible (background, nuisance, signal, or nuisance and signal), so the overall distribution of data points is:

f(x_{t})=\begin{cases}f_{NS}(x_{t};\theta_{i},\theta^{N}_{j})&\text{ if }\exists i,j:t\in[s_{i};e_{i}]\cap[s^{N}_{j},e^{N}_{j}]\\ f_{S}(x_{t};\theta_{i})&\text{ if }\exists i:t\in[s_{i},e_{i}],t\notin\cup_{j}[s^{N}_{j},e^{N}_{j}]\\ f_{N}(x_{t};\theta^{N}_{j})&\text{ if }\exists j:t\in[s^{N}_{j},e^{N}_{j}],t\notin\cup_{i}[s_{i},e_{i}]\\ f_{B}(x_{t})&\text{ otherwise}\end{cases}

(8)

We add two more conditions to ensure identifiability: 1) the nuisance process evolves more slowly than the signal process, so $\min(e^{N}_{j}-s^{N}_{j})>\max(e_{i}-s_{i})$ ; 2) signal segments are either entirely contained within a nuisance segment, or entirely out of it:

\forall i,j,\text{ either }[s_{i},e_{i}]\subset(s^{N}_{j},e^{N}_{j}),\text{ or }[s_{i},e_{i}]\cap[s^{N}_{j},e^{N}_{j}]=\emptyset

(9)

We adapt the changepoint detection method as follows. Define $X_{S}=\bigcup_{i}x_{s_{i}:e_{i}},X_{N}=\bigcup_{j}x_{s^{N}_{j}:e^{N}_{j}}$ , a penalty $p^{\prime}(m)=\beta^{\prime}m$ for the number of nuisance segments, and cost functions $C^{NS},~{}C^{S},~{}C^{N},~{}C^{0}$ corresponding to each of the distributions in (8). Then the full cost of a model with $k$ true segments and $m$ nuisance segments is:

F(n;\{(s_{i},e_{i})\}_{k},k,\{(s^{N}_{j},e^{N}_{j})\}_{m},m,\bm{\theta})=C^{0}(x_{0:n}\setminus(X_{S}\cup X_{N}))+\sum_{i=0}^{k}C^{S}(x_{s_{i}:e_{i}}\setminus X_{N};\theta_{i})+\\ +\sum_{j=0}^{m}\left(C^{N}(x_{s^{N}_{j}:e^{N}_{j}}\setminus X_{S};\theta_{j})+\sum_{i=0}^{k}C^{NS}(x_{s^{N}_{j}:e^{N}_{j}}\cap x_{s_{i}:e_{i}};\theta_{i},\theta_{j})\right)+\beta k+\beta^{\prime}m.

(10)

Note that this cost can also be expressed, using $k_{j}$ as the number of signal segments that overlap a nuisance segment $j$ and $k_{0}=k-\sum k_{j}$ , as:

F(n;\{(s_{i},e_{i})\}_{k},k,\{(s^{N}_{j},e^{N}_{j})\}_{m},m,\bm{\theta})=C^{0}(x_{0:n}\setminus(X_{S}\cup X_{N}))+\sum_{i=0}^{k_{0}}C^{S}(x_{s_{i}:e_{i}}\setminus X_{N};\theta_{i})+\\ +\beta k_{0}+\sum_{j=0}^{m}(C^{\prime}(x_{s^{N}_{j}:e^{N}_{j}})+\beta k_{j}+\beta^{\prime}).

Condition (9) ensures that $C^{\prime}$ and $C^{S}$ are independent (no points or $\theta_{i}$ are shared between them), so $F$ is block-additive over (signal or nuisance) segments and can be minimised by dynamic programming.

4.2 Proposed method

To minimise this cost, we propose the method outlined in Algorithm 2. Its outer loop proceeds over the data to identify segments by the usual OP approach. However, to evaluate the cost $C^{\prime}(x_{a:b})$ , which allows segment and nuisance overlaps, a changepoint detection problem with unknown background level must be solved over $x_{a:b}$ . This is achieved by an inner loop using Algorithm 1. If $C^{\prime}=C^{S}$ , this method would reduce to a standard detector of epidemic changepoints, with the exception that segments are separated into two types depending on their length; this is very similar to the collective and point anomaly (CAPA) detector in Fisch et al. (2018).

The cost $C^{\prime}$ is minimised using Algorithm 1, so by Theorem 2, the number of positions of true segments will be estimated consistently given accurate assignment of the $m$ nuisance positions $s^{N}_{j},e^{N}_{j}$ . However, the latter event is subject to a complex set of assumptions on relative signal strength, position, and duration of the segments. Therefore, we do not attempt to describe these in full here, but instead investigate the performance of the method by extensive simulations in Section 5.1.3.

Algorithm 2 Adaptive segmentation by changepoint detection

1: Input:

C^{0},C^{S},C^{N},C^{NS},\beta,\beta^{\prime},x_{0:n}

2: Initialize

F(0)=0,\theta_{0}=x_{0}

, lists of tuples

chp_{S},chp_{N},chp_{NS}

3: If

\theta_{0}

not known: estimate using an appropriate robust estimator

4: for

t\in 1,\dots,n

5: for

t^{\prime}\in 0,\dots,t-l-1

6: Apply Algorithm 1 to data

x_{t^{\prime}+1:t}

, store returned cost as

C^{\prime}

and changepoints as

chp_{NS}(t^{\prime})

F_{N}(t^{\prime})=F(t^{\prime})+C^{\prime}+\beta^{\prime}

8: end for

F_{B}=F(t-1)+C^{0}(x_{t})

10:

F_{S}=\min_{1\leq k\leq l}F(t-k)+C^{S}(x_{t-k+1:t})+\beta

11:

F_{N}=\min F_{N}(t^{\prime})

12: Assign

F(t)=\min\{F_{B},F_{S},F_{N}\}

13: if

F(t)=F_{B}

then

14: Assign

chp_{S}(t)=chp_{S}(t-1),chp_{N}(t)=chp_{N}(t-1)

15: else if

F(t)=F_{S}

then

16: Assign

chp_{S}(t)=chp_{S}(t-k)\cup(t-k,t)

with

k

determined by

\operatorname*{argmin}F_{S}

chp_{N}(t)=chp_{N}(t-1)

17: else

18: Assign

chp_{S}(t)=chp_{S}(t^{\prime})\cup chp_{NS}(t^{\prime})

chp_{N}(t)=chp_{N}(t^{\prime})\cup(t^{\prime},t)

with

t^{\prime}

determined by

\operatorname*{argmin}F_{N}

19: end if

20: end for

21: Output: changepoint positions

chp_{S}(n)=\{(s_{i},e_{i})\}

and

chp_{N}(n)=\{(s^{N}_{j},e^{N}_{j})\}

Algorithm 2 is stated assuming that a known or estimated value of the parameter $\theta_{0}$ , corresponding to the background level without the nuisance variations, is available. In practice, it may be known when there is a technical noise floor or a meaningful baseline that can be expected after removing the nuisance changes. Alternatively, $\theta_{0}$ may be estimated by a robust estimator if a sufficient fraction of background points can be assumed; in either case, the (known or estimated) value is substituted to reduce the computational cost. Otherwise, the method can be modified to estimate $\theta_{0}$ simultaneously with segmentation, using a principle similar to Algorithm 1.

4.3 Pruning

In the proposed method, the estimation of parameter $\mu_{j}$ (the mean of segment $j$ ) is sensitive to the segment length, therefore the cost $C^{\prime}$ does not necessarily satisfy the additivity condition (3), and so it cannot be guaranteed that PELT-like pruning will be exact. However, we can establish a local pruning scheme that retains the exact optimum with probability $\rightarrow 1$ as $n\rightarrow\infty$ .

Proposition 2.

Assume data $x_{0:n}$ is generated from a Gaussian epidemic changepoint model, and that the distance between changepoints is bounded by some function $A(n)$ :

\forall i,j:s^{N}_{j+1}-e^{N}_{j}>A(n),|s_{i}-s^{N}_{j}|>A(n).

At time $t$ , the solution space is pruned by removing:

\bm{k}_{pr,t}=\{k:F(k-1)+C^{\prime}(x_{k:t})\geq\min_{m}F(m-1)+C^{\prime}(x_{m:t})+\alpha\log(n)^{1+\delta}\}.

(11)

Here $m\in(t-A(n);t],k\in(t-A(n);t],k\neq m$ . Then $\forall\epsilon>0$ , there exist constants $B,n_{0}$ , such that when $n>n_{0}$ , the true nuisance segment positions are retained with high probability:

P(\forall j:s^{N}_{j}\notin\bigcup_{t}\bm{k}_{pr,t},e^{N}_{j}\notin\bigcup_{t}\bm{k}_{pr,t})\geq 1-Bn^{-\epsilon}.

The proof is given in Appendix C. The assumed distance bound serves to simplify the detection problem: within each window $(t-A(n),t]$ , at most 1 true changepoint may be present, and the initial part of Algorithm 2 is identical to a standard epidemic changepoint detector. It can be shown that other candidate segmentations in the pruning window are unlikely to have significantly lower cost than the one associated with $s^{N}_{j},e^{N}_{j}$ , and therefore $s^{N}_{j},e^{N}_{j}$ are likely to be retained in pruning.

As this scheme only prunes within windows of size $A(n)$ , and setting large $A(n)$ may violate the true data generation model, it is less efficient than PELT-type pruning. However, assuming that the overall estimation of nuisance and signal changepoints is consistent, Proposition 2 extends to standard pruning over the full dataset. We show that this holds empirically in Section 5.1.3.

5 Experiments

In this section, we present the results of simulations used to evaluate the performance of the methods, and demonstrations on real-world datasets.

The methods proposed in this paper are implemented in R 3.5.3, and available on GitHub at https://github.com/jjuod/changepoint-detection. The repository also includes the code for recreating the simulations and data analyses presented here.

5.1 Simulations

5.1.1 Algorithm 1 estimates the background parameter consistently

In the first experiment we tested the performance of the algorithm with a fixed background parameter. Three different sets of data were generated and used:

Scenario 1. Gaussian data with one signal segment: $n$ data points were drawn as $x_{t}\sim\mathcal{N}(\theta_{t},1)$ , with a segment of $\theta_{t}=3$ at times $t\in(0.3n;0.5n]$ , and $0$ otherwise (background level).

Scenario 2. Gaussian data with multiple signal segments: $n$ data points were drawn as $x_{t}\sim\mathcal{N}(\theta_{t},1)$ , with:

\theta_{t}=\begin{cases}-1\text{ when }t\in(0.2n;0.3n]\cup(0.7n;0.8n]\\ 1\text{ when }t\in(0.5n;0.6n]\\ 0\text{ otherwise }.\end{cases}

Scenario 3. Heavy tailed data with one signal segment: $n$ data points were drawn from the generalized $t$ distribution as $x_{t}\sim T(3)+\theta_{t}$ , with $\theta_{t}=2$ at times $t\in(0.2n;0.6n]$ and $0$ otherwise (background level).

To evaluate the background level ( $\theta_{0}$ ) estimation by Algorithm 1, we generated time series for each of the three scenarios with values of $n$ between 30 and 750, with 500 replications for each $n$ . Note that the segment positions remain the same for all sample sizes, so the increase in $n$ could be interpreted as denser sampling of the time series.

Each time series was analysed using Algorithm 1 to estimate $\theta_{0}$ . The maximum segment length was set to $l=0.5n$ and the penalty to $\beta=3\log(n)^{1.1}$ . The cost was computed using the Gaussian log-likelihood function with known variance. This means that the cost function was mis-specified in Scenario 3, and so provides a test of the robustness of the algorithm (although the variance was set to $3$ , as expected for $T(3)$ ).

For comparison, in each replication we also retrieved the median of the entire time series $x_{1:n}$ , which is used as the estimate of $\theta_{0}$ by Fisch et al. (2018). For scenarios (1) and (2), we also computed the quantiles of $\mathcal{N}(0,\sigma^{2}/\sqrt{n_{B}})$ as the oracle efficiency limit based on the CLT, with $n_{B}$ the total number of background points: $n_{B}=0.8n$ for (1) and $0.7n$ for (2).

Figure 1 shows the results for all three scenarios. It can be seen that Algorithm 1 produced consistent estimates given sufficiently long time series. When the signal-to-noise ratio was large (scenario 1), most of the estimated background values were accurate even at small $n$ , and at larger $n$ our estimator reached efficiency and accuracy similar to the oracle estimator. Predictably, the median of the non-segmented data produced biased estimates, which nonetheless may be preferrable in some cases, as our estimate showed large variability at low $n$ in scenarios 2 and 3. At $n>400$ , our estimator provides the lowest total error in all tested scenarios, even with the mis-specified cost function in scenario 3.

Refer to caption — Figure 1: Consistency of the background parameter values estimations. Time series simulated in three different scenarios were analysed by Algorithm 1 (shown in red). Lines are the inter-quartile range (solid) and 2.5-97.5% range (faint) of the background parameter estimates observed in 500 replications at each $n$ . For comparison, we also show the ranges of estimates obtained by the median (blue) and mean of true background points (green; an oracle estimator).

5.1.2 Segment positions estimated by Algorithm 1 are accurate and consistent

The same simulation setup was used to evaluate the consistency of estimated segment number and positions. Data were generated and analysed with Algorithm 1 as in Section 5.1.1, but we also retrieved the changepoint positions that were estimated in step 12. This corresponds to the online usage of the algorithm, in which segmentation close to the start of the data is based on early $\theta_{0}$ estimates and may therefore itself be less precise. We extracted the mean number of segments reported by the algorithm and also calculated the true positive rate (TPR) as the fraction of simulations in which the algorithm reported at least 1 changepoint within $0.05n$ points of each true changepoint.

The results show that the TPR approaches 1 for all three scenarios (Table 1). As expected, detecting a strong segment was easier than multiple weak ones: segmentation in scenario 1 was accurate even at $n=30$ , while $n>500$ was needed to reach $\geq 90\%$ TPR in scenario 2. In scenario 3, the algorithm correctly detected changes at the true segment start and end (TPR $\approx 100\%$ ), but the algorithm tended to fit the segment as multiple ones, likely due to the heavy tails of the $t$ distribution. Notably, skipping step 13 had very little impact on the performance of the algorithm, suggesting that the faster online version can be safely used.

Table 1: Consistency of the number and position of estimated changepoints. Time series simulated in three different scenarios were analysed using Algorithm 1. For each

n

, 500 replications were performed. The mean number of reported segments and the TPR (fraction of iterations when a changepoint was detected within

0.05n

of each true changepoint) are shown. The number of true segments was 1, 3, and 1 for scenarios 1, 2, 3 respectively. The same time series were also analysed in an online manner, i.e., without Step 13 of the algorithm; the results are shown on the right of the table.

		(Full algorithm)		(no step 13)
Scenario	$n$	Mean # segm.	TPR	Mean # segm.	TPR
1	30	1.112	0.932	1.138	0.924
	90	1.040	0.998	1.056	0.998
	180	1.054	1.000	1.060	1.000
	440	1.032	1.000	1.040	1.000
	750	1.028	1.000	1.034	1.000
2	30	0.552	0.000	0.584	0.000
	90	1.096	0.008	1.050	0.004
	180	1.762	0.086	1.720	0.080
	440	2.900	0.824	2.854	0.782
	750	3.010	0.992	3.008	0.988
3	30	0.644	0.142	0.662	0.132
	90	1.426	0.592	1.416	0.570
	180	1.930	0.878	1.932	0.866
	440	2.798	0.998	2.798	0.996
	750	3.692	1.000	3.664	1.000

5.1.3 Algorithm 2 recovers true signal segments under interference

To evaluate Algorithm 2, we generated time series with both signal and nuisance segments under two different scenarios:

Scenario 1. Gaussian data with a signal segment overlapping a nuisance segment: $n$ data points were drawn as $x_{t}\sim\mathcal{N}(\theta^{S}_{t}+\theta^{N}_{t},1)$ , with $\theta^{S}_{t}=2$ when $t\in(0.3n;0.5n]$ , 0 otherwise, and $\theta^{N}_{t}=2$ when $t\in(0.2n;0.7n]$ , 0 otherwise.

Scenario 2. Gaussian data with a nuisance segment and two non-overlapping signal segments: $n$ data points were drawn as $x_{t}\sim\mathcal{N}(\theta^{S}_{t}+\theta^{N}_{t},1)$ , with:

	$\displaystyle\theta^{N}_{t}=1\text{ when }t\in(0.2n;0.4n],~{}0\text{ otherwise}$
	$\displaystyle\theta^{S}_{t}=\begin{cases}3\text{ when }t\in(0.5n;0.6n]\\ -3\text{ when }t\in(0.7n;0.8n]\\ 0\text{ otherwise }\end{cases}$

Time series were generated for $n$ between 30 and 240, in 500 replications at each $n$ . Each series was analysed by three methods. Algorithm 2 (proposed) was run with penalties $\beta=\beta^{\prime}=3\log(n)^{1.1}$ as before. For classical detection of epidemic changes in mean, we used the R package anomaly (Fisch et al., 2018); the implementation allows separating segments of length 1, but we treated them as standard segments and set all penalties to $3\log(n)^{1.1}$ . In both methods, background parameters were set to $\mu_{0}=0,\sigma_{0}=1$ , and maximum segment length was $l=0.33n$ in scenario 1 and $l=0.15n$ in scenario 2. As an example of a different approach, we included the narrowest-over-threshold detector implemented in R package not, with default parameters. This is a non-epidemic changepoint detector that was shown to outperform most comparable methods (Baranowski et al., 2019). Since it does not include the background-signal distinction, we define signal segments as regions between two successive changepoints where the mean exceeds $\mu_{0}\pm\sigma_{0}$ .

As before, evaluation metrics were the mean number of segments and the TPR. For the proposed method, we required accurate segment type detection as well, i.e., a true positive is counted when the detector reported a nuisance start/end within $0.05n$ of each nuisance changepoint and a signal start/end within $0.05n$ of each signal changepoint. In scenario 1, we also extracted the estimate of $\theta$ corresponding to the detected signal segment (or segment closest to $(0.3n;0.5n]$ if multiple were detected). The average of these over the 500 replications is reported as $\hat{\theta}^{S}$ .

To evaluate the effects of pruning, Algorithm 2 was applied without pruning, or with global pruning as in (11), with $m\in(0;t-l)$ at each $t$ . Only 5 out of 5000 runs (500 iterations $\times$ 10 settings) showed any differences between the globally-pruned and non-pruned methods (Supplementary Table S1 in Appendix D), so we only present results obtained with pruning from here on.

We observed that the proposed method (with pruning) successfully detected true signal segments in both scenarios (Figure 2 and Table 2). The number of nuisance detections was accurate in scenario 1, and slightly underestimated in favour of more signal segments in scenario 2, most likely because the simulated nuisance length was close to the cutoff of $0.15n$ .

Table 2: The true positive rate of changepoint estimation by the proposed Algorithm 2 (pruned), and two alternative detectors: anomaly and not. Data simulated for two scenarios, in 500 replications for each

n

. The true positive rate is the fraction of iterations when a changepoint of correct type was detected within

0.05n

of each true changepoint (types are signal, S, and nuisance, N).

Scenario	$n$	Proposed (S)	Proposed (N)	anomaly	not
1	30	0.444	0.548	0.382	0.428
	60	0.742	0.752	0.532	0.724
	100	0.938	0.930	0.870	0.912
	160	0.984	0.972	0.982	0.994
	240	1.000	0.986	1.000	0.998
2	30	0.700	0.110	0.934	0.152
	60	0.924	0.238	0.962	0.872
	100	0.986	0.412	0.998	0.988
	160	0.998	0.640	1.000	1.000
	240	0.998	0.764	1.000	1.000

As expected, when nuisance segments are not included in the model (anomaly and not methods), they were identified as multiple changepoints; as a result, the number of segments was over-estimated up to 3-fold. These models are also unable to capture the signal-specific change in mean $\theta^{S}$ : anomaly estimated $\hat{\theta}^{S}=3.99$ , and not estimated $\hat{\theta}^{S}=3.95$ in scenario 1 at $n=240$ . These values correspond to the sum of the signal and nuisance effects. While the estimation is accurate and could be used to recover the value of interest by post-hoc analysis, our proposed method estimated the signal-specific change directly, as $\hat{\theta}^{S}=2.00$ in scenario 1 at $n=240$ .

5.2 Real-world Data

5.2.1 ChIP-seq

As an example application of the algorithms proposed in this paper, we demonstrate peak detection in chromatin immunoprecipitation sequencing (ChIP-seq) data. The goal of ChIP-seq is to identify DNA locations where a particular protein of interest binds, by precipitating and sequencing bound DNA. This produces a density of binding events along the genome that then needs to be processed to identify discrete peaks. Typically, some knowledge of the expected peak length is available to the researcher, although with considerable uncertainty (Rashid et al., 2011). Furthermore, the background level may contain local shifts of various sizes, caused by sequencing bias or true structural variation in the genome (Zhang et al., 2008). The method proposed in this paper is designed for such cases, and can potentially provide more accurate and more robust detection.

We used datasets from two ChIP-seq experiments, investigating histone modifications in human immune cells. Broad Institute H3K27ac data was obtained from UCSC Genome Browser, GEO accession GSM733771, as mean read coverage in non-overlapping windows of 25 bp. While ground truth is not available for this experiment, we also retrieved an input control track for the same cell line from UCSC (GEO accession GSM733742). We analysed a window near the centromere of chromosome 1, between 120,100,000 to 120,700,000 bp. This window contains nuisance variation in the background level, as seen in the input control (Figure 3, top). To improve runtime, data was downsampled to approximately 1000 points, each corresponding to mean coverage in a window of 500 bp. All read coordinates in this paper correspond to hg19 genome build.

The second dataset was the UCI/McGill dataset, obtained from https://archive.ics.uci.edu/ml/datasets/chipseq. This dataset was previously used for evaluating peak detection algorithms (Hocking et al., 2017, 2018). Mean read coverage at 1 bp resolution is provided, as well as peak annotations based on visual inspection. The annotations are weak labels in the sense that they indicate peak presence or absence in a region, not their exact positions, to acknowledge the uncertainty when labelling visually. From this data, we used the H3K36me3 modification in monocyte sample ID McGill0104, AM annotation. Around the labels $L=\{(s_{i},e_{i})\}$ in each chromosome, we extracted read coverage for the window between $s-(e-s)$ and $e+(e-s)$ bp, $s=\min s_{i}$ , $e=\max e_{i}$ , and downsampled to about 1000 points as before. Based on visual inspection of the coverage, we chose a group of labels in chromosome 12, which provides a variety of annotations and a structure that appears to contain both nuisance shifts and signal peaks.

The ChIP-seq datasets were analysed by the method proposed here (Algorithm 2), an epidemic detector from package anomaly, and the non-epidemic detector not. The length of the signal segments was limited to 50 kbp in anomaly and the proposed method. As an estimate of global mean $\mu_{0}$ , the median of the unsegmented data was used, and $\sigma_{0}$ estimated by the standard deviation of non-segmented data. As before, penalties were set to $3\log(n)^{1.1}$ . Only segments with estimated $\theta>\mu_{0}$ are shown, as we are a priori interested only in regions of increased binding.

We also used GFPOP, implemented in R package PeakSegDisk (Hocking et al., 2018): this detector has been developed specifically for ChIP-seq data processing, and models a change in Poisson rate parameter. This method does not include a single background level, but enforces alternating up-down constraints. It is intended as a supervised method, with the penalty value $\lambda$ chosen based on the best segmentation provided by the training data. Therefore, for the Broad Institute dataset, we repeated the segmentation with $\lambda\in\{10^{1},10^{2},10^{3},10^{4},10^{5},10^{6}\}$ , and show the $\lambda$ that produces between 2 and 10 segments.

To evaluate the results quantitatively, we calculated the SIC based on each method’s segmentation. We used Gaussian likelihood with parameters matching the detection status (i.e., estimated mean for the points in each segment $\hat{\mu}_{0}$ for points outside segments, and $\hat{\sigma}_{0}$ as the standard deviation for all points). The number of parameters was set to $3k+2$ : a mean and two endpoints for each of the $k$ segments reported, and two for the background parameters.

In the H3K27ac data, all methods detected the three most prominent peaks, but produced different results for smaller peaks and more diffuse change areas (Figure 3, bottom). Both PeakSegDisk and anomaly marked a broad segment in the area around 120,600,000 bp. Based on comparison with the control data, this change is spurious, and it exceeds the 50 kbp bound for target segments. While this bound was provided to the anomaly detector, it does not include an alternative way to model these changes, and therefore still reports one or more shorter segments. In contrast, our method accurately modelled the area as a nuisance segment with two overlapping sharp peaks.

Using not, the data was partitioned into 10 segments. By defining segments with low mean ( $\theta<\hat{\mu}_{0}+\hat{\sigma}_{0}$ ) as background, we could reduce this to 4 signal segments; while this removed the spurious background change, it also discarded the shorter change around 120,200,000 bp, which fits the definition of a signal peak ( $<50$ kbp) and was retained by the proposed method. This data illustrates that choosing the post-processing required for most approaches is not trivial, and can have a large impact on the results. In contrast, the parameters required for our method have a natural interpretation and may be known a priori or easily estimated, and the outputs are provided in a directly useful form. This data illustrates that choosing the post-processing options is not trivial, and can have a large impact on the results. In contrast, the parameters required for our method have a natural interpretation and may be known a priori or easily estimated, and the outputs are provided in a directly useful form.

In the UCI data, segment detections also generally matched the visually determined labels. However, our method produced the most parsimonious models to explain the changes. Figure 4 shows such an example from chromosome 12, where our method reported two nuisance segments and a single sharp peak around 62,750,000 bp. The nuisance segments correspond to broad regions of mean shift, which were also detected by anomaly and not, but using 6 and 16 segments, respectively. Notably, PeakSeg differed considerably: as this method does not incorporate a single background level, but requires segments to alternate between background and signal, the area around 62,750,000 bp was defined as background, despite having a mean of $4.5~{}\hat{\mu}_{0}$ . In total, 12 segments were reported by this method. This shows that the ability to separate nuisance and signal segments helps produce more parsimonious models, and in this way minimises the downstream efforts such as experimental replication of the peaks.

The visual annotations provided for this region are shown in the first row in Figure 4. Note that they do not distinguish between narrow and broad peaks (single annotations in this sample range up to 690 kb in size). Furthermore, comparison with such labels does not account for finer segmentation, coverage in the peak area, or the number of false alarms outside it. For these reasons we are unable to use the labels in a quantitative way.

Quantitative comparison of the segmentations by SIC also favours our proposed method in both datasets. In the Broad dataset, SICs corresponding to segmentations reported by PeakSeg, not and anomaly were 4447.4, 4532.2, and 4311.6, respectively, while the segmentation produced by our model had an SIC of 4285.8. The smallest criterion values in UCI data were also produced by our method (4664.5), closely followed by anomaly (4664.8), while not segmentation resulted in an SIC of 4705.7 and PeakSeg 6886.7. This suggests that in addition to the practical benefits of separating unwanted segments, the nuisance-signal structure provides a better fit to these datasets than models that allow only one type of segments.

5.2.2 European mortality data

The recent pandemic of coronavirus disease COVID-19 prompted a renewed interest in early outbreak detection and quantification. In particular, analysis of mortality data provided an important resource for guiding the public health responses to it (e.g. Baud et al., 2020; Yuan et al., 2020; Dowd et al., 2020).

We analysed Eurostat data of weekly deaths in Spain over a three year period between 2017 and 2020. Data was retrieved from https://ec.europa.eu Data Explorer. Besides the impact of the pandemic, mortality data contains normal seasonal variations, particularly in older age groups. We use the 60-64 years age group in which these trends are visually clear and thus provide a ground truth.

Data were analysed with the four methods introduced earlier. For the proposed and anomaly methods, we used the median and standard deviation of the first 52 weeks of the dataset as estimates of $\mu_{0}$ and $\sigma_{0}$ , respectively. Penalties in both were set to $3\log(n)^{1.1}$ , as previously. The maximum length of signal segments was set to 10 weeks, to separate seasonal effects. In addition, not was used with default parameters, defining signal as regions where the mean exceeds $\mu_{0}\pm\sigma_{0}$ , and PeakSegDisk with a basic grid search to select the penalty as before.

The results of the four analyses are shown in Figure 5. Three of the methods, anomaly, PeakSeg, and Algorithm 2, detected a sharp peak around the pandemic period. However, anomaly and PeakSeg also marked one winter period as a signal segment, while ignoring the other two. Four segments were created by not, including a broad peak continuing well past the end of the pandemic spike. In contrast, the proposed method marked the pandemic spike sharply, while also labelling all three winter periods as nuisance segments. The resulting detection using our method is again parsimonious and flexible: if only short peaks are of interest, our method reports those with lower false alarm rate than the other methods, but broader segments are also marked accurately and can be retrieved if relevant.

As in the ChIP-seq data, comparing the results by SIC identifies our method as optimal for this dataset. The values corresponding to PeakSeg, not and anomaly models were 1629.3, 1648.2, and 1626.3 respectively, while Algorithm 2 produced a segmentation with SIC 1568.2. Note that the SIC penalizes both signal and nuisance segments, so in this case our model still appears optimal despite having more parameters.

6 Discussion

In this paper, we have presented a pair of algorithms for improving detection of epidemic changepoints. Similarly to stochastic gradient descent, the iterative updating of the background estimate in Algorithm 1 leads to fast convergence while allowing large fraction of non-background points. This is utilised in Algorithm 2 to analyse nuisance-signal overlaps.

The computational complexity of both algorithms presented here is $\mathcal{O}(n)$ in the best case, which is similar to state-of-the-art pruned algorithms (Killick et al., 2012; Hocking et al., 2018). However, note that this is stated in the number of required evaluations of $C$ . It is usually implicitly assumed that this function can be evaluated and minimised over $\theta$ recursively, so that the total number of operations may also be linear. This would not be achievable with methods that estimate the background level strictly offline, such as aPELT-profile (Zhao and Yau, 2019). Therefore, development of Algorithm 1 was essential to create the overlap detector.

One of major practical benefits of the proposed model is the ability to separate non-target segments. We anticipate that this will greatly improve downstream processing, effectively reducing the false alarm rate or the manual load if the detections are reviewed. Despite that, it is difficult to evaluate this benefit at present: while there are recent datasets with annotations specifically for testing changepoint detection (Hocking et al., 2017; van den Burg and Williams, 2020), they are based on labelling all visually apparent changes. In future work, we expect to provide further application-specific comparisons that would measure the impact of separating and neutralising the nuisance process.

7 Acknowledgements

This work was supported by the Royal Society of New Zealand Marsden Fund and Te Pūnaha Matatini, a New Zealand Centre of Research Excellence.

References

Aminikhanghahi and Cook (2016) Aminikhanghahi, S. and Cook, D. J. (2016) A survey of methods for time series change point detection. Knowledge and Information Systems, 51, 339–367.
Baranowski et al. (2019) Baranowski, R., Chen, Y. and Fryzlewicz, P. (2019) Narrowest-over-threshold detection of multiple change points and change-point-like features. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 81, 649–672.
Baud et al. (2020) Baud, D., Qi, X., Nielsen-Saines, K., Musso, D., Pomar, L. and Favre, G. (2020) Real estimates of mortality following COVID-19 infection. The Lancet Infectious Diseases, 20, 773.
Bottou (1998) Bottou, L. (1998) Online algorithms and stochastic approximations. In Online Learning and Neural Networks (ed. D. Saad). Cambridge, UK: Cambridge University Press. Revised, oct 2012.
Bracker and Smith (1999) Bracker, K. and Smith, K. L. (1999) Detecting and modeling changing volatility in the copper futures market. Journal of Futures Markets, 19, 79–100.
van den Burg and Williams (2020) van den Burg, G. J. J. and Williams, C. K. I. (2020) An evaluation of change point detection algorithms. Preprint arXiv:2003.06222.
Dowd et al. (2020) Dowd, J. B., Andriano, L., Brazel, D. M., Rotondi, V., Block, P., Ding, X., Liu, Y. and Mills, M. C. (2020) Demographic science aids in understanding the spread and fatality rates of COVID-19. Proceedings of the National Academy of Sciences, 117, 9696–9698.
Fisch et al. (2018) Fisch, A. T. M., Eckley, I. A. and Fearnhead, P. (2018) A linear time method for the detection of point and collective anomalies. Preprint arXiv:1806.01947.
Fryzlewicz (2014) Fryzlewicz, P. (2014) Wild binary segmentation for multiple change-point detection. The Annals of Statistics, 42, 2243–2281.
Harvey et al. (2019) Harvey, N. J. A., Liaw, C., Plan, Y. and Randhawa, S. (2019) Tight analyses for non-smooth stochastic gradient descent. In Proceedings of the Thirty-Second Conference on Learning Theory (eds. A. Beygelzimer and D. Hsu), vol. 99 of Proceedings of Machine Learning Research, 1579–1613. Phoenix, USA: PMLR.
Hochenbaum et al. (2017) Hochenbaum, J., Vallis, O. S. and Kejariwal, A. (2017) Automatic anomaly detection in the cloud via statistical learning. Preprint arXiv:1704.07706.
Hocking et al. (2017) Hocking, T. D., Goerner-Potvin, P., Morin, A., Shao, X., Pastinen, T. and Bourque, G. (2017) Optimizing ChIP-seq peak detectors using visual labels and supervised machine learning. Bioinformatics, 33, 491–499.
Hocking et al. (2018) Hocking, T. D., Rigaill, G., Fearnhead, P. and Bourque, G. (2018) Generalized functional pruning optimal partitioning (gfpop) for constrained changepoint detection in genomic data. Preprint arXiv:1810.00117.
Jackson et al. (2005) Jackson, B., Scargle, J., Barnes, D., Arabhi, S., Alt, A., Gioumousis, P., Gwin, E., San, P., Tan, L. and Tsai, T. T. (2005) An algorithm for optimal partitioning of data on an interval. IEEE Signal Processing Letters, 12, 105–108.
Killick et al. (2012) Killick, R., Fearnhead, P. and Eckley, I. A. (2012) Optimal detection of changepoints with a linear computational cost. Journal of the American Statistical Association, 107, 1590–1598.
Lau and Tay (2019) Lau, T. S. and Tay, W. P. (2019) Quickest change detection in the presence of a nuisance change. IEEE Transactions on Signal Processing, 67, 5281–5296.
Levin and Kline (1985) Levin, B. and Kline, J. (1985) The cusum test of homogeneity with an application in spontaneous abortion epidemiology. Statistics in Medicine, 4, 469–488.
Li et al. (2016) Li, S., Cao, Y., Leamon, C., Xie, Y., Shi, L. and Song, W. (2016) Online seismic event picking via sequential change-point detection. In 2016 54th Annual Allerton Conference on Communication, Control, and Computing (Allerton). IEEE.
Maidstone et al. (2016) Maidstone, R., Hocking, T., Rigaill, G. and Fearnhead, P. (2016) On optimal multiple changepoint algorithms for large data. Statistics and Computing, 27, 519–533.
Olshen et al. (2004) Olshen, A. B., Venkatraman, E. S., Lucito, R. and Wigler, M. (2004) Circular binary segmentation for the analysis of array-based DNA copy number data. Biostatistics, 5, 557–572.
Page (1954) Page, E. S. (1954) Continuous inspection schemes. Biometrika, 41, 100.
Pearce et al. (2000) Pearce, D., Hirsch, H.-g. and Ericsson Eurolab Deutschland Gmbh (2000) The aurora experimental framework for the performance evaluation of speech recognition systems under noisy conditions. In ISCA ITRW ASR2000, 29–32.
Purkayastha (1998) Purkayastha, S. (1998) Simple proofs of two results on convolutions of unimodal distributions. Statistics & Probability Letters, 39, 97–100.
Rashid et al. (2011) Rashid, N. U., Giresi, P. G., Ibrahim, J. G., Sun, W. and Lieb, J. D. (2011) ZINBA integrates local covariates with DNA-seq data to identify broad and narrow regions of enrichment, even within amplified genomic regions. Genome Biology, 12, R67.
Rigaill (2015) Rigaill, G. (2015) A pruned dynamic programming algorithm to recover the best segmentations with 1 to ${K}_{max}$ change-points. Journal de la Société Française de Statistique, 156, 180–205.
Shewhart (1930) Shewhart, W. A. (1930) Economic quality control of manufactured product. Bell System Technical Journal, 9, 364–389.
Texier et al. (2016) Texier, G., Farouh, M., Pellegrin, L., Jackson, M. L., Meynard, J.-B., Deparis, X. and Chaudet, H. (2016) Outbreak definition by change point analysis: a tool for public health decision? BMC Medical Informatics and Decision Making, 16.
Truong et al. (2020) Truong, C., Oudre, L. and Vayatis, N. (2020) Selective review of offline change point detection methods. Signal Processing, 167, 107299.
Vaisman et al. (2010) Vaisman, L., Zariffa, J. and Popovic, M. R. (2010) Application of singular spectrum-based change-point analysis to EMG-onset detection. Journal of Electromyography and Kinesiology, 20, 750–760.
Yuan et al. (2020) Yuan, J., Li, M., Lv, G. and Lu, Z. K. (2020) Monitoring transmissibility and mortality of COVID-19 in europe. International Journal of Infectious Diseases, 95, 311–315.
Zhang et al. (2008) Zhang, Y., Liu, T., Meyer, C. A., Eeckhoute, J., Johnson, D. S., Bernstein, B. E., Nussbaum, C., Myers, R. M., Brown, M., Li, W. and Liu, X. S. (2008) Model-based analysis of ChIP-seq (MACS). Genome Biology, 9, R137.
Zhao and Yau (2019) Zhao, Z. and Yau, C. Y. (2019) Alternating pruned dynamic programming for multiple epidemic change-point estimation. Preprint arXiv:1907.06810.
Zheng et al. (2019) Zheng, C., Eckley, I. A. and Fearnhead, P. (2019) Consistency of a range of penalised cost approaches for detecting multiple changepoints. Preprint arXiv:1911.01716.
Zidek and van Eeden (2003) Zidek, J. V. and van Eeden, C. (2003) Uncertainty, entropy, variance and the effect of partial information, vol. 42 of Lecture Notes–Monograph Series, 155–167. Institute of Mathematical Statistics.

Appendix A Proof of Theorem 1 (Convergence of Algorithm 1)

Bottou (1998) analyses the case of an online algorithm iteratively minimising some function $f(x,w)$ (where $x$ represents the complete data and $w$ the parameters). Data points $\{x_{t}\}$ arrive sequentially, and at each iteration an estimate $w_{t}$ of the location of the minimum $w^{*}$ is obtained using some update function $H(x,w)$ and learning rate $\gamma_{t}$ as:

w_{t+1}=w_{t}-\gamma_{t}H(x_{t+1},w_{t}).

(12)

This updating mechanism gives rise to stochastic gradient descent if $\mathbb{E}H(x_{t+1},w_{t})=\nabla_{w}f(x,w)$ , but for the following argument this is not required.

To make the link with Algorithm 1 explicit, the update equation applied by this algorithm can be written as:

w_{t+1}=w_{t}+\gamma_{t}(x_{t+1}-w_{t}).

Then $w^{*}=\theta_{0}$ (i.e., the background mean that is to be estimated), and we ask whether the sequence of updates converges $w_{t}\rightarrow w^{*}$ . It was shown by Bottou (1998) that this occurs almost surely if the following three conditions are met:

1.

“convexity” – a single optimum $w^{*}$ exists and the expected value of the updates always points towards it:

$\forall\epsilon>0,\inf_{(w-w^{*})^{2}>\epsilon}(w-w^{*})\mathbb{E}H(x,w)>0;$ (13)
2.

learning rate convergence:

$\sum_{i=1}^{\infty}\gamma_{t}=\infty\text{ , }\sum_{i=1}^{\infty}\gamma^{2}_{t}<\infty;$ (14)
3.

bounded variance of the updates:

$\mathbb{E}H(x,w)^{2}\leq A+B(w-w^{*})^{2}\text{ , }A,B\geq 0.$ (15)

Thus, proof of convergence of our algorithm reduces to showing that these requirements are satisfied. We start with the assumption that the global mean of segment points is also $\theta_{0}$ , and then relax this requirement.

The following lemma will be needed:

Lemma 1.

Let $f$ be a unimodal distribution, symmetric around a point $\mu$ (so that $f(x_{1})<f(x_{2})$ when $x_{1}<x_{2}\leq\mu$ and $f(x_{1})>f(x_{2})$ when $\mu\leq x_{1}<x_{2}$ ), such as a Gaussian. Consider a truncated random variable $X$ with pdf:

g(x)=\begin{cases}0&\text{if }x<m-a\\ \frac{f(x)}{P(m-a\leq x\leq m+a)}&\text{if }m-a\leq x\leq m+a\\ 0&\text{if }x>m+a\end{cases}

for some $a>0,m$ . Then $\inf_{m\neq\mu}(\mathbb{E}X-m)(\mu-m)>0$ .

Proof.

$\displaystyle\mathbb{E}X-m$	$\displaystyle=\int_{m-a}^{m+a}(x-m)f(x)dx$
	$\displaystyle=\int_{-a}^{0}yf(y+m)dy+\int_{0}^{a}yf(y+m)dy$
	$\displaystyle=\int_{0}^{a}y(f(m+y)-f(m-y))dy.$	(16)

When $m+a<\mu$ , $f$ is increasing throughout the integration range, and $\mathbb{E}X-m>0$ ; the opposite is true for $m-a>\mu$ . If $m-a<m<\mu<m+a$ , split the integral in (A) as:

\mathbb{E}X-m=\int_{0}^{\mu-m}y(f(m+y)-f(m-y))dy+\int_{\mu-m}^{a}y(f(m+y)-f(m-y))dy.

The first integral covers the range where $f$ is increasing, and thus is positive. Since $\mu-m>0$ , $|m+y-\mu|<|m-y-\mu|$ for $y>0$ , and $f(m+y)>f(m-y)$ by symmetry of $f$ around $\mu$ and monotonicity, so the second interval is positive as well. Similarly, $\mathbb{E}X-m<0$ for $m-a<\mu<m<m+a$ . ∎

A.1 When the global mean of segments matches the background mean

Consider the case that the background points are independent draws from $\mathcal{N}(w^{*},\sigma^{2})$ , so that the points within each segment are $\mathcal{N}(\theta_{i},\sigma^{2})$ , with $\sigma^{2}$ known, and $\theta\sim\mathcal{N}(w^{*},\tau^{2})$ . Let $w_{t}$ be the value of the background mean estimated by Algorithm 1 after processing $t$ data points. In this case $w_{t}\xrightarrow{a.s.}w^{*}$ .

Proof.

Denote the true class of the next data point $x_{t+1}$ by $\delta_{t+1}$ (1 for background points, 0 for signal). Algorithm 1 estimates this as:

\hat{\delta}_{t+1}=\begin{cases}1&\text{if }F(t)+C^{0}(x_{t+1};w_{t})<\min_{1\leq k\leq l}F(t+1-k)+C(x_{t+2-k:t+1})+\beta\\ 0&\text{otherwise}.\end{cases}

Initially, assume for simplicity that the true maximum segment length is 1 (and so only $k=1$ is tested). When $\hat{\delta}_{t+1}=1$ , the background estimate is updated as:

w_{t+1}=w_{t}+\frac{1}{\sum_{i=1}^{t}\hat{\delta}_{i}+1}(x_{t+1}-w_{t})

(otherwise $w_{t+1}=w_{t}$ ). So $\gamma_{t}=1/(\sum_{i=1}^{t}\hat{\delta}_{i}+1)$ , and hence the learning rate convergence conditions (14) are satisfied.

Substituting in the costs based on a one-dimensional Gaussian pdf $\phi$ , $\hat{\delta}_{t+1}=1$ if:

	$\displaystyle-\log\phi(x_{t+1};w_{t},\sigma^{2})<-\log\phi(x_{t+1};x_{t+1},\sigma^{2})+\beta$
$\displaystyle\Rightarrow$	$\displaystyle\frac{1}{\sigma\sqrt{2\pi}}\exp{\frac{-(x_{t+1}-w_{t})^{2}}{2\sigma^{2}}}>\frac{1}{\sigma\sqrt{2\pi}}e^{-\beta}$
$\displaystyle\Rightarrow$	$\displaystyle\|x_{t+1}-w_{t}\|<\sqrt{2\beta\sigma^{2}}$
$\displaystyle\Rightarrow$	$\displaystyle x_{t+1}\in(w_{t}-\sqrt{2\beta\sigma^{2}};w_{t}+\sqrt{2\beta\sigma^{2}}).$	(17)

The distribution of true segment data points $f(x_{t+1}|\delta_{t+1}=0)$ is symmetric with mean $w^{*}$ by assumption (as is easily verified for the Gaussian case). The background point distribution is also automatically symmetric. Thus, the overall distribution of the points used to update the $w_{t}$ estimate is a truncation of a symmetric unimodal distribution. In the present case, it is a truncated normal with limits $(w_{t}-\sigma\sqrt{2\beta},w_{t}+\sigma\sqrt{2\beta})$ , based on (17); more generally it is a truncated variant of the parent distribution with symmetric limits of the form $w_{t}\pm a$ , and parent mean $w^{*}$ (that the acceptance set is an interval follows from the unimodality of $f$ ).

This means that $f(x_{t+1}|\hat{\delta}_{t+1}=1)$ satisfies the requirements for Lemma 1 with $\mu=w^{*}$ , which implies the “convexity” condition (13):

\inf_{(w-w^{*})^{2}>\epsilon}(w-w^{*})\mathbb{E}H(x_{t+1},w)=\inf_{(w-w^{*})^{2}>\epsilon}(w-w^{*})(w-\mathbb{E}(x_{t+1}|\hat{\delta}_{t+1}=1))>0

Following a similar approach – conditioning on $\delta_{t+1}$ – and using the law of total variance it can be shown that the variance of $\mathbb{E}(w-x_{t+1})$ is finite, as required for the condition (15), and so $w_{t}\xrightarrow{a.s.}w^{*}$ .

Remark. So far, we assumed that segment length $k=1$ . If segments occur and are tested in non-overlapping windows of any fixed size $k\geq 2$ , the result is similar: $\forall j\in[t-k+1;t],\hat{\delta}_{j+1}=1$ if:

-\sum_{i=t-k+1}^{i=t}\log\phi(x_{i+1};w_{t-k},\sigma^{2})<-\sum_{i=t-k+1}^{i=t}\log\phi(x_{i+1};\bar{x}_{k},\sigma^{2})+\beta,

where $\bar{x}_{k}=\sum_{i=t-k+1}^{t}x_{i+1}/k$ . This can be expressed as truncation limits for accepted $\bar{x}_{k}$ , analogously to (17):

$\displaystyle\beta$	$\displaystyle>-\sum_{i=t-k+1}^{t}\log\phi(x_{i+1};w_{t-k},\sigma^{2})+\sum_{i=t-k+1}^{i=t}\log\phi(x_{i+1};\bar{x}_{k},\sigma^{2})$
	$\displaystyle=\sum_{i=t-k+1}^{t}\frac{(x_{i+1}-w_{t-k})^{2}}{2\sigma^{2}}-\frac{(x_{i+1}-\bar{x}_{k})^{2}}{2\sigma^{2}}$
	$\displaystyle=\frac{1}{2\sigma^{2}}\sum_{i=t-k+1}^{t}(w_{t-k}-\bar{x}_{k})^{2}$
	$\displaystyle\Rightarrow\sqrt{2\beta\sigma^{2}/k}>\|w_{t-k}-\bar{x}_{k}\|.$	(18)

The pdf of $\bar{x}_{k}$ is a $k$ -fold convolution of $f(x)$ . Since $f$ , as shown earlier, is symmetric unimodal, so is their convolution, and hence the distribution of $\bar{x}_{k}$ , with a mean $\mu=w^{*}$ (Purkayastha, 1998). In the special case when $f$ is the normal pdf, this can also be shown directly from Gaussian properties. Then Lemma 1 implies condition (13), and the rest of the proof follows as before.

When all segment positions (overlapping or not) are tested, the background point acceptance rule is:

\hat{\delta}_{t}=1\text{ if }x_{t}\in\bigcap_{1\leq k\leq l}S_{k},\text{ with }S_{k}=\{x_{t}:F(t-1)+C^{0}(x_{t})<F(t-k)+C(x_{t-k+1:t})+\beta\}.

As demonstrated earlier, $S_{1}=(w_{t}-\sqrt{2\beta\sigma^{2}},w_{t}+\sqrt{2\beta\sigma^{2}})$ . Define three sets of $x_{t}$ based on which rules they pass: $X_{1}=x_{t}\in S_{1}$ , $X_{a}=x_{t}\in S_{1}\cap S_{\geq 2}$ , $X_{r}=x_{t}\in S_{1}\setminus S_{\geq 2}$ , and let $P_{a}$ and $P_{r}$ be the probabilities of the corresponding $x_{t}$ sets. Clearly, $X_{1}=X_{a}\cup X_{r}$ . We are interested in the mean of the points accepted as background, i.e. $\mathbb{E}X_{a}$ . Assume w.l.o.g. $\mu=0,\sigma=1,w_{t}>\mu$ , as the other case is symmetric. We will now show that for sufficiently large $n$ , $\mathbb{E}X_{a}<w_{t}$ , satisfying (13).

Using the conditional mean formula:

	$\displaystyle\mathbb{E}X_{1}$	$\displaystyle=P_{r}\mathbb{E}X_{r}+P_{a}\mathbb{E}X_{a}$
	$\displaystyle\mathbb{E}X_{a}$	$\displaystyle=\mathbb{E}X_{1}/(1-P_{r})-P_{r}\mathbb{E}X_{r}/(1-P_{r}).$		(19)

Assume for now $\mathbb{E}X_{r}=\mu=0$ . Then, to obtain $\mathbb{E}X_{a}<w_{t}$ , we need:

\mathbb{E}X_{1}/w_{t}<1-P_{r}.

(20)

Denote $p=\sqrt{2\beta}$ , which is an increasing function of $n$ , and consider the growth of both sides of (20) as $n$ increases. For the Gaussian or other distributions in the exponential family with mean $\mu$ , truncated to a symmetric region $(-a,a)$ , it is known that $Var(X|S;\mu)=\frac{d}{d\mu}\mathbb{E}(X|S;\mu)$ (Zidek and van Eeden, 2003). Then (denoting $S^{\prime}_{1}=(-p;+p)$ ):

	$\displaystyle\mathbb{E}X_{1}$	$\displaystyle=w_{t}+\mathbb{E}(X\|S^{\prime}_{1};-w_{t})$
		$\displaystyle=w_{t}+\mathbb{E}(X\|S^{\prime}_{1};0)+\int_{0}^{w_{t}}\frac{d\mathbb{E}(X\|S^{\prime}_{1};-a)}{da}da$
		$\displaystyle=w_{t}+0-\int_{0}^{w_{t}}\frac{d\mathbb{E}(X\|S^{\prime}_{1};a)}{da}da$
		$\displaystyle<w_{t}-w_{t}\min_{0\leq a\leq w_{t}}Var(X\|S^{\prime}_{1};a).$

Hence:

\mathbb{E}X_{1}/w_{t}<1-Var(X|S^{\prime}_{1};w_{t})=1-\int_{w_{t}-p}^{w_{t}+p}x^{2}f(x)dx,

where $f(x)$ is the pdf of $x_{t}$ given $x\in S_{1}$ . Hence, this side grows with $p$ as $-x^{2}f(x)$ .

To analyse $P_{r}$ , we first simplify the background acceptance condition. For any $k\geq 2$ , by definition of $F$ and additivity of $C$ , we have:

	$\displaystyle F(t-k)+C(x_{t-k+1:t})+\beta$	$\displaystyle=F(t-k)+C(x_{t-k+1:t-1};\bar{x}_{k})+C(x_{t};\bar{x}_{k})+\beta$
		$\displaystyle=F(t-k)+C(x_{t-k+1:t-1};\bar{x}_{-t})+(k-1)d(\bar{x}_{k},\bar{x}_{-t})+d(\bar{x}_{k},x_{t})+\beta$
		$\displaystyle=F(t-1)+A(t-1)+(k-1)d(\bar{x}_{k},\bar{x}_{-t})+d(\bar{x}_{k},x_{t}),$

where $d$ is some distance function, $\bar{x}_{-t}$ is the mean of points $x_{t-k+1:t-1}$ , and $A(t)\geq 0$ is a constant depending only on $x_{1:t}$ .

It is also helpful to note that $F(t-1)\leq F(t-k)+C^{0}(x_{t-k+1:t-1})$ , hence:

$\displaystyle A(t-1)$	$\displaystyle=F(t-k)+C(x_{t-k+1:t-1};\bar{x}_{-t})+\beta-F(t-1)$
	$\displaystyle\geq C(x_{t-k+1:t-1};\bar{x}_{-t})+\beta-C^{0}(x_{t-k+1:t-1})$
	$\displaystyle\geq\beta-(k-1)d(\bar{x}_{-t},w_{t-1}).$	(21)

This corresponds to the case when all $x_{t-k+1:t-1}$ were identified as background.

Using the Gaussian cost, i.e. $d(a,b)=(a-b)^{2}/2$ , and recursive formula for the mean, the acceptance condition for $k$ becomes:

	$\displaystyle F(t-1)+C^{0}(x_{t})$	$\displaystyle<F(t-1)+A(t-1)+\frac{k-1}{2k^{2}}(x_{t}-\bar{x}_{-t})^{2}+\frac{(k-1)^{2}}{2k^{2}}(x_{t}-\bar{x}_{-t})^{2}$
	$\displaystyle\Rightarrow(x_{t}-w_{t-1})^{2}$	$\displaystyle<2A(t-1)+\frac{k-1}{k}(x_{t}-\bar{x}_{-t})^{2}.$		(22)

By substituting in the value of $A(t-1)$ from (A.1), we obtain the following lower bound for $P(x\in S_{k}|x)$ :

	$\displaystyle P(x\in S_{k}\|x)$	$\displaystyle\geq P\left((x_{t}-w_{t-1})^{2}<2\beta-(k-1)(\bar{x}_{-t}-w_{t-1})^{2}+\frac{k-1}{k}(x_{t}-\bar{x}_{-t})^{2}\right)$
		$\displaystyle=P\left((x_{t}-w_{t-1})^{2}-\frac{k-1}{k}(x_{t}-\bar{x}_{-t})^{2}+(k-1)(\bar{x}_{-t}-w_{t-1})^{2}<2\beta\right)$
		$\displaystyle\geq P\left((x_{t}-w_{t-1})^{2}+(k-1)(\bar{x}_{-t}-w_{t-1})^{2}<2\beta\right)$
		$\displaystyle\geq 1-\mathbb{E}\left((x_{t}-w_{t-1})^{2})+(k-1)(\bar{x}_{-t}-w_{t-1})^{2}\right)/\left(2\beta\right)$
	$\displaystyle P(x\notin S_{k}\|x)$	$\displaystyle\leq\mathbb{E}\left((x_{t}-w_{t-1})^{2})+(k-1)(\bar{x}_{-t}-w_{t-1})^{2}\right)/\left(2\beta\right).$

Thus, we have the following bound for $P_{r}$ at any $k$ :

	$\displaystyle P_{r}$	$\displaystyle=\int_{w_{t}-p}^{w_{t}+p}P(x\notin S_{k}\|x)f(x)dx$
		$\displaystyle\leq\int_{w_{t}-p}^{w_{t}+p}O(x^{2}/p^{2})f(x)dx.$		(23)

Therefore, as $N$ increases, $1-P_{r}$ grows faster than $1-\mathbb{E}X_{1}/w_{t}=1-\int_{w_{t}-p}^{w_{t}+p}x^{2}f(x)dx$ . This means that $\exists p_{0}$ , and thus $\exists n_{0}$ , such that for $n>n_{0}$ , and thus $p>p_{0}$ , (20) holds.

So far we assumed $\mathbb{E}X_{r}=\mu$ . Clearly, for larger values of $\mathbb{E}X_{r}$ , $\mathbb{E}X_{a}$ is even smaller and (13) is satisfied.

For the case when $\mathbb{E}X_{r}<\mu$ , consider the worst case scenario $\mathbb{E}X_{r}=w_{t}-p$ (this is the bound to $X_{1}$ , and thus to $X_{r}$ , imposed by $S_{1}$ ). Similarly to (19), we need:

		$\displaystyle\mathbb{E}X_{a}<w_{t}$
	$\displaystyle\Rightarrow$	$\displaystyle\mathbb{E}X_{1}/(1-P_{r})-P_{r}(w_{t}-p)/(1-P_{r})<w_{t}$
	$\displaystyle\Rightarrow$	$\displaystyle\mathbb{E}X_{1}/w_{t}-P_{r}+P_{r}p/w_{t}<1-P_{r}$
	$\displaystyle\Rightarrow$	$\displaystyle P_{r}p/w_{t}<1-\mathbb{E}X_{1}/w_{t}$
	$\displaystyle\Rightarrow$	$\displaystyle P_{r}p/w_{t}<\int_{w_{t}-p}^{w_{t}+p}x^{2}f(x)dx.$

However, based on (23), $P_{r}p/w_{t}\leq\int_{w_{t}-p}^{w_{t}+p}O(x^{2}/p)f(x)dx$ , so again $\exists n>n_{0}$ such that $\mathbb{E}X_{a}<w_{t}$ , and condition (13) holds.

And so overall $\mathbb{E}(x_{t}|\hat{\delta}_{t}=1)$ satisfies condition (13). Since the distribution of accepted $x_{t}$ still has bounded support imposed by $S_{1}$ , condition (15) still holds, and the learning rate condition (14) holds as before, implying convergence.

∎

A.2 When the global mean of segments does not match the background mean

Consider now $\theta\sim\mathcal{N}(\mu,\tau^{2})$ for some $\mu\neq w^{*}$ , in particular $\mu>w^{*}$ , so that the overall mean of segment points is $\mathbb{E}(\theta)>w^{*}$ . Then if $w^{*}<w_{t}<\mathbb{E}(\theta)$ , any segment points that were misclassified as background will (on average) push the estimates away from the background mean, in violation of the “convexity” condition (13).

We assume that each segment point is followed by no less than $n$ background points. Then, as $n\rightarrow\infty$ , $w_{t}\xrightarrow{a.s}w^{*}$ . For every finite $n$ , $\exists\epsilon>0$ such that $P(|w_{t}-w^{*}|>\epsilon)=0$ .

Proof.

Suppose that a misclassification at time $T$ is followed by $n$ correctly classified background points: $\delta_{T}=0$ , $\delta_{t}=1$ for $t\in[T+1;T+n]$ , $\hat{\delta}_{t}=1$ for $t\in[T;T+n]$ . For the points $t\in[T+1;T+n]$ , almost sure convergence of $w_{t}$ was established above, i.e. for all $\epsilon>0$ , there exists a $t_{0}$ such that $\forall t:n\geq t\geq t_{0},P(|w_{T+t}-w^{*}|<\epsilon)=1$ . Therefore, given $n\geq t_{0}$ :

	$\displaystyle P(\|w_{T+n}-w^{}\|<\|w_{T-1}-w^{}\|)=1$
$\displaystyle\Rightarrow$	$\displaystyle\begin{cases}P(w_{T+n}-w_{T-1}<0)=1,\text{ if }w_{T-1}-w^{}>0\\ P(w_{T+n}-w_{T-1}>0)=1,\text{ if }w_{T-1}-w^{}<0\end{cases}$
$\displaystyle\Rightarrow$	$\displaystyle\inf_{w_{T-1}\neq w^{}}(w_{T-1}-w^{})\mathbb{E}(w_{T+n}-w_{T-1})<0.$	(24)

Indexing the segment-background cycles by $i$ , denote the first estimate of that segment by $w^{\prime}_{i}$ , so the set of these estimates are:

\{w^{\prime}_{i}\}=\{w_{1},\dots,w_{T-2-n},w_{T-1},w_{T+n},w_{T+1+2n},\dots\}.

The elements of this sequence can be expressed recursively as:

w^{\prime}_{i+1}=w^{\prime}_{i}-\gamma^{\prime}_{i}H^{\prime}(\{x^{\prime}_{i}\},w^{\prime}_{i}),

with $\{x^{\prime}_{i}\}=\{x_{t}:i(n+1)\leq t\leq i(n+1)+n\}$ .

From (A.2), $\mathbb{E}(w^{\prime}_{i}-w^{\prime}_{i+1})=\gamma^{\prime}_{i}\mathbb{E}H^{\prime}(\{x^{\prime}_{i}\},w^{\prime}_{i})$ is “convex” as defined in (13), and because $\gamma^{\prime}_{i}>0$ so is $\mathbb{E}H^{\prime}(\{x^{\prime}_{i}\},w^{\prime}_{i})$ .

Let $\gamma^{\prime}_{i}=\frac{1}{i(n+1)+1}$ . Then:

	$\displaystyle H^{\prime}(\{x^{\prime}_{i}\},w^{\prime}_{i})$	$\displaystyle=\sum_{t=i(n+1)}^{i(n+1)+n}\gamma_{t}(w_{t-1}-x_{t})/\gamma^{\prime}_{i}$
		$\displaystyle=\sum_{t=i(n+1)}^{i(n+1)+n}\frac{i(n+1)+1}{t+1}(w_{t-1}-x_{t})$
		$\displaystyle<(n+1)(w_{t-1}-x_{t}).$

So for $n<\infty$ , conditions (14)–(15) are satisfied as well, and $w^{\prime}_{i}\xrightarrow{a.s.}w^{*}$ . (When $n\to\infty$ , the convergence conditions are satisfied directly without using the sequence $\{w^{\prime}_{i}\}$ .) ∎

A.3 Martingale Approach

We can also describe the update process over the background points using martingales. The algorithm estimates are random variables $w_{t}$ ; let $\{\mathcal{W}_{t}\}$ be the sequence of $\sigma$ -algebras such that for each $t$ , $w_{t}$ is measurable with respect to $\mathcal{W}_{t}$ . Using Lemma 1, and assuming $w^{*}<w_{t}$ again, within each cycle the estimates comprise a supermartingale $\mathbb{E}(w_{t+1}|\mathcal{W}_{t})<w_{t}$ over the points $T\leq t<\min(T+n,FHT_{w}(w^{*}))$ , here $FHT_{x}(a)=\inf\{t:x_{t}\leq a\}$ is the first hitting time of the process realisation $\{x_{t}\}$ to value $a$ .

Consider again the problematic case when the global mean does not match the background mean and misclassification pushes the estimate away from the background mean, i.e. $w^{*}<w_{T-1}<w_{T}<\mathbb{E}(\theta)$ . In order for $w^{\prime}_{i}$ to converge, we need the perturbed estimates to return to a value below $w_{T-1}$ in each cycle. At the extremes, we have:

	$\displaystyle FHT_{w}(w_{T})=T$	starting position
	$\displaystyle FHT_{w}(w_{T-1})<\infty$	$\displaystyle\text{for sufficiently large $n$, because }w_{t}\xrightarrow{a.s.}w^{*}.$

Clearly, $FHT_{w}(w_{T-1})\leq FHT_{w}(w^{*})$ . However, the number of background points $n$ required to satisfy $FHT_{w}(w_{T-1})<T+n$ will depend on five factors: the distribution $f_{B}$ and penalty $p$ (since they determine the distribution of update values $H$ ), the size of estimated background set at time $T$ (as it determines the relevant $\gamma_{t}$ ), and $w_{T-1}$ and $w_{T}$ .

In practice, $n$ is bounded by the available data, so there is a non-zero probability that, over the segment-background cycles indexed by $i$ :

\max_{i}FHT_{w}(w_{T_{i}-1})>T_{i}+n.

In that case, define $b=\min\{a:\forall i,FHT_{w}(a)\leq T_{i}+n\}$ ; the final estimate of $w^{\prime}_{i}$ will be bound by $[w^{*},b]$ . As $n$ increases, $P(|w^{*}-b|>\epsilon)\rightarrow 0$ for any $\epsilon>0$ .

Similar reasoning applies when $\mathbb{E}(\theta)<w_{T}<w_{T-1}<w^{*}$ .

Appendix B Proof of Theorem 2 (Consistency)

Proof.

Some general consistency results for changepoint detection by penalised cost methods are given in Fisch et al. (2018). In particular, an equivalent of our Theorem 2 is established for any algorithm that detects changepoints by exactly minimising a cost $F(n;\{s_{i},e_{i}\},\theta,\hat{\mu},\hat{\sigma})$ , where $\hat{\cdot}$ marks robust estimates of background parameters. While the original statement uses the median and interquartile range of $x_{0:n}$ for $\hat{\mu}$ and $\hat{\sigma}$ , the proof only requires that the estimates satisfy certain upper bounds on deviations from the true values. Therefore, we will first show that the online estimates produced by Algorithm 1 are within these bounds, and then follow the rest of the proof from Fisch et al. (2018).

Noting again that Algorithm 1 is effectively a stochastic gradient descent procedure, with each data point seen precisely once, we can use the error bound on estimates produced by such algorithms as provided in Theorem 7.5 of Harvey et al. (2019):

Theorem 3 ((Harvey et al., 2019)).

Let function $f(w)$ be 1-strongly convex and 1-Lipschitz. A stochastic gradient algorithm for minimising this function runs for $T$ cycles, and at each cycle updates the estimate as in (12) with $\gamma_{t}=1/t$ , $\mathbb{E}H=\nabla f(w)$ . Then:

P\left(\lVert w_{t}-w^{*}\rVert^{2}\leq O\left(\frac{\log(1/\delta)}{t}\right)\right)\geq 1-\delta.

Using $\delta=n^{-\epsilon}$ , and assuming without loss of generality that $\sigma_{0}=1$ , we can establish an upper bound on the error of background parameters estimated by Algorithm 1 after $n$ cycles:

P\left((\hat{\mu}-\mu_{0})^{2}\leq O\left(\frac{\log(n^{\epsilon})}{n}\right)\right)=P\left(|\hat{\mu}-\mu_{0}|\leq O\left(\sqrt{\epsilon}\sqrt{\frac{\log n}{n}}\right)\right)\geq 1-n^{-\epsilon}

P\left(|\hat{\sigma}^{2}-\sigma^{2}_{0}|\leq O\left(\sqrt{\epsilon}\sqrt{\frac{\log n}{n}}\right)\right)\geq 1-n^{-\epsilon}.

Application of Boole’s inequality leads to:

P\left(|\hat{\mu}-\mu_{0}|\leq D_{1}\sigma_{0}\sqrt{\frac{\log(n)}{n}},\left|\frac{\hat{\sigma}^{2}}{\sigma^{2}_{0}}-1\right|\leq D_{2}\sqrt{\frac{\log(n)}{n}}\right)\geq 1-C_{1}n^{-\epsilon},

(25)

for some constants $C_{1},D_{1},D_{2}$ and sufficiently large $n$ . Since the objective function $f$ in our algorithm is the Gaussian log-likelihood (i.e., the updates $\mathbb{E}H$ approximate its gradient), for any given segmentation it is 1-strongly convex. For other functions, overall consistency can still be achieved similarly, but the convergence rate may be slower than $n^{-\epsilon}$ .

Having established the bound on estimate errors, we can use Lemma 9 from Fisch et al. (2018) and the proof method reported there.

First, introduce an event $E$ based on a combination of bounds limiting the behaviour of Gaussian data $x_{1:n}$ , which for any $\epsilon>0$ , occurs with probability $P(E)>1-C_{2}n^{-\epsilon}$ , with some constant $C_{2}$ and sufficiently large $n$ (Lemmas 1 and 2 in Fisch et al. (2018)). Conditional on this event, the following lemma holds for the epidemic cost $F$ defined as in (5):

Lemma 2 ((Fisch et al., 2018)).

Let $\{\tau\}$ be the set of true segment positions $\{(s_{i},e_{i})\}$ , and $\theta$ the vector of true segment means. Assume $E$ holds, and some $\hat{\mu},\hat{\sigma}$ are available for which the event in (25) holds. Then, there exist constants $C_{3}$ and $n_{1}$ such that when $n>n_{1}$ ,

F(n;\{\tau\},\theta,\hat{\mu},\hat{\sigma})-F(n;\{\tau\},\theta,\mu_{0},\sigma_{0})<C_{4}\log n.

This lemma, together with results established for classical changepoint detection, can be used to show that the cost of any inconsistent solution will exceed the cost based on true segment positions and parameters (Proposition 8 in Fisch et al. (2018)):

Proposition 3 ((Fisch et al., 2018)).

Define $\{\tau^{\prime}\}$ to be any set of segments $\{(s_{i},e_{i})\}$ that does not satisfy the consistency event in (7). Let $\tilde{\theta}=\operatorname*{argmin}_{\theta}F(n;\theta)$ be the parameters estimated by minimising the cost for a given segmentation (i.e. the vector of means and/or variances of $x_{s_{i}:e_{i}}$ for each $i$ ). Assume $E$ holds. Then there exist constants $C_{4}$ and $n_{2}$ such that, when $n>n_{2}$ :

F(n;\{\tau^{\prime}\},\tilde{\theta},\hat{\mu},\hat{\sigma})\geq F(n;\{\tau\},\theta,\hat{\mu},\hat{\sigma})+C_{3}\log(n)^{1+\delta/2}.

See the original publication for a detailed proof of these results.

Finally, for a given set of changepoints, using fitted maximum-likelihood parameters by definition results in minimal cost:

F(n;\{\tau\},\theta,\hat{\mu},\hat{\sigma})\geq F(n;\{\tau\},\tilde{\theta},\hat{\mu},\hat{\sigma}).

Thus, when Proposition 3 holds, we have:

F(n;\{\tau^{\prime}\},\tilde{\theta},\hat{\mu},\hat{\sigma})>F(n;\{\tau\},\tilde{\theta},\hat{\mu},\hat{\sigma}),

and an exact minimisation algorithm will always find a solution in the consistent set. The overall probability of the events required for Proposition 3 is a combination of $P(E)$ , established before, and (25), which by Boole’s inequality is:

P>1-C_{5}n^{-\epsilon},

for any $\epsilon>0$ , $n>n_{3}$ and some constants $n_{3},C_{5}$ . ∎

Appendix C Proof of Proposition 2 (Pruning of Algorithm 2)

Proof.

Denote the true start and end of a nuisance segment as $s_{j},e_{j}$ . Consider the case when $s_{j}\in(t-A(n);t]$ . Pruning at time $t$ will not remove this point (i.e. $s_{j}\notin\bm{k}_{pr,t}$ ) iff:

C^{0}(x_{t-A(n):s_{j}-1})+C^{N}(x_{s_{j}:t})<C^{0}(x_{t-A(n):m-1})+C^{N}(x_{m:t})+\alpha\log(n)^{1+\delta}

with $m$ such that the right hand side is minimised and $m\neq s_{j}$ .

Denote by $C(x_{a:b};\hat{\mu},\hat{\sigma})$ the Gaussian cost calculated with MLE estimates of the parameters (i.e. mean and variance of $x_{a:b}$ ). Note that since $C^{0}(x_{a:b})=C(x_{a:b};\mu_{0},\sigma_{0})$ and $C^{N}(x_{a:b})=C(x_{a:b};\hat{\mu},\sigma_{N})$ , the required event to preserve $s_{j}$ can be stated as

C(x_{t-A(n):s_{j}-1};\mu_{0},\sigma_{0})+C(x_{s_{j}:t};\hat{\mu},\sigma_{N})-\\ C(x_{t-A(n):m-1};\mu_{0},\sigma_{0})-C(x_{m:t};\hat{\mu},\sigma_{N})<\alpha\log(n)^{1+\delta}

We can establish the probability of this using the following bound (Proposition 4 in Fisch et al. (2018)):

Lemma 3.

Let $x_{1:n}$ be piecewise-Gaussian data. Choose any subset $x_{i:j},1\leq i<j\leq n$ , with a true changepoint at $s$ , i.e., we have $x_{t}\sim\mathcal{N}(\mu_{1},\sigma_{1})$ for $t\in[i;s-1]$ , and $x_{t}\sim\mathcal{N}(\mu_{2},\sigma_{2})$ for $t\in[s;j]$ . Then, for any candidate changepoint $\tau$ and any $\epsilon>0$ , there exist constants $B,n_{0},K_{1}$ such that:

C(x_{i:s-1};\mu_{1},\sigma_{1})+C(x_{s:j};\mu_{2},\sigma_{2})-C(x_{i:\tau-1};\hat{\mu},\hat{\sigma})-C(x_{\tau:j};\hat{\mu},\hat{\sigma})\leq K_{1}\log(n)

is true for all $i,j$ with $P\geq 1-Bn^{-\epsilon}$ when $n>n_{0}$ .

Now take $i=s_{j}-A(n)+1,j=s_{j}+A(n)-1$ . Note that there is one and only one changepoint within $x_{i:j}$ because of the required distance between changepoints. Applying Lemma 3 to such $x_{i:j}$ states that, conditional on an event with probability $P\geq 1-Bn^{-\epsilon}$ , the following is true for all $t\in[s_{j},s_{j}+A(n))$ :

	$\displaystyle C(x_{t-A(n):s_{j}-1};\mu_{0},\sigma_{0})+C(x_{s_{j}:t};\hat{\mu},\sigma_{N})-C(x_{t-A(n):m-1};\mu_{0},\sigma_{0})-C(x_{m:t};\hat{\mu},\sigma_{N})$
	$\displaystyle\leq C(x_{t-A(n):s_{j}-1};\mu_{0},\sigma_{0})+C(x_{s_{j}:t};\mu_{N},\sigma_{N})-C(x_{t-A(n):m-1};\hat{\mu},\hat{\sigma})-C(x_{m:t};\hat{\mu},\hat{\sigma})$
	$\displaystyle\leq K_{1}\log(n)<\alpha\log(n)^{1+\delta},$

where we also used the fact that $\hat{\mu},\hat{\sigma}=\operatorname*{argmin}_{\mu_{,}\sigma}C(x;\mu,\sigma)$ .

Therefore, with the same probability, $s_{j}\notin\bigcup_{t=s_{j}}^{s_{j}+A(n)-1}\bm{k}_{pr,t}$ . Also, $s_{j}\notin\bigcup_{t=s_{j}+A(n)}^{n}\bm{k}_{pr,t}$ because then $s_{j}\leq t-A(n)$ and is not considered in the pruning scheme, and clearly $s_{j}\notin\bigcup_{t=1}^{s_{j}-1}\bm{k}_{pr,t}$ . The case for $e_{j}$ follows by symmetry, and obviously no true changepoint can be pruned out if $s_{j},e_{j}=\emptyset$ , so the overall probability of retaining a true changepoint remains at $P\geq 1-Bn^{-\epsilon}$ . ∎

Appendix D Supplementary Table

Table S1: Comparison of Algorithm 2 detections without pruning (Full) or with global pruning (Pruned). Of the 5000 simulation runs, all runs where the two options produced any difference in segments were identified. All detections from these runs were extracted and are shown here. Segment types are N – nuisance, S – signal.

			(Full)			(Pruned)
Scenario	$n$	Run number	Segm. type	Start	End	Segm. type	Start	End
1	30	88	N	7	21	S	7	11
	30	88	S	13	15	N	12	21
	30	88	S	16	20	S	16	20
	100	84	N	21	69	S	21	50
	100	84	S	51	68	S	51	69
2	30	378	N	10	18	S	10	12
	30	378	S	13	15	S	16	18
	30	393	N	2	12	N	3	12
	30	393	S	4	7	S	4	7
	30	393	S	16	18	S	16	18
	30	393	S	22	24	S	22	24
	240	445	N	55	95	N	55	95
	240	445	S	121	144	S	121	144
	240	445	N	170	218	S	170	192
	240	445	S	193	214	S	215	218

	$\displaystyle\mathbb{E}X_{1}$	$\displaystyle=w_{t}+\mathbb{E}(X\|S^{\prime}_{1};-w_{t})$
		$\displaystyle=w_{t}+\mathbb{E}(X\|S^{\prime}_{1};0)+\int_{0}^{w_{t}}\frac{d\mathbb{E}(X\|S^{\prime}_{1};-a)}{da}da$
		$\displaystyle=w_{t}+0-\int_{0}^{w_{t}}\frac{d\mathbb{E}(X\|S^{\prime}_{1};a)}{da}da$
		$\displaystyle<w_{t}-w_{t}\min_{0\leq a\leq w_{t}}Var(X\|S^{\prime}_{1};a).$

Epidemic changepoint detection in the presence of nuisance changes

Abstract

1 Introduction

2 Background

3 Detection of changepoints with unknown background parameter

3.1 Convergence

Theorem 1.

3.2 Consistency

Theorem 2.

3.3 Pruning

Proposition 1.

Proof.

4 Detecting changepoints with a nuisance process

4.1 Problem setup

4.2 Proposed method

4.3 Pruning

Proposition 2.

5 Experiments

5.1 Simulations

5.1.1 Algorithm 1 estimates the background parameter consistently

5.1.2 Segment positions estimated by Algorithm 1 are accurate and consistent

5.1.3 Algorithm 2 recovers true signal segments under interference

5.2 Real-world Data

5.2.1 ChIP-seq

5.2.2 European mortality data

6 Discussion

7 Acknowledgements

References

Appendix A Proof of Theorem 1 (Convergence of Algorithm 1)

Lemma 1.

Proof.

A.1 When the global mean of segments matches the background mean

Proof.

A.2 When the global mean of segments does not match the background mean

Proof.

A.3 Martingale Approach

Appendix B Proof of Theorem 2 (Consistency)

Proof.

Theorem 3 ((Harvey et al., 2019)).

Lemma 2 ((Fisch et al., 2018)).

Proposition 3 ((Fisch et al., 2018)).

Appendix C Proof of Proposition 2 (Pruning of Algorithm 2)

Proof.

Lemma 3.

Appendix D Supplementary Table

Epidemic changepoint detection
in the presence of nuisance changes