# Call for Papers - The BabyLM Challenge: Sample-efficient pretraining
on a developmentally plausible corpus
https://babylm.github.io/

Alex Warstadt
ETH Zürich Leshem Choshen
IBM Research Aaron Mueller
Johns Hopkins University Ethan Wilcox
ETH Zürich Adina Williams
Meta AI Chengxu Zhuang
MIT

###### Abstract

We present the call for papers for the BabyLM Challenge: Sample-efficient
pretraining on a developmentally plausible corpus. This shared task is
intended for participants with an interest in small scale language modeling,
human language acquisition, low-resource NLP, and cognitive modeling. In
partnership with CoNLL and CMCL, we provide a platform for approaches to
pretraining with a limited-size corpus sourced from data inspired by the input
to children. The task has three tracks, two of which restrict the training
data to pre-released datasets of 10M and 100M words and are dedicated to
explorations of approaches such as architectural variations, self-supervised
objectives, or curriculum learning. The final track only restricts the amount
of text used, allowing innovation in the choice of the data, its domain, and
even its modality (i.e., data from sources other than text is welcome). We
will release a shared evaluation pipeline which scores models on a variety of
benchmarks and tasks, including targeted syntactic evaluations and natural
language understanding.

##  1 Motivation

Huge efforts have been put into optimizing LM pretraining at massive scales in
the last several years (Raffel et al., 2020; Brown et al., 2020; Chowdhery et
al., 2022; Hoffmann et al., 2022). While growing parameter counts often get
the most attention, datasets have also grown by orders of magnitude. These
increasingly larger pretraining datasets are visualized, to scale, in Figure
1. At the same time, there has been almost no progress in pretraining at
smaller human-like data scales.

Focusing on scaled-down pretraining has several potential benefits: First,
small-scale pretraining can be a sandbox for developing novel techniques that
improve data efficiency. These techniques have the potential to then scale up
to larger datasets commonly seen in applied NLP, and could be used to enhance
current approaches to modeling low-resource languages. Second, improving our
ability to train LMs on the same types and quantities of data that humans
learn from will give us greater access to more plausible cognitive models of
humans and help us understand what allows humans to acquire language so
efficiently (Keller, 2010; Dupoux, 2018). That is, even model failure can help
in developing hypotheses about the differences between human and LM language
learning.

The goal of this shared task will be to incentivize researchers with an
interest in pretraining and/or cognitive modeling to focus their efforts on
optimizing pretraining given data limitations inspired by human development.
Additionally, we hope to democratize research on pretraining—which is
typically thought to be practical only for large industry groups—by drawing
attention to open problems that can be addressed on a university budget.

Figure 1: Data Scale: Modern Language Models are trained on data multiple
orders of magnitude larger than the amount available to a typical human child.
Image based off Fig. 1 from Warstadt and Bowman (2022)

##  2 Key Dates

HTML]E9F0E9  • January 2023: Training data released • March 2023: Evaluation
pipeline released • July 15, 2023: Results due • August 1, 2023: Paper
submissions due • Date TBA: Presentation at CoNLL

|  | # Words |
---|---|---|---
Dataset | Domain | Strict-small | Strict | Proportion
CHILDES (MacWhinney, 2000) | Child-directed speech | 0.44M | 4.21M | 5%
British National Corpus (BNC),1 dialogue portion | Dialogue | 0.86M | 8.16M | 8%
Children’s Book Test (Hill et al., 2016) | Children’s books | 0.57M | 5.55M | 6%
Children’s Stories Text Corpus2 | Children’s books | 0.34M | 3.22M | 3%
Standardized Project Gutenberg Corpus (Gerlach and Font-Clos, 2018) | Written English | 0.99M | 9.46M | 10%
OpenSubtitles (Lison and Tiedemann, 2016) | Movie subtitles | 3.09M | 31.28M | 31%
QCRI Educational Domain Corpus (QED; Abdelali et al., 2014) | Educational video subtitles | 1.04M | 10.24M | 11%
Wikipedia3 | Wikipedia (English) | 0.99M | 10.08M | 10%
Simple Wikipedia4 | Wikipedia (Simple English) | 1.52M | 14.66M | 15%
Switchboard Dialog Act Corpus (Stolcke et al., 2000) | Dialogue | 0.12M | 1.18M | 1%
_Total_ | – | 9.96M | 98.04M | 100%

Table 1: The datasets we release for the Strict and Strict-small tracks of the
BabyLM Challenge. We present the number of words in the training set of each
corpus that we include. 1http://www.natcorp.ox.ac.uk
2https://www.kaggle.com/datasets/edenbd/children-stories-text-corpus
3https://dumps.wikimedia.org/enwiki/20221220/
4https://dumps.wikimedia.org/simplewiki/20221201/

##  3 Tracks

This shared task includes three tracks: Strict, Strict-small, and Loose.

The Strict and Strict-small tracks require that submissions are trained
exclusively on a fixed dataset, which we provide. The main difference between
these tracks is the size of the dataset ($\sim$10M words vs. $\sim$100M
words). Both datasets contain child-directed speech, transcribed speech from
multiple sources, children’s books, and Wikipedia, among other datasets. The
Strict-small dataset is an approximately 10% uniform subsample of the Strict
dataset. See §4 for a full description of the fixed datasets. Winners will be
determined based on performance on the shared evaluation set.

The Loose track relaxes these restrictions. Submissions must still be trained
on a maximum of 100M words, and will be tested on the shared evaluation set.
However, they are permitted to use unlimited non-linguistic data or text which
differs from the restricted shared task. Training on additional text is
allowed without limits if that text is generated by a model trained following
the above restrictions. For this track, winners will be selected holistically
based on evaluation performance, relevance to the shared task goals, potential
impact, and novelty.

##  4 Dataset

We distribute a developmentally plausible pretraining dataset inspired by the
input to children.111Clicking on the following link will download the dataset
(240MB zipped, 700MB unzipped):
https://github.com/babylm/babylm.github.io/raw/main/babylm_data.zip
Submissions must use only this training data to be considered for the
Strict(-small) tracks, but may use different data for the Loose track. The
dataset has two key properties:

  * •

Under 100M words: Children are exposed to 2M-7M words per year (Gilkerson et
al., 2017). Choosing the beginning of adolescence (age 12) as a cutoff, the
dataset should be between 24M-84M words.

  * •

Mostly transcribed speech: Most of the input to children is spoken. Thus, we
include a higher proportion of transcribed speech in our dataset.

The datasets we release are mixed domain, taken from multiple sources. Table 1
summarizes the composition of the datasets.

##  5 Evaluation

We will distribute a shared evaluation pipeline based in Google Colab. Colab
provides access to relatively small GPUs; this will allow users from various
research settings of varying resources to efficiently evaluate their
submissions. Our evaluation code will also be public, such that those wishing
to use their own computational resources may do so. More details about the
evaluation pipeline and the set of tasks will be released subsequently.

The pipeline assumes all models can be loaded and queried in HuggingFace’s
transformers library (Wolf et al., 2020).222While discouraged, participants
whose models are not compatible with the transformers library can still
conduct the necessary evaluation through their own pipeline. Additionally, all
models must be able to score a sequence—e.g., assign a log-likelihood or
pseudo log-likelihood Wang and Cho (2019); Salazar et al. (2020)—and must be
able to be fine-tuned to perform classification tasks. Models do not need to
be able to generate sequences. Submissions must include model outputs for each
of the core evaluations in a format that we specify in our evaluation
pipeline.

We choose evaluations that represent the core interests of this shared task,
focusing on efficiency and applied NLP, as well as cognitive science,
linguistics and language acquisition. Especially good performance in one but
not both of these areas may be acknowledged with a special award.

###  5.1 Baselines

We will also release a series of baseline models with the evaluation pipeline.
To train these, we simply take the hyperparameters from a series of
established large language models and train them from scratch on our fixed
datasets. We use hyperparameters from OPT (decoder-only; Zhang et al., 2022),
RoBERTa (encoder-only; Liu et al., 2019), and T5 (encoder-decoder; Raffel et
al., 2020). These are not meant to be strong baselines, but rather to provide
a naïve starting point for improving language models for this domain.

##  6 Submissions

HTML]E9F0E9  What you Need to Submit • A link where we can download the model
• A .zip of predictions (from our eval pipeline) • A short description of the
approaches taken • If Loose track: a link where we can download any additional
data

Although scaled-down pretraining is more accessible to research groups with
limited resources, pretraining is still expensive from a computational,
energy, and financial perspective. To help groups plan for total costs, we
will release an estimate of the resources required to pretrain on 10M words
and 100M words. For the Loose track, evaluation of submissions may take into
consideration computational efficiency as part of the holistic evaluation.

##  7 FAQs

#### Can papers be submitted to multiple tracks?

Yes. For example, a single paper can describe models which are submitted
separately to the Strict and Strict-small tracks.

#### Can I submit a paper about my work?

Yes, we encourage all participants to submit their reports, which will be
published in the proceedings of CoNLL. You may also describe any additional
experiments beyond those required for the shared task evaluation.

#### Can I submit additional evaluation metrics?

Yes, if you wish to submit your own evaluation metrics, along with model
performance, alongside our standardized evaluation results these can be
considered as part of the holistic evaluation in the Loose track.

#### What training regimes are permitted?

For the Strict/Strict-small tracks, any kind of training objective/regime is
permitted, as long as the data restrictions are followed. Pretrained models
may not be used for any purpose such as reranking or data augmentation.

We do however require for evaluation purposes that the model provides a
function to score a sequence—e.g., log-likelihood for autoregressive models or
pseudo-log-likelihood for masked language models—without the need for
additional fine-tuning.

#### Are there any limits on hyperparameters?

No. In the Loose track, parameter efficiency and training efficiency may be
considered along with other factors in ranking submissions.

#### Are there any limits on the number of epochs?

No. We put no restrictions on the number of epochs, for several reasons:
First, from an engineering perspective, training LMs with SGD tends to require
multiple epochs at these scales to achieve peak performance. Second, from a
cognitive perspective, humans have a memory of linguistic experience, and can
continue to access and learn from these memories. Third, we try not to make a
stand on implementations to allow the most freedom for innovation.

##  8 Organizing Committee

Leshem Choshen | Aaron Mueller
---|---
Ryan Cotterell | Alex Warstadt
Kundan Krishna | Ethan Wilcox
Tal Linzen | Adina Williams
Haokun Liu | Chengxu Zhuang

Questions? Feel free to contact us at the following email addresses:
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
## References

  * Abdelali et al. (2014) Ahmed Abdelali, Francisco Guzman, Hassan Sajjad, and Stephan Vogel. 2014.  The AMARA corpus: Building parallel language resources for the educational domain.  In _Proceedings of the Ninth International Conference on Language Resources and Evaluation (LREC’14)_ , Reykjavik, Iceland. European Language Resources Association (ELRA).
  * Brown et al. (2020) Tom Brown, Benjamin Mann, Nick Ryder, Melanie Subbiah, Jared D Kaplan, Prafulla Dhariwal, Arvind Neelakantan, Pranav Shyam, Girish Sastry, Amanda Askell, Sandhini Agarwal, Ariel Herbert-Voss, Gretchen Krueger, Tom Henighan, Rewon Child, Aditya Ramesh, Daniel Ziegler, Jeffrey Wu, Clemens Winter, Chris Hesse, Mark Chen, Eric Sigler, Mateusz Litwin, Scott Gray, Benjamin Chess, Jack Clark, Christopher Berner, Sam McCandlish, Alec Radford, Ilya Sutskever, and Dario Amodei. 2020.  Language models are few-shot learners.  In _Advances in Neural Information Processing Systems_ , volume 33, pages 1877–1901. Curran Associates, Inc.
  * Chowdhery et al. (2022) Aakanksha Chowdhery, Sharan Narang, Jacob Devlin, Maarten Bosma, Gaurav Mishra, Adam Roberts, Paul Barham, Hyung Won Chung, Charles Sutton, Sebastian Gehrmann, Parker Schuh, Kensen Shi, Sasha Tsvyashchenko, Joshua Maynez, Abhishek Rao, Parker Barnes, Yi Tay, Noam Shazeer, Vinodkumar Prabhakaran, Emily Reif, Nan Du, Ben Hutchinson, Reiner Pope, James Bradbury, Jacob Austin, Michael Isard, Guy Gur-Ari, Pengcheng Yin, Toju Duke, Anselm Levskaya, Sanjay Ghemawat, Sunipa Dev, Henryk Michalewski, Xavier Garcia, Vedant Misra, Kevin Robinson, Liam Fedus, Denny Zhou, Daphne Ippolito, David Luan, Hyeontaek Lim, Barret Zoph, Alexander Spiridonov, Ryan Sepassi, David Dohan, Shivani Agrawal, Mark Omernick, Andrew M. Dai, Thanumalayan Sankaranarayana Pillai, Marie Pellat, Aitor Lewkowycz, Erica Moreira, Rewon Child, Oleksandr Polozov, Katherine Lee, Zongwei Zhou, Xuezhi Wang, Brennan Saeta, Mark Diaz, Orhan Firat, Michele Catasta, Jason Wei, Kathy Meier-Hellstern, Douglas Eck, Jeff Dean, Slav Petrov, and Noah Fiedel. 2022\.  PaLM: Scaling language modeling with pathways.  _Computing Research Repository_ , arXiv:2204.02311.
  * Dupoux (2018) Emmanuel Dupoux. 2018.  Cognitive science in the era of artificial intelligence: A roadmap for reverse-engineering the infant language-learner.  _Cognition_ , 173:43–59.
  * Gerlach and Font-Clos (2018) Martin Gerlach and Francesc Font-Clos. 2018.  A standardized Project Gutenberg corpus for statistical analysis of natural language and quantitative linguistics.  _Computing Research Repository_ , arXiv:1812.08092.
  * Gilkerson et al. (2017) Jill Gilkerson, Jeffrey A Richards, Steven F Warren, Judith K Montgomery, Charles R Greenwood, D Kimbrough Oller, John HL Hansen, and Terrance D Paul. 2017\.  Mapping the early language environment using all-day recordings and automated analysis.  _American journal of speech-language pathology_ , 26(2):248–265.
  * Hill et al. (2016) Felix Hill, Antoine Bordes, Sumit Chopra, and Jason Weston. 2016.  The Goldilocks principle: Reading children’s books with explicit memory representations.
  * Hoffmann et al. (2022) Jordan Hoffmann, Sebastian Borgeaud, Arthur Mensch, Elena Buchatskaya, Trevor Cai, Eliza Rutherford, Diego de Las Casas, Lisa Anne Hendricks, Johannes Welbl, Aidan Clark, et al. 2022.  Training compute-optimal large language models.  _arXiv preprint arXiv:2203.15556_.
  * Keller (2010) Frank Keller. 2010.  Cognitively plausible models of human language processing.  In _Proceedings of the ACL 2010 Conference Short Papers_ , pages 60–67, Uppsala, Sweden. Association for Computational Linguistics.
  * Lison and Tiedemann (2016) Pierre Lison and Jörg Tiedemann. 2016.  OpenSubtitles2016: Extracting large parallel corpora from movie and TV subtitles.  In _Proceedings of the Tenth International Conference on Language Resources and Evaluation (LREC’16)_ , pages 923–929, Portorož, Slovenia. European Language Resources Association (ELRA).
  * Liu et al. (2019) Yinhan Liu, Myle Ott, Naman Goyal, Jingfei Du, Mandar Joshi, Danqi Chen, Omer Levy, Mike Lewis, Luke Zettlemoyer, and Veselin Stoyanov. 2019.  RoBERTa: A robustly optimized bert pretraining approach.  _Computing Research Repository_ , arXiv:1907.11692.
  * MacWhinney (2000) Brian MacWhinney. 2000.  _The CHILDES project: The database_ , volume 2.  Psychology Press.
  * Raffel et al. (2020) Colin Raffel, Noam Shazeer, Adam Roberts, Katherine Lee, Sharan Narang, Michael Matena, Yanqi Zhou, Wei Li, and Peter J. Liu. 2020.  Exploring the limits of transfer learning with a unified text-to-text transformer.  _Journal of Machine Learning Research_ , 21(140):1–67.
  * Salazar et al. (2020) Julian Salazar, Davis Liang, Toan Q. Nguyen, and Katrin Kirchhoff. 2020.  Masked language model scoring.  In _Proceedings of the 58th Annual Meeting of the Association for Computational Linguistics_ , pages 2699–2712, Online. Association for Computational Linguistics.
  * Stolcke et al. (2000) Andreas Stolcke, Klaus Ries, Noah Coccaro, Elizabeth Shriberg, Rebecca Bates, Daniel Jurafsky, Paul Taylor, Rachel Martin, Marie Meteer, and Carol Van Ess-Dykema. 2000.  Dialogue act modeling for automatic tagging and recognition of conversational speech.  _Computational Linguistics_ , 26(3):339–371.
  * Wang and Cho (2019) Alex Wang and Kyunghyun Cho. 2019.  BERT has a mouth, and it must speak: BERT as a Markov random field language model.  In _Proceedings of the Workshop on Methods for Optimizing and Evaluating Neural Language Generation_ , pages 30–36, Minneapolis, Minnesota. Association for Computational Linguistics.
  * Warstadt and Bowman (2022) Alex Warstadt and Samuel R Bowman. 2022.  What artificial neural networks can tell us about human language acquisition.  In _Algebraic Structures in Natural Language_ , pages 17–60. CRC Press.
  * Wolf et al. (2020) Thomas Wolf, Lysandre Debut, Victor Sanh, Julien Chaumond, Clement Delangue, Anthony Moi, Pierric Cistac, Tim Rault, Rémi Louf, Morgan Funtowicz, Joe Davison, Sam Shleifer, Patrick von Platen, Clara Ma, Yacine Jernite, Julien Plu, Canwen Xu, Teven Le Scao, Sylvain Gugger, Mariama Drame, Quentin Lhoest, and Alexander M. Rush. 2020.  Transformers: State-of-the-art natural language processing.  In _Proceedings of the 2020 Conference on Empirical Methods in Natural Language Processing: System Demonstrations_ , pages 38–45, Online. Association for Computational Linguistics.
  * Zhang et al. (2022) Susan Zhang, Stephen Roller, Naman Goyal, Mikel Artetxe, Moya Chen, Shuohui Chen, Christopher Dewan, Mona Diab, Xian Li, Xi Victoria Lin, Todor Mihaylov, Myle Ott, Sam Shleifer, Kurt Shuster, Daniel Simig, Punit Singh Koura, Anjali Sridhar, Tianlu Wang, and Luke Zettlemoyer. 2022.  OPT: Open pre-trained transformer language models.  _Computing Research Repository_ , arXiv:2205.01068.


# A Log-Linear Non-Parametric Online Changepoint Detection Algorithm based on
Functional Pruning

Gaetano Romano, Idris Eckley, Paul Fearnhead

###### Abstract

Online changepoint detection aims to detect anomalies and changes in real-time
in high-frequency data streams, sometimes with limited available computational
resources. This is an important task that is rooted in many real-world
applications, including and not limited to cybersecurity, medicine and
astrophysics. While fast and efficient online algorithms have been recently
introduced, these rely on parametric assumptions which are often violated in
practical applications. Motivated by data streams from the telecommunications
sector, we build a flexible nonparametric approach to detect a change in the
distribution of a sequence. Our procedure, NP-FOCuS, builds a sequential
likelihood ratio test for a change in a set of points of the empirical
cumulative density function of our data. This is achieved by keeping track of
the number of observations above or below those points. Thanks to functional
pruning ideas, NP-FOCuS has a computational cost that is log-linear in the
number of observations and is suitable for high-frequency data streams. In
terms of detection power, NP-FOCuS is seen to outperform current nonparametric
online changepoint techniques in a variety of settings. We demonstrate the
utility of the procedure on both simulated and real data.

##  1 Introduction

One of the contemporary challenges in time-series analysis is to detect
changes in some measurable properties of a process. This task finds its roots
in a plethora of applications spanning numerous fields including engineering
(Alvarez-Montoya et al.,, 2020; Henry et al.,, 2010), neuroscience (Jewell et
al.,, 2020), genomics Nicolas et al., (2009) and astrophysics (Fridman,, 2010;
Fuschino et al.,, 2019). In the previous decade, we saw many offline
changepoint procedures appearing in the statistical literature. A common
approach for detecting a single change is to scan over all possible locations
for a change, and apply a generalised likelihood ratio (GLR) test for a
change, with evidence for a change being the maximum of the GLR test
statistics Fearnhead and Fryzlewicz, (2022). The GLR procedure has been proven
to be asymptotically optimal (Basseville et al.,, 1993), as it searches over
all possible parameters of the unknown pre and post-change distributions. This
approach includes, as a special case, the popular CUSUM method for detecting a
change in mean, and can be extended to detect multiple changepoints by using
binary segmentation methods Scott and Knott, (1974); Fryzlewicz and Rao,
(2014); Kovács et al., (2020) or by maximising a penalised likelihood
Fearnhead and Rigaill, (2020); see Cho and Kirch, (2021) for a recent review
of this area.

In contrast to the offline setting, where the data is first collected and
later analyzed _a posteriori_ , one of the modern challenges is to detect a
change within a data stream in real time. We find many applications in need of
an online procedure. These include control of industrial processes (Pouliezos
and Stavrakakis,, 2013), or monitoring of computer networks (Tartakovsky et
al.,, 2005), social networks (Chen,, 2019) and telecommunication devices
(Austin et al.,, 2023).

Many challenges arise in analysing data online. An online procedure should be
able to process observations in real-time, as quickly as they arrive, in order
to avoid memory overflow and result in a delayed evaluation. This can be
particularly valuable in settings where limited computational resources are
available, or in the high-frequency domain.

As noted by (Ross,, 2015), one way of performing an online analysis is to
collect data in batches and analyse those through an offline algorithm.
However, such an approach can be sensitive to the batch size. If this is too
small, then we may be unable to detect small changes, while if it is large it
can lead to delayed detection of bigger changes. Alternatively, observations
can be processed on the go in a sequential fashion. At each new observation, a
decision is made on whether a change has occurred, or not, based on the new
data point and on past information.

The sequential CUSUM approach, and more generally, sequential GLR approaches,
demonstrate excellent statistical properties Wang and Xie, (2022); Yu et al.,
(2020). However, they can be computationally inefficient. A naive computation
of the GLR test, in fact, involves, at time $n$, considering $\mathcal{O}(n)$
possible locations for a change. I.e., the algorithm has a computational cost
per iteration that increases linearly, which is impracticable for online
applications.

Recently, Romano et al., (2023) presented the FOCuS procedure, a fast
algorithm to perform the sequential CUSUM test, that decreases the
computational complexity from $\mathcal{O}(n)$ to $\mathcal{O}(\log n)$ per
iteration. This is to our knowledge the fastest way to solve the sequential
CUSUM test exactly. For example, the expected cost per iteration of FOCuS at a
given time, say at iteration one million, is roughly equal to the cost of
iterating approximately 20 objects stored in memory to find the global maximum
of our statistics – and thus it is suitable for high-frequency online
applications.

As previously mentioned, the CUSUM and GLR tests in general rely on parametric
assumptions that are often hard to meet in real-world applications. For
instance, with an industrial collaborator, we found that certain data streams
arising in the telecommunications sector do not follow common distributions,
nor do they satisfy the usual parametric assumptions. In monitoring
operational metrics from network devices, for example, we often deal with
nonstandard distributions with multiple modes, outliers or heavy tails. In
many cases, the nature of the change is unknown _a priori_ , see for example
the sequences from Figure 8. In Section 4.2 we present one example of a
contemporary telecommunication application, where the aim is to monitor the
operational performance of devices on an optical cable network. In those
scenarios, a test for Gaussian change-in-mean would, in fact, be prone to
either a missed detection, in the case of a change of a different nature, or
to false positives, in case of misspecification of the underlying noise
process.

To better outline the limitations of a parametric approach in online
changepoint detection, we show a simple introductory example. Let us compare
the Gaussian FOCuS from Romano et al., (2023) with its non-parametric
counterpart, NP-FOCuS, the methodology that will be introduced in this paper.
In Figure 1 we study both the statistics on 3 different simulated changepoint
scenarios. Thresholds were tuned to achieve comparable average run lengths of
2000 observations under the null (the expected number of observations until we
detect a change if a change is not present). This way, after placing a true
change at 1000, we can have a fair comparison focusing on the detection alone.
In Figure 1a, the simple Gaussian change-in-mean case, we notice how the
Gaussian procedure achieves the fastest detection. This is because the
procedure’s assumptions hold. However in Figure 1b, the second example, the
sequence now shows a change in variance. In this setting, the Gaussian change-
in-mean has no power to detect the change, which is missed. Lastly, in Figure
1c, we find the same series as in the first example, but some observations are
shifted up by 3. Under this scenario, to achieve the same run length of 2000
observations under the null, the threshold for the Gaussian FOCuS procedure
needs to be inflated. This results in a slower detection than its non-
parametric counterpart.

Figure 1: Detection delays of Gaussian FOCuS (in red) and NP-FOCuS (in green)
for 3 sequences of length 2000 with a change at $t=1000$. In a, a Gaussian
change-in-mean example, with a jump size of 0.25. In b, a change-in-variance
from variance $\sigma^{2}=1$ to $\sigma^{2}=0.25$. Lastly, in c, we find the
same scenario in a with some of the observations shifted by a value of 3 at
uniformly random sampled locations (blue dots).

Whilst there are numerous offline non-parametric approaches in the current
literature (including and not limited to Pettitt,, 1979; Zou et al.,, 2014;
Haynes et al.,, 2017; Matteson and James,, 2014; Chen and Chu,, 2023), it is
only more recently that there has been substantial interest on online
approaches. This includes Ross and Adams, (2012), who propose a control-chart
approach that allows for the detection of changes in scale or location, or
both, when the underlying distribution of the process is unknown, and Shin et
al., (2022) who propose a general framework for detecting changes in
measurable properties of a process. More generally there are methods that
choose a number of univariate summaries of data and monitor for a change in
any of these summaries Kurt et al., (2020); Keriven et al., (2020), kernel-
based methods Flynn and Yoo, (2019); Huang et al., (2014) or related methods
that use information on distances between data points Chen, (2019); Chu and
Chen, (2022); Chen et al., (2020). These approaches can be applied in
multivariate settings, but their statistical efficiency depends on choosing
appropriate choice of summary, kernel or distance that had power to detect the
type of change that occurs. These methods are often not invariant to
transformations of the data. By comparison, methods based on ranking
information (e.g. Hawkins and Deng,, 2010), are simple to apply for univariate
data and are invariant to any monotone transformation of the data. This is
also the case for methods that look for a change in the distribution function
for univariate data, which is the approach we take in this paper.

One example of such a method is the NUNC algorithm Austin et al., (2023),
which tests for a change in the empirical cumulative distribution function
(eCDF) within a rolling window. As with any other moving window method,
however, NUNC’s performance is extremely dependent on the size of the window.
For example, a window too small would end up missing changes of a smaller
size, whilst a window too big might result in longer detection delays and an
increased computational cost.

Our procedure, NP-FOCuS, tries to mitigate such issues by building a GLR test
(extending ideas from Romano et al.,, 2023) for the non-parametric procedure
NUNC of Austin et al., (2023). Our procedure maps the change-in-distribution
problem into a Bernoulli GLR test. This is achieved by evaluating the eCDF at
a fixed point and keeping track of the number of observations above or below
that point. We then extend the FOCuS algorithm so that it can apply to the
Bernoulli GLR test, and show theoretically that it has the same strong
computational properties as the original FOCuS algorithm. We perform this test
for a grid of quantiles, and then merge the test statistic values. In practice
our final procedure monitors both the sum of the test statistics across the
quantiles, which can detect small shifts in the distribution, and also the
maximum of the test statistics, which can detect a larger change in just one
part of the distribution. Our approach only assumes that the observations are
i.i.d.. However, we also demonstrate empirically that even if such a
hypothesis is violated we still retain the strong performance over competitor
methods.

The paper is structured as follows. In Section 2 we introduce the NP-FOCuS
procedure, starting from related work, to the introduction of a novel
functional pruning recursion. In Section 3 we provide some guarantees on the
computational complexity of NP-FOCuS. In Section 4.1 we evaluate empirically
performances of NP-FOCuS, with a comparison to other methods in a simulation
study. In Section 4.2 we evaluate the method on some real-word data streams.

Open-source code, making the methods developed in this paper accessible, will
be made available via an appropriate repository in due course.

##  2 A Functional Representation for the Nonparametric Cost Function

###  2.1 Overview of Problem and Approach

We aim to build an online non-parametric changepoint detector by monitoring
the empirical Cumulative Density Function (eCDF) of a stream of observations.
This approach has been shown to work well in an offline setting by Zou et al.,
(2014); Haynes et al., (2017).

Let $p\in\mathbb{R}$ and let $F_{1:n}(p)$ be the unknown CDF for a series of
independent real-valued observations $y_{1},\dots,y_{n}$, and let
$\hat{F}_{1:n}(p)$ be its eCDF:

$\hat{F}_{1:n}(p)=\frac{1}{n}\sum_{t=1}^{n}\left[\mathbb{I}(y_{t}\leq
p)\right],$

where $\mathbb{I}$ is an indicator function. For an independent stream of
observations, for a fixed $p$, the eCDF follows a binomial distribution:

$n\hat{F}_{1:n}(p)\sim\mathrm{Binom}(n,\theta),$

with $\theta=F_{1:n}(p)$. A simple approach to detect a change in
$\hat{F}_{1:n}(p)$ would be to record which data points are above $p$ and
build a likelihood-ratio test for the Bernoulli data. That is, for a fixed
$p$, we know that the likelihood for a segment is given by

$\begin{split}&\mathcal{L}(y_{\tau_{1}+1:\tau_{2}};\
\theta,p)=(\tau_{2}-\tau_{1})\times\\\
&\times\left[\hat{F}_{\tau_{1}+1:\tau_{2}}(p)\log\theta+(1-\hat{F}_{\tau_{1}+1:\tau_{2}}(p))\log(1-\theta)\right].\end{split}$
(1)

The GLR test for a change in parameter from $\theta_{0}$ to $\theta_{1}$ is

$\max_{0\leq\tau<n}\left[-\mathcal{L}(y_{\tau+1:n};\theta_{0},p)+\mathcal{L}(y_{\tau+1:n};\theta_{1},p)\right].$
(2)

By writing $x_{t}=\mathbb{I}(y_{t}\leq p)$, this can be rewritten as

$\displaystyle\mathcal{L}(y_{\tau+1:n};\theta_{1},p)-\mathcal{L}(y_{\tau+1:n};\theta_{0},p)$
$\displaystyle=$
$\displaystyle\sum_{t=\tau+1}^{n}\left[x_{t}\log\theta_{1}+(1-x_{t})\log(1-\theta_{1})\right]$
$\displaystyle-$
$\displaystyle\sum_{t=\tau+1}^{n}[x_{t}\log\theta_{0}+(1-x_{t})\log(1-\theta_{0})]$

If both $\theta_{0}$ and $\theta_{1}$ are known, then the sequential procedure
of Page, (1954) can be used to calculate the GLR test (2). In Section 2.2 we
show how to extend the FOCuS algorithm of Romano et al., (2023) to calculate
the GLR test if either only $\theta_{0}$ is known, or if neither $\theta_{0}$
nor $\theta_{1}$ is known. We will call this procedure Ber-FOCuS.

In practice, we would want to check changes in the entire eCDF rather than
just at a single quantile value $p$. To achieve this we follow Zou et al.,
(2014) and Haynes et al., (2017) who suggest summing up the test statistic
across a grid of values $p_{1},\dots,p_{M}$. Specifically, we approximate our
eCDF on a fixed grid of $M$ quantile values $p_{1},\dots,p_{M}$ computed over
a probation period. To obtain NP-FOCuS we run the Ber-FOCuS routine
independently for each quantile of our eCDF. To construct a global statistic
we aggregate the quantile statistics through the sum and the maximum. In
Section 2.5 we will formally describe the full procedure.

###  2.2 Detecting change-in-rate in a Bernoulli process

We focus on detecting a change in the rate parameter $\theta$ in a univariate
stream of a Bernoulli process. Let $x_{n}\in\\{0,1\\}\sim\text{Ber}(\theta)$
be a realization of a Bernoulli random variable with parameter $\theta$, with
$\theta$ subject to a change. Assume the pre-change rate parameter,
$\theta_{0}$, is known. For every observation we can obtain evidence for a
post-change rate parameter $\theta_{1}$ via the likelihood ratio test:

$\displaystyle g(x_{n},\theta_{1})$
$\displaystyle=-2\log\left(\frac{\theta_{0}^{x_{n}}(1-\theta_{0})^{1-x_{n}}}{\theta_{1}^{x_{n}}(1-\theta_{1})^{1-x_{n}}}\right)=$
(3)
$\displaystyle=2\left[x_{n}\log\left(\frac{\theta_{1}}{\theta_{0}}\right)+(1-x_{n})\log\left(\frac{1-\theta_{1}}{1-\theta_{0}}\right)\right].$
(4)

At time $n$ we have observed $x_{1},\dots,x_{n}$, we can then employ (4) and
derive the test statistics:

$\mathcal{Q}_{n,\theta_{1}}=\max_{0\leq\tau<n}\
\sum_{t=\tau+1}^{n}g(x_{t},\theta_{1}),$

where we maximise for all possible time-points of the change, $\tau$. Such a
statistic is known as the sequential Page-CUSUM statistic (Page,, 1954), and
can be solved efficiently in constant time per iteration through the
recursion:

$\mathcal{Q}_{n,\theta_{1}}=\max\left\\{0,\
\mathcal{Q}_{n-1,\theta_{1}}+g(x_{n},\theta_{1})\right\\}.$ (5)

This statistic is straightforward to compute for a known value of $\theta_{0}$
and a fixed value of $\theta_{1}$. As we do not know $\theta_{1}$, one would
compute the statistics over a fixed grid of values for $\theta_{1}$ and choose
the maximum value of those to maximise the power of a detection. However, this
would introduce a trade-off between power and computational complexity.
Alternatively, Xie et al., (2023) recently suggested a way of approximating
the Page-CUSUM cost by using a plug-in value for $\theta_{1}=\hat{\theta}$
into the equation above, that is calculated recursively based on only recent
data point.

Our approach is to calculate the GLR test statistic exactly for both the case
where $\theta_{1}$ is unknown. The idea, based on Romano et al., (2023), is to
solve (5) simultaneously for all values of $\theta$ through the functional
representation of the sequential test statistics $Q_{n}(\theta)$ for a post-
change rate parameter $\theta$. That is, we have $Q_{0}(\theta)=0$ and:

$Q_{n}(\theta)=\max\left\\{0,\ Q_{n-1}(\theta)+g(x_{n},\theta)\right\\}$ (6)

for $n\geq 1$. The value of our test statistic will be
$\mathcal{Q}_{n}=\max_{\theta}Q_{n}(\theta)$. In this optimization lies the
major computational contribution of NP-FOCuS. At each iteration of the
algorithm we explicitly compute and store the full functional representation
of our cost $Q_{n}(\theta)$, which will be a piecewise function (see Figure
2). By maximizing over all possible values of $\theta\in\ (0,1)$, which is
simple as we can maximise each piecewise part of the function and take the
maximum of these maxima, we can find the highest possible value for the Page-
CUSUM statistics, and avoid altogether any choice of window $w$ or specify a
post-change parameter.

Figure 2: The functional cost $Q_{n}(\theta)$ at iteration $n$. Dotted lines
illustrate (i) $\text{arg}\max_{\theta}Q_{n}(\theta)$ and (ii) the
$\max_{\theta}Q_{n}(\theta)$, the highest possible value of our test
statistics.

At a given iteration $Q_{n}(\theta)$ is obtained from a set of component
functions of the form:

$q^{(n,\tau)}(\theta)=a^{(n,\tau)}\log\left(\frac{\theta}{\theta_{0}}\right)+b^{(n,\tau)}\log\left(\frac{1-\theta}{1-\theta_{0}}\right),$
(7)

with $a^{(n,\tau)}=\sum_{t=\tau}^{n}x_{t}$,
$b^{(n,\tau)}=\sum_{t=\tau}^{n}(1-x_{t})$. This curve represents the LR
statistic for a change at $\tau$.

Whilst there are $n$ such possible curves, each one relative to a different
candidate changepoint $\tau$, in practice $Q_{n}(\theta)$ can be defined in
term of a much smaller number of curves. In other words, for a large $n$, only
a small proportion of these curves will contribute to the functional cost,
i.e. will be the components of our piecewise function $Q_{n}(\theta)$, such
that for some values of $\theta$ we have that
$q^{(n,\tau)}(\theta)=Q_{n}(\theta)$. In our example, we plot in Figure 3 the
the curves contributing to Figure 2. If there are no values of $\theta$ for
which a curve contributes to the optimal cost, we can drop it, a mechanism
that is known as pruning that leads to the speed-ups in the algorithm.
Therefore, Ber-FOCuS presents a way of iteratively and efficiently updating
$Q_{n}(\theta)$, figuring out which are its components’ functions and finding
its global maximum to solve 6 exactly.

Figure 3: The components function of the functional cost $Q_{n}(\theta)$ at
iteration $n$. The functions are plotted in different colours as dotted lines.
The values for which they are optimal, e.g. they dominate all other lines and
contribute to $Q_{n}(\theta)$, are solid. These curves are in line with those
of the bottom right plot of in Fig. 4, illustrating the last step of our
algorithm for $\theta>\theta_{0}=0.5$.

If $k$ is the number of functions that contribute to $Q_{n}(\theta)$, we can
represent $Q_{n}(\theta)$ on a computer efficiently as an ordered set of $k$
triples $\mathrm{Q}_{n}=\\{q_{i}=(a_{i}^{(n,\tau_{i})},\
b_{i}^{(n,\tau_{i})},l_{i}^{(n,\tau_{i})}),\ \ i=1,\dots,\ k\\}$. Each triple
represents a curve $q^{(n,\tau_{i})}(\theta)$, and contains the coefficients
of $\log(\theta/\theta_{0})$ and of $\log((1-\theta)/(1-\theta_{0}))$, and a
value $l_{i}$ which is the left-most value of $\theta$ for which the
associated function is contributing to the optimal cost, i.e.
$\theta\in[l_{i}^{(\tau_{i},n)},l_{i+1}^{(\tau_{i+1},n)})$,
$q_{i}^{(n,\tau_{i})}(\theta)=Q_{n}(\theta)$. In this interval our curve will
be greater than all the other curves and the zero line.

For brevity, for the rest of this section, we will fix the iteration. We can
simplify the notation by dropping the superscript $(n,\tau)$, and hence we
will denote our triples and functions simply as $q_{1},\ldots,q_{k}$. This is
possible as $n$ is fixed and $\tau$ is redundant information: by construction
$\tau_{i}=n-(a_{i}+b_{i})$ as we are dealing with binary data. We will assume
that the triples are ordered in such a way that
$\theta_{0}=l_{1}<\dots<l_{k}<l_{k+1}=1$, so that for all
$\theta\in[l_{i},l_{i+1})$,
$Q_{n}(\theta)=a_{i}\log(\theta/\theta_{0})+b_{i}\log((1-\theta)/(1-\theta_{0})$.

A description of the full Bernoulli FOCuS (Ber-FOCuS) procedure is given in
Algorithm 1, with a graphical illustration of the procedure in Figure 4. For
each iteration, we have three main steps. Step 1, the update step, updates the
values of the coefficients of stored functions according to (4). I.e., we
update the count of positive or negative observations since the introduction
of a curve. This is a simple increment of either coefficient $a_{i}$ or
$b_{i}$. Step 2, which will be detailed in Section 2.3, is the functional
pruning step and focuses on identifying which pieces are no longer optimal
(and will never be in the future). Finally, the optimization step, where we
maximise the total cost $Q_{n}(\theta)$ over $\theta$, to get our test
statistics. This is equivalent to taking the maximum of all maximums of our
functions. As we will see in Section 3 to reconstruct the optimal cost we need
only to store a small subset of candidate changepoints: at the $n^{th}$
iteration we expect to store on average $k\approx\log(n)$ functions.

Data: $x_{n}\in\\{0,1\\}$ the data point at time $n$.

Input: $Q_{n-1}(\theta)$ the cost function from the previous iteration.

1

2$\tilde{Q}(\theta)\xleftarrow{}Q_{n-1}(\theta)+g(x_{n},\theta)$

$Q_{n}(\theta)\xleftarrow{}\max\left\\{0,\ \tilde{Q}(\theta)\right\\}$ ;

// see Algorithm 2

3

4$\mathcal{Q}_{n}\xleftarrow{}\max_{\theta}Q_{n}(\theta)$

return _$Q_{n}(\theta)$ for the next iteration, $\mathcal{Q}_{n}$ as the test
statistic._

Algorithm 1 Ber-FOCuS (one iteration) Figure 4: One iteration of the NP-FOCuS
procedure illustrated, for $\theta>\theta_{0}$. In this example, we assume
$\theta_{0}=0.5$. Top-left: the cost function from the previous iteration
$Q_{n-1}(\theta)$. The total cost is given by the maximum of 4 functions: we
can see how each of those is optimal for some values of $\theta$. Labels refer
to the time $\tau$ at which each function was introduced. Top-right: the
update step. By adding the new observation $x_{n}$, we update the coefficients
$a_{i,\tau}$ or $b_{i,\tau}$ of each function. Although the shape of our
function has changed, the points where quadratics intersect are unchanged.
Bottom-left, the pruning step. Here we compare the functions with the zero
line (dashed line). We find the lowest value of $\theta$ for which each
function is optimal. We notice how for the grey function (1000) there is no
such value of $\theta$ for which the function is optimal, as the intersection
with its closest quadratic falls below the zero line. We can therefore safely
prune the function. Bottom-right: $Q_{n}(\theta)$, after we add a new piece
(in red).

###  2.3 Step 2: Efficient Pruning of the Bernoulli cost

Functional pruning consists of restricting the set of values to consider as
candidate changepoints on account of the information they carry up to time
$n$. That is, at each iteration we can stop considering functions introduced
at times if they will never be optimal for our cost $Q_{n}(\theta)$. To
achieve this, Romano et al., (2023) propose an efficient pruning rule, we can
adapt to our setting. As in Romano et al., (2023), we propose to update
$Q_{n}(\theta)$ separately for $\theta>\theta_{0}$ and $\theta<\theta_{0}$. We
will report the update rule for $\theta>\theta_{0}$ alone as, by symmetry, the
update for $\theta<\theta_{0}$ can be found by applying the same algorithm
with the data inverted, i.e. replacing $x_{t}$ with $(1-x_{t})$.

If we consider Ber-FOCuS algorithm of Algorithm 1, then at the start of an
iteration we will have a set of triples. The first step, the update step, only
affects the $a_{i}$ and $b_{i}$ coefficients. The $l_{i}$ components, which
specify the intervals of $\theta$ that each curve is optimal, are unchanged.
This follows as the differences between the curves are unaffected by the
update, as each curve is changed by adding the same function
$g(x_{n},\theta)$, and it is the difference between curves that determines
where one if greater than all others.

The pruning of curves only occurs in step 2, and is as a result of taking the
maximum of $\tilde{Q}(\theta)$ and the zero line. The idea is that we need to
find the values of $\theta$ for which the zero line is optimal. Define
$\textsf{root}(q_{i})$ to be the largest value of $\theta>\theta_{0}$ such
that $q_{i}(\theta)=0$, if such a value exists. It is straightforward to show
that, for a given function $q_{i}$ that contributes to $\tilde{Q}(\theta)$ for
$\theta>\theta_{0}$, if $\textsf{root}(q_{i})<1$ exists, then the zero-line is
better on the interval:

$\left[\textsf{root}(q_{i}),\ 1\right).$

Considering all functions, this implies that if $\max_{\tau}\
\textsf{root}(q_{i})<1$ then the zero line is globally optimal on $[l^{*},1)$,
where

$l^{*}=\max_{\tau}\ \textsf{root}(q_{i}).$

Therefore, the triple defining the zero line will be given by $(0,0,l^{*})$.
Given that the quadratics are ordered, it means that any quadratic with
$l_{i}<l^{*}$ can therefore be removed. Finally, the ordering of the
quadratics means that $l_{i}<l_{j}$ if and only if $i<j$. For this reason, we
can start checking from the last quadratic $q_{k}$ and stop as soon as the
pruning condition is not met. The pruning procedure is summarised in Algorithm
2. Following the update, for pruning, given that at any point in time
$l_{i}=\textsf{root}(q_{i}-q_{i-1})\leq l^{*}$ and that $q_{i}(\theta_{0})=0$,
we know that whether $q_{i}(l_{i})<0$ there will be no such value of $\theta$
for which $q_{i}$ will be optimal. In fact, we know that $q_{i}$ will be
greater than $q_{i-1}$ only starting from $l_{i}$, and if $q_{i}(l_{i})<0$
falls below the zero line, then $q_{i-1}(\theta)<q_{i}(\theta)<0\ \forall\
\theta\in(l_{i},1)$ and therefore $q_{i}$ will be never be optimal and can be
safely pruned. The bottom-left plot of Figure 4 can be of aid in understanding
the pruning and the following theorems.

Input: $\tilde{Q}$, an ordered set of triples $\\{q_{i}=(a_{i},\ b_{i},l_{i})\
\forall\ i=1,\dots,\ k\\}$

1

2$i\xleftarrow{}k$;

3 while  _$q_{i}(l_{i})\leq 0$ and $i\geq 1$_ do

4

5 $i\xleftarrow{}i-1$;

6

7 end while

8

9if  _$i\neq k$_ then

$\mathrm{Q}\xleftarrow{}\mathrm{Q}\setminus\\{q_{i+1},\ \dots,q_{k}\\}$;

// drop pieces

10 $k\xleftarrow{}i$

11 end if

12

13if  _k = 0_ then

$q_{k+1}\xleftarrow{}(0,0,\theta_{0})$;

// add first piece

14

15 else

$q_{k+1}\xleftarrow{}(0,0,\textsf{root}(q_{k}))$;

// add new piece

16

17 end if

18$\mathrm{Q}\xleftarrow{}\\{\mathrm{Q},q_{k+1}\\}$

return _$\mathrm{Q}$_

Algorithm 2 Algorithm for $\max\\{0,\ \tilde{Q}(\theta)\\}$ for
$\theta>\theta_{0}$

We therefore learn that the algebraically ideal version of the algorithm
relies on computing the value $\textsf{root}(q_{k})$ in step 9 at each
iteration. However, given the shape of our cost function, we do not have a
closed form for computing this root. A naive implementation would thus require
obtaining such a value numerically.

Fortunately, there is an alternative way of performing the same pruning
exactly that avoids computing the root value. As we will learn in the
following Theorem and Lemma, this will be achieved by studying the link
between the value of the roots of our pieces and their argmax.

###### Theorem 1

Assume $q_{i}$ and $q_{j}$ are two curves that contribute to $Q_{n}(\theta)$
for $\theta>\theta_{0}$, with $i>j$. Then,
$\textsf{root}(q_{i})<\textsf{root}(q_{j})$ if and only if
$a_{i}/(a_{i}+b_{i})<a_{j}/(a_{j}+b_{j})$.

Proof. We begin with some properties of the curves. A curve $q_{i}$ has
$\arg\max q_{i}(\theta)=a_{i}/(a_{i}+b_{i})$. By construction, if $i>j$ then
$a_{i}\leq a_{j}$ and $b_{i}\leq b_{j}$. Finally
$q_{i}(\theta_{0})=q_{j}(\theta_{0})=0$.

We now link the $\arg\max$ of our functions to their roots. First we can
rescale each curve without changing the location of its maximum of its root.
For curve $q_{i}$ define $\tilde{a}_{i}=a_{i}/(a_{i}+b_{i}),\
1-\tilde{a_{i}}=b_{i}/(a_{i}+b_{i})$ and write

$\tilde{q}_{i}(\theta)=\tilde{a}\log\left(\frac{\theta}{\theta_{0}}\right)+(1-\tilde{a})\log\left(\frac{1-\theta}{1-\theta_{0}}\right).$
(8)

We then have that

$\textsf{root}(\tilde{q}_{i})\geq\textsf{root}(\tilde{q}_{j})\iff\tilde{a}_{i}\geq\tilde{a}_{j},$
(9)

This follows by noting that $\tilde{q}_{i}(\theta)$ is unimodal as
$\frac{\partial\tilde{q}_{i}(\theta)}{\partial\theta}$ has one zero, and that
$\forall\
\theta\in[0,1],\frac{\partial\tilde{q}_{i}(\theta)}{\partial\theta}\geq\frac{\partial\tilde{q}_{j}(\theta)}{\partial\theta}\iff\tilde{a}_{i}\geq\tilde{a}_{i}$.

As $q_{i}$ has the same root as $\tilde{q}_{i}$, the result follows
immediately. $\Box$

We will use this result to simplify the pruning procedure in the following
Lemma.

###### Lemma 2

The condition for pruning $q_{i}(l_{i})\leq 0$ at step 2 of Algorithm 2 is
implied by $a_{i}/(a_{i}+b_{i})<a_{i-1}/(a_{i-1}+b_{i-1})$ for $i>1$ and
$a_{i}/(a_{i}+b_{i})<\theta_{0}$ for $i=1$.

Proof. To link the pruning condition and the results from Theorem 1 we need to
show that for $i>1$ whether $q_{i}(l_{i})\leq 0$ we find
$\textsf{root}(q_{i})<\textsf{root}(q_{i-1})$. We start by noting that $l_{i}$
is the root of the function $q^{*}(\cdot)=q_{i}(\cdot)-q_{i-1}(\cdot)$, i.e.
$l_{i}=\textsf{root}(q^{*})$, and that $q^{*}(\theta)\geq 0$ on
$[\theta_{0},l_{i}]$.

Then, if $q_{i}(l_{i})\leq 0$ we have

$q_{i-1}(\theta)=q_{i}(\theta)+q^{*}(\theta)>q_{i}(\theta)\ \text{for}\
\theta\in(\theta_{0},l_{i}),$

and as $\textsf{root}(q_{i-1})\in(\theta_{0},l_{i})$, it follows that

$q_{i}(\textsf{root}(q_{i-1}))<q_{i-1}(\textsf{root}(q_{i-1}))=0.$

By concavity of $q_{i}$ then $\textsf{root}(q_{i})<\textsf{root}(q_{i-1})$.
The argument for $i=1$ follows simply by noting that in this case
$l_{1}=\textsf{root}(q_{1})<\theta_{0}\iff a_{1}/(a_{1}+b_{1})<\theta_{0}$.
$\Box$

Hence we can swap the pruning condition involving the root-finding with the
simpler condition: $a_{i}/(a_{i}+b_{i})<\max\\{\theta_{0},\
a_{i}/(a_{i}+b_{i})\\}$. Given that this novel pruning only involves the
$a_{i},b_{i}$ values, it is no longer necessary to find and store the $l_{i}$
value at each iteration. This result relies, in part, on the functions,
$q_{i}(.)$, being unimodal. Similar results apply for other functions, and
these have been used to generalise the FOCuS algorithm to other one-parameter
exponential family models (see Ward et al.,, 2023).

###  2.4 Extension to the $\theta_{0}$ Unknown Case

We next extend our recursion to the case where both $\theta_{0}$ and
$\theta_{1}$ are unknown. Assume we are observing a stream of Bernoulli random
variables $x_{1},\dots,x_{n}$ distributed as a $Ber(\theta_{0})$ under the
null and as a $Ber(\theta_{1})$ under the alternative. The likelihood ratio
test statistic is

$\displaystyle\mathcal{Q}_{n}$
$\displaystyle=-\max_{\theta\in\mathbb{R}}\sum_{t=1}^{n}h(x_{t},\theta)+$ (10)
$\displaystyle+\underset{\begin{subarray}{c}\tau\in\\{1,\dots,n-1\\}\\\
\theta_{0},\theta_{1}\in\mathbb{R}\end{subarray}}{\max}\left\\{\sum_{t=1}^{\tau}h(x_{t},\theta_{0})+\sum_{t=\tau+1}^{n}h(x_{t},\theta_{1})\right\\},$
(11)

where

$h(x_{t},\theta)=x_{t}\log\theta+(1-x_{t})\log(1-\theta).$

Solving this directly for all possible values of $\tau$, via directly storing
the partial sums, will result in a procedure that has computational complexity
of $\mathcal{O}(n)$ per iteration and $\mathcal{O}(n)$ in memory. We employ a
functional pruning approach and solve (11) exactly through the following
recursion.

###### Proposition 3

Let $Q_{0}(\theta)=0$ and for $n=1,2,\dots$ let $Q_{n}$ be defined by the
recursion:

$Q_{n}(\theta)=\max\left\\{\max_{\theta}\sum_{t=1}^{n}h(x_{t},\theta),\
{Q}_{n-1}(\theta)+h(x_{n},\theta)\right\\}.$

Then, at time $n$,
$\mathcal{Q}_{n}=\max_{\theta}\sum_{t=1}^{n}h(x_{t},\theta)-\max_{\theta}Q_{n}(\theta).$

The proof to Proposition 3 is found in Appendix Proof of Proposition 3.

This recursion can be solved using the same ideas as from Romano et al.,
(2021). As before $Q_{n}(\theta)$ will be the maximum of a set of curves. Each
curve will be of the form

$a\log(\theta)+b\log(1-\theta)+c,$

where the coefficients will depend on the changepoint location associated with
that curve. Importantly the set of changepoint locations that we need to store
curves correspond to locations whose curves would be stored for the
$\theta_{0}$ known case for some $\theta_{0}$. The only difference that
$\theta_{0}$ has on the changepoint locations whose curves are kept is that
when considering positive changes we only need consider
$\theta_{1}>\theta_{0}$, and for negative changes we need only consider
$\theta_{1}<\theta_{0}$. As we need to keep curves that are optimal for some
value of $\theta_{0}$ this means that for positive changes we need to prune
the same curves as for the $\theta_{0}$ known case by for
$\theta_{0}\rightarrow 0$. For negative changes we prunes as for the
$\theta_{0}$ known case but for $\theta_{0}\rightarrow 1$.

The algorithm to implement Ber-FOCuS in the $\theta_{0}$ unknown for
$\theta_{1}>\theta_{0}$ is almost identical to Algorithm 2 except we now store
a 4-tuple $(a_{i},b_{i},c_{i},l_{i})$ for each curve, so in step 9 the new
curve is described by the 4-tuple $(0,0,c,\textsf{root}(q_{k}))$ with
$c=\max_{\theta}\sum_{t=1}^{n}h(x_{t},\theta)$. As before, in practice we can
implement the resulting algorithm without storing $l_{i}$.

###  2.5 Aggregating the Bernoulli Traces and NP-FOCuS

Having covered how to independently monitor the $M$ Bernoulli processes for
all quantiles, we describe how to obtain a global statistic for NP-FOCuS.

Our approach is to consider two aggregation functions. One is to take the
maximum of the statistics, and the other is to take the sum, as in e.g., Mei,
(2010).

NP-FOCuS: Let $Q^{1}_{n}(\theta),\dots,Q^{M}_{n}(\theta)$ be the costs at time
$n$ for the $M$ Bernoulli sequences $\mathbb{I}(y_{n}\leq
p)\,p\in\\{p_{1},\dots,p_{M}\\}$, and let the Ber-FOCuS statistics:

$\mathcal{Q}^{m}_{n}=\max_{\theta}Q^{m}_{n}(\theta)$

Then, we will detect a changepoint at time n whether:

$\sum_{m=1}^{M}\mathcal{Q}^{m}_{n}\geq\xi^{sum}\ \text{or}\
\max_{m\in\\{1,\dots,M\\}}\mathcal{Q}^{m}_{n}\geq\xi^{max},$

with $\xi^{sum},\ \xi^{max}\in\mathbb{R}$. A formal description of the NP-
FOCuS algorithm is reported in Algorithm 3.

Data: $y_{n}\in\mathbb{R}$ the data point at time $n$.

Input: $\xi^{sum},\ \xi^{max},\\{p_{1},\dots,p_{M}\\}$,
$\\{Q^{1}_{n-1}(\theta),\dots,Q^{M}_{n-1}(\theta)\\}$.

1

2for  _$m=1,\dots,M$_ do

$Q^{m}_{n}(\theta),\mathcal{Q}^{m}_{n}\xleftarrow{}\text{Ber-
FOCuS}(x_{n}=\mathbb{I}(y_{n}<p_{m});\ Q^{m}_{n-1}(\theta))$;

// Algorithm 1

3

4 end for

5$\mathcal{S}^{sum}\xleftarrow{}\sum_{m=1}^{M}\mathcal{Q}^{m}_{n}$;

6 $\mathcal{S}^{max}\xleftarrow{}\max_{m\in\\{1,\dots,M\\}}$;

7 if  _$\mathcal{S}^{sum}_{n}\geq\xi^{sum}\ or\
\mathcal{S}^{max}_{n}\geq\xi^{max}$_ then

8 return _$n$ as a stopping point_.

9 end if

10

return _$\\{Q^{1}_{n}(\theta),\dots,Q^{M}_{n}(\theta)\\}$ for the next
iteration._

Algorithm 3 NP-FOCuS (one iteration)

One drawback of having to monitor two streams from two different aggregators
comes with the procedure requiring two thresholds. However, there is an
advantage in monitoring both streams simultaneously, as there are change
scenarios where we expect one to perform better than the other. The sum
statistic should have greater power to detect small shifts in the distribution
that affect all or many quantiles. The maximum statistic will have greater
power to detect larger shifts that affect e.g. just the tail of the
distribution. By setting either $\xi^{sum}$ or $\xi^{max}$ to infinity, our
method would result to considering a test just based on the maximum or the
sum.

###  2.6 Tuning Strategy and Choice of Quantiles

We now describe how to tune the initial parameters of the NP-FOCuS procedure.
Quantile values can be computed either on training observations or on a
probation period, in case of no training data. Rather than taking evenly
spaced quantiles, we build geometrically spaced quantiles following the
approach of Haynes et al., (2017). That is, for fixed $M$, we take
$p_{1},\dots,p_{M}$ in such a way that $p_{m}$ is the empirical quantile with
probability

$\left\\{1+(2n-1)\exp\left[-\frac{(2m-1)}{M}\log(2n-1)\right]\right\\}^{-1}.$

This is to give more importance to values in the tail of the distribution of
the data. In case of a change in the tails of a distribution, for instance,
this would allow for a quicker detection as the traces of the more-extreme
quantiles are more likely to pass the threshold for the max statistic.

As we will find out from Section 3, the choice of $M$ affects the
computational complexity of the procedure. Even if a higher number of
quantiles should lead to better statistical power, in practice there is not
much of a gain in picking values of $M$ greater than 15. This can be checked
from the elbow plot in Figure 5, where we measure the average detection delay
as a function of $M$. Simulations were performed as described in Section 4.1.

Figure 5: Average Detection Delay in function of the number of quantiles M.
Vertical bars represent the standard deviation across the various replicates.

Finally, we tune the thresholds $\xi^{sum}$ or $\xi^{max}$ via a Monte Carlo
approach. This is common practice in the literature when using fixed
thresholds. In particular, we follow the approach of Chen et al., (2022,
Section 4.1). We generate, say, $T$ many sequences as long as the desired
average run length $N$ under the null. This can be achieved either by
simulating the process directly or bootstrapping from training data. We then
compute both the sum and the max statistics on all the $T$ sequences and pick
the max on each. The idea is that the stopping times $S^{sum},S^{max}$, the
times in which respectively the sum or the max statistics go over the
thresholds $\xi^{sum},\ \xi^{max}$, are approximately exponentially
distributed. This means that we want to choose the thresholds so that the
probability of stopping after time $N$, if there is no change, is $\exp(-1)$.

Based on our simulations we pick the estimated thresholds $\tilde{\xi}^{sum},\
\tilde{\xi}^{max}$ to be the smallest values such that the proportion of times
we do not detect a change by time $N$, based on the sum and max statistics
respectively, is less than $\exp(-1)$. This approach would give thresholds
with the approximate average run length we require if we were using only the
max or only the sum statistic. As we are using both, we now fix the ratio of
$\tilde{\xi}^{sum}$ to $\tilde{\xi}^{max}$ but use the same procedure to scale
their values to get a composite test with the correct average run length.

##  3 Computational Complexity of the NP-FOCuS Procedure

The computational and memory complexity per iteration of Algorithm 3 is $M$
times the cost of running Ber-FOCuS. The computational cost of Ber-FOCuS
consists of the pruning step and the update and maximisation steps.

The pruning step has a computational cost that is, on average $O(1)$ per
iteration. This is because at each iteration we will consider some number, $m$
say of curves – and this will result in at least $m-1$ curves being removed.
As we only add one curve per iteration, and a curve can only be removed once,
this limits the average number of curves to be considered per iteration to be
at most 2.

The other steps have a computational cost that is proportional to the number
of curves that are kept. Furthermore, empirically, these step dominate the
computational cost of the algorithm. We can bound the expected cost of these
steps with the following result which bounds the number of curves that are
stored.

###### Theorem 4

Let $x_{1},\dots,x_{n},\dots$ be a realization of an independent Bernoulli
process centred on $\theta_{0}$. Let the number of functions stored by Ber-
FOCuS for $\theta>\theta_{0}$ at iteration $n$ be $\\#\mathcal{I}_{1:n}$. Then
if there is no change prior to $n$

$E(\\#\mathcal{I}_{1:n})\leq(\log(n)+1),$

while if there is one change prior to $n$ then

$E(\\#\mathcal{I}_{1:n})\leq 2(\log(n/2)+1).$

The proof of this theorem is found in Appendix Proof of Theorem 4. By
symmetry, the theorem extends to values of $\theta<\theta_{0}$. From the
theorem, we learn that, at each iteration, we expect to check $k=\log(n)+1$
curves in case a change has not occurred, and $k=2\log(n/2)+1$ curves in the
case where a change has been occurred but is undetected.

##  4 Empirical Evaluation of the Method

###  4.1 Simulation study

We performed a simulation study to assess performances of NP-FOCuS and other
online changepoint procedures in a variety of scenarios, illustrated in Figure
6. Those scenarios were chosen to benchmark procedures over a range of
different challenges in online changepoint detection. Specifically, present
scenarios for Cauchy change-in-scale, Gaussian change-in-mean, change-in-mode
in a mixture of two Gaussian distributions, change-in-mean in an Ornstein-
Uhlenbeck process, decay in a sinusoidal process with random noise, and a
change-in-tails scenario – we cover in details how to generate these in the
Supplementary Materials, Section .2. All the scenarios show i.i.d sequences
with the exception of the OU process and the Sinusoidal process. These two
were added to account for changes in scale and location in presence of strong
temporal dependency, a well-known challenge in changepoint detection which is
present in many real-world applications (see Romano et al.,, 2023; Cho and
Fryzlewicz,, 2020; Hallgren et al.,, 2021).

Figure 6: Example sequences from the six different scenarios considered for
the simulation study. Solid grey lines demarcate a changepoint.

We compare NP-FOCuS with the online non-parametric procedure NUNC (Austin et
al.,, 2023), the method from Ross, (2021), the NEWMA procedure with RFF (see
Section IV of Keriven et al.,, 2020) and the recent online kernel CUSUM
procedure from Wei and Xie, (2022) with Gaussian Kernel (that we call Wei-
CUSUM). As with other methods, Wei-CUSUM procedure assumes independent data.
However, whereas the other methods can be applied to non-independent data by
increasing the threshold for detecting a change, this is impossible for such a
procedure. In fact, in line with other Kernel based online statistics, such as
KCUSUM from Flynn and Yoo, (2019), if the data is serially dependent then the
contribution to the statistic for each new data point will on average by
positive even if there is no change – and thus the method will be prone to
false positives regardless of how large the threshold for the test is chosen.
Thus we exclude Wei-CUSUM from the sinusoidal and OU scenarios. Lastly we add
to the comparison the FOCuS procedure for Gaussian change-in-mean from Romano
et al., (2023). A robust bi-weight loss was employed in this case to account
for the presence of outliers in some scenarios.

Each experiment consisted of 100 replicates. Initial parameters were tuned
according to Section 2.6. Thresholds were selected to achieve an average run
length of $10,000$ observations. Other parameters – such as quantiles, the
value of the bi-weight loss parameter or the NEWMA bandwidth – were obtained
over a probation period consisting of the first $100$ observations (in
Appendix .3 we study the effect of the length of the probation on NP-FOCuS
performances). The NEWMA method involves monitoring the difference of two
exponentially weighted moving averages of features of the data, and it can
take a substantial number of observations for this difference to stabilise. To
account for this we only start monitoring the NEWMA statistic after a burn-in
period of 1000 observations to avoid false-positives whilst its statistics
stabilise. For Wei-CUSUM, we gave oracle knowledge of the pre-change
distribution, and tuned the scale parameter of the Gaussian kernel based on
the variance of the pre-change distribution.

We measure performances in terms of detection delay on sequences generated to
have a change at $\tau=1500$. In Figure 7, we report the proportions of
experiments where a change was detected by time step $t$. Prior to the change,
the lines show the proportion of false positives, while following the change,
they show the number of true detections by a given time step: the perfect
online procedure would achieve a detection within 1 observation from the
change in all sequences. Additionally, we report results in terms of average
detection delay and false positive rate in Tables 1(a) and 1(b).

Figure 7: Proportions of changes detected within $t$ observations following
the change in six different change scenarios. The change is denoted by the
vertical dotted line at $t=1000$.

scenario | FOCuS | Wei-CUSUM | NEWMA | NPFOCuS | NUNC | Ross
---|---|---|---|---|---|---
OUmean |  $>$1000 |  | 225.92 | 87.99 | 995.35 |  $>$1000
cauchy |  $>$1000 | 35.03 | 24.52 | 33.98 | 242.72 |  $>$1000
gauss | 22.08 | 27.2 | 35.71 | 22.26 | 680.64 | 238.87
multim |  $>$1000 | 65.44 | 42.54 | 44.86 | 943.52 |  $>$1000
sinusoidal |  $>$1000 |  | 988.15 | 165.8 | 423.91 |  $>$1000
tails |  $>$1000 | 745.03 |  $>$1000 | 46.97 |  $>$1000 |  $>$1000
(a) Average detection delay. In bold, the best result by row.

scenario | FOCuS | Wei-CUSUM | NEWMA | NPFOCuS | NUNC | Ross
---|---|---|---|---|---|---
OUmean | 0.00 |  | 0.01 | 0.00 | 0.00 | 0.00
cauchy | 0.01 | 0.01 | 0.00 | 0.01 | 0.01 | 0.00
gauss | 0.00 | 0.01 | 0.00 | 0.01 | 0.01 | 0.03
multim | 0.00 | 0.02 | 0.01 | 0.03 | 0.01 | 0.02
sinusoidal | 0.12 |  | 0.00 | 0.03 | 0.13 | 0.00
tails | 0.02 | 0.01 | 0.00 | 0.00 | 0.02 | 0.00
(b) False positive rate.

Table 1: Average detection delay and false positive rate across the various
algorithms for all scenarios of our simulation study.

We learn that, overall NP-FOCuS shows good performances in terms of
statistical power. In the Gaussian case, unsurprisingly, the simple FOCuS with
Gaussian loss has best performances overall, immediately followed by NP-FOCuS.
Additionally, we found that NEWMA achieves faster detections both in the
multimodal scenario, and in the Cauchy scenario, immediately followed by NP-
FOCuS and Wei-CUSUM with similar performances. However, upon further testing,
we noticed that the performance of NEWMA tends to degrade if the change occurs
earlier, even if we use a shorter burn-in period (see Appendix .3 for further
details). Lastly, Kernel based procedures such as Wei-CUSUM, in addition to
requiring the ability to sample from the pre-change distribution, are
sensitive to the choice of the kernel.

We find NP-FOCuS to be more robust to strong dependence in the signal than
other procedures. The most challenging scenario for NP-FOCuS is the sinusoidal
scenario, where we start to consistently estimate the change with a delay of
200 observations. This, as it can be seen from Figure 6, corresponds to the
point where there is a clear difference in scale. In this scenario, for
comparison, other methods are either prone to false positives, or missed
detections.

###  4.2 Monitoring Power Attenuation on Optical Lines

We now evaluate NP-FOCuS on a real-world application, that of monitoring power
attenuation on an optical communication line, a metric that is associated with
the operational performance of the line. We present in Figure 8 four examples
of such a metric. Under normal circumstances, the power attenuation should
show a stationary behaviour over time, as seen in 8a. Any change in
distribution might therefore be of interest to the engineers and be an object
of further investigation. Efficient techniques are required as there is a need
of monitoring a large number of instances over time.

Each time-series has recordings over 80 days. The time-series were first
classified by a domain expert into stationary series, and series with a change
that would want to be flagged to an engineer. Thresholds were tuned over 85
non-changing sequences to achieve a false positive rate of 0.01, as described
in Section 2.6, with quantiles being trained over a probation period of 2
weeks in each sequence.

We draw a comparison of NP-FOCuS stopping times with both Gaussian FOCuS and
NUNC. Threshold and other tuning parameters were adjusted on the same NP-FOCuS
training instances to achieve comparable false positive rates. All 3
algorithms do not report a detection on the example without a change. In the
sequences with a change, NP-FOCuS reported faster detection times. Overall,
Gaussian FOCuS showed the slowest detection delay, even for the example where
the change appears to be a change-in-location. This is attributable to the
fact that the training data is heavy-tailed, and as mentioned in the example
in the introduction this leads to slower detections caused by over-inflated
thresholds. Comparing the nonparametric online procedures, NP-FOCuS and NUNC,
they show similar detection times in the b and d sequences, with NP-FOCuS
resulting in slightly faster detections. However, in line with the previous
numerical study, NUNC struggles in detecting changes in the tails of a
distribution.

Figure 8: Four examples of the measured power attenuation over time. Red
segments at the bottom of the sequences point to detection times of FOCuS for
Gaussian change-in-mean (red), the presented procedure NP-FOCuS (green) and
finally NUNC (blue). Respectively, we show, in a, a sequence with no change,
in b, a sequence with a shift in some modes, in c, a sequence with changes in
the behaviours of the tails, and finally in d an abrupt change-in-location.

##  5 Remarks and Conclusions

We presented a nonparametric approach for online changepoint detection. The
procedure keeps track of the eCDF of a data stream over a set of fixed
quantiles. Keeping track of which observations are above or below each
quantile maps the problem into a Bernoulli change-in-rate problem. We derived
a functional pruning sequential LR test to perform such a task efficiently in
almost linear time.

While we have provided theoretical guarantees about the computational aspects
of our procedure, we have not given any results on the average run length or
the expected detection delay. Obtaining tight results is challenging for our
approach due to additional complexity of allowing for any post change (and
potentially pre-change) parameter, and due to combining of information across
different quantiles. For theoretical results for the related off-line
changepoint procedure, including showing that in the limit where we increase
the number of quantiles, we are able to have power to detect any change in
distribution, see Zou et al., (2014).

The only assumption required by the NP-FOCuS, that of i.i.d. observations, can
lead to a drop in performances when dealing with time-dependent sequences. In
presence of strong dependence, to achieve comparable run lengths to the i.i.d.
case, thresholds are to be inflated, resulting in a slower detection delay. An
approach that can account for the time dependence, such as those of Romano et
al., (2021); Cho and Fryzlewicz, (2020), could be an avenue for further
research. Relying on fixed quantiles also means that the procedure needs some
training observations to understand the range of the data. This can be a
problem in particular when restarting the procedure. If the quantiles are not
known, in fact, those are needed to be re-estimated at every detected change.
This means that no detection will be made for the entire duration of the
probation period. One simple solution would be to store a rolling window of
observations of the length of the probation period. Then, as soon as a
detection is made, the procedure could be restarted with quantiles computed
over those stored observations.

##  6 Acknowledgments

The authors would like to thank Trevor Burbridge at BT for the helpful
conversations that helped motivate the presented procedure. We are grateful to
the Associate Editor and reviewers for their thoughtful comments and
suggestions. The authors acknowledge the financial support of the EPSRC and
Lancaster University grants EP/N031938/1, EP/R004935/1 and Mathematical
Sciences RA (EP/W522612/1).

## References

  * Abramson, (2012) Abramson, J. S. (2012).  Some minorants and majorants of random walks and Lévy processes.  PhD thesis, UC Berkeley.
  * Alvarez-Montoya et al., (2020) Alvarez-Montoya, J., Carvajal-Castrillón, A., and Sierra-Pérez, J. (2020).  In-flight and wireless damage detection in a uav composite wing using fiber optic sensors and strain field pattern recognition.  Mechanical Systems and Signal Processing, 136:106526.
  * Andersen, (1955) Andersen, E. S. (1955).  On the fluctuations of sums of random variables ii.  Mathematica Scandinavica, pages 195–223.
  * Austin et al., (2023) Austin, E., Romano, G., Eckley, I. A., and Fearnhead, P. (2023).  Online non-parametric changepoint detection with application to monitoring operational performance of network devices.  Computational Statistics & Data Analysis, 177:107551.
  * Basseville et al., (1993) Basseville, M., Nikiforov, I. V., et al. (1993).  Detection of abrupt changes: theory and application, volume 104\.  Prentice Hall Englewood Cliffs.
  * Chen, (2019) Chen, H. (2019).  Sequential change-point detection based on nearest neighbors.  The Annals of Statistics, 47(3):1381–1407.
  * Chen and Chu, (2023) Chen, H. and Chu, L. (2023).  Graph-based change-point analysis.  Annual Review of Statistics and Its Application, 10:475–499.
  * Chen et al., (2020) Chen, Y., Wang, T., and Samworth, R. J. (2020).  High-dimensional, multiscale online changepoint detection.  arXiv preprint arXiv:2003.03668.
  * Chen et al., (2022) Chen, Y., Wang, T., and Samworth, R. J. (2022).  High-dimensional, multiscale online changepoint detection.  Journal of the Royal Statistical Society Series B: Statistical Methodology, 84(1):234–266.
  * Cho and Fryzlewicz, (2020) Cho, H. and Fryzlewicz, P. (2020).  Multiple change point detection under serial dependence: Wild energy maximisation and gappy schwarz criterion.  arXiv:2011.13884.
  * Cho and Kirch, (2021) Cho, H. and Kirch, C. (2021).  Data segmentation algorithms: Univariate mean change and beyond.  Econometrics and Statistics.
  * Chu and Chen, (2022) Chu, L. and Chen, H. (2022).  Sequential change-point detection for high-dimensional and non-euclidean data.  IEEE Transactions on Signal Processing, 70:4498–4511.
  * Fearnhead and Fryzlewicz, (2022) Fearnhead, P. and Fryzlewicz, P. (2022).  Detecting a single change-point.  arXiv preprint arXiv:2210.07066.
  * Fearnhead and Rigaill, (2020) Fearnhead, P. and Rigaill, G. (2020).  Relating and comparing methods for detecting changes in mean.  Stat, 9(1):e291.
  * Flynn and Yoo, (2019) Flynn, T. and Yoo, S. (2019).  Change detection with the kernel cumulative sum algorithm.  In 2019 IEEE 58th Conference on Decision and Control (CDC), pages 6092–6099. IEEE.
  * Fridman, (2010) Fridman, P. (2010).  A method of detecting radio transients.  Monthly Notices of the Royal Astronomical Society, 409(2):808–820.
  * Fryzlewicz and Rao, (2014) Fryzlewicz, P. and Rao, S. S. (2014).  Multiple-change-point detection for auto-regressive conditional heteroscedastic processes.  Journal of the Royal Statistical Society: Series B: Statistical Methodology, pages 903–924.
  * Fuschino et al., (2019) Fuschino, F., Campana, R., Labanti, C., Evangelista, Y., Feroci, M., Burderi, L., Fiore, F., Ambrosino, F., Baldazzi, G., Bellutti, P., et al. (2019).  Hermes: An ultra-wide band x and gamma-ray transient monitor on board a nano-satellite constellation.  Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 936:199–203.
  * Hallgren et al., (2021) Hallgren, K. L., Heard, N. A., and Adams, N. M. (2021).  Changepoint detection in non-exchangeable data.  arXiv preprint arXiv:2111.05054.
  * Hawkins and Deng, (2010) Hawkins, D. M. and Deng, Q. (2010).  A nonparametric change-point control chart.  Journal of Quality Technology, 42(2):165–173.
  * Haynes et al., (2017) Haynes, K., Eckley, I. A., and Fearnhead, P. (2017).  Computationally efficient changepoint detection for a range of penalties.  Journal of Computational and Graphical Statistics, 26(1):134–143.
  * Henry et al., (2010) Henry, D., Simani, S., and Patton, R. J. (2010).  Fault detection and diagnosis for aeronautic and aerospace missions.  In Fault tolerant flight control, pages 91–128. Springer.
  * Huang et al., (2014) Huang, S., Kong, Z., and Huang, W. (2014).  High-dimensional process monitoring and change point detection using embedding distributions in reproducing kernel Hilbert space.  IIE Transactions, 46(10):999–1016.
  * Jewell et al., (2020) Jewell, S. W., Hocking, T. D., Fearnhead, P., and Witten, D. M. (2020).  Fast nonconvex deconvolution of calcium imaging data.  Biostatistics, 21(4):709–726.
  * Keriven et al., (2020) Keriven, N., Garreau, D., and Poli, I. (2020).  Newma: A new method for scalable model-free online change-point detection.  IEEE Transactions on Signal Processing, 68:3515–3528.
  * Kovács et al., (2020) Kovács, S., Li, H., Bühlmann, P., and Munk, A. (2020).  Seeded binary segmentation: A general methodology for fast and optimal change point detection.  arXiv preprint arXiv:2002.06633.
  * Kurt et al., (2020) Kurt, M. N., Yılmaz, Y., and Wang, X. (2020).  Real-time nonparametric anomaly detection in high-dimensional settings.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 43(7):2463–2479.
  * Matteson and James, (2014) Matteson, D. S. and James, N. A. (2014).  A nonparametric approach for multiple change point analysis of multivariate data.  Journal of the American Statistical Association, 109(505):334–345.
  * Mei, (2010) Mei, Y. (2010).  Efficient scalable schemes for monitoring a large number of data streams.  Biometrika, 97(2):419–433.
  * Nicolas et al., (2009) Nicolas, P., Leduc, A., Robin, S., Rasmussen, S., Jarmer, H., and Bessières, P. (2009).  Transcriptional landscape estimation from tiling array data using a model of signal shift and drift.  Bioinformatics, 25(18):2341–2347.
  * Page, (1954) Page, E. S. (1954).  Continuous inspection schemes.  Biometrika, 41(1/2):100–115.
  * Pettitt, (1979) Pettitt, A. N. (1979).  A non-parametric approach to the change-point problem.  Journal of the Royal Statistical Society: Series C (Applied Statistics), 28(2):126–135.
  * Pouliezos and Stavrakakis, (2013) Pouliezos, A. and Stavrakakis, G. S. (2013).  Real time fault monitoring of industrial processes, volume 12.  Springer Science & Business Media.
  * Romano et al., (2023) Romano, G., Eckley, I. A., Fearnhead, P., and Rigaill, G. (2023).  Fast online changepoint detection via functional pruning cusum statistics.  Journal of Machine Learning Research, 24(81):1–36.
  * Romano et al., (2021) Romano, G., Rigaill, G., Runge, V., and Fearnhead, P. (2021).  Detecting abrupt changes in the presence of local fluctuations and autocorrelated noise.  Journal of the American Statistical Association, pages 1–16.
  * Ross, (2015) Ross, G. J. (2015).  Parametric and nonparametric sequential change detection in r: The cpm package.  Journal of Statistical Software, 66:1–20.
  * Ross, (2021) Ross, G. J. (2021).  Nonparametric detection of multiple location-scale change points via wild binary segmentation.  arXiv preprint arXiv:2107.01742.
  * Ross and Adams, (2012) Ross, G. J. and Adams, N. M. (2012).  Two nonparametric control charts for detecting arbitrary distribution changes.  Journal of Quality Technology, 44(2):102–116.
  * Scott and Knott, (1974) Scott, A. J. and Knott, M. (1974).  A cluster analysis method for grouping means in the analysis of variance.  Biometrics, 30(3):507–512.
  * Shin et al., (2022) Shin, J., Ramdas, A., and Rinaldo, A. (2022).  E-detectors: a nonparametric framework for online changepoint detection.  arXiv preprint arXiv:2203.03532.
  * Tartakovsky et al., (2005) Tartakovsky, A. G., Rozovskii, B. L., and Shah, K. (2005).  A nonparametric multichart cusum test for rapid intrusion detection.  In Proceedings of Joint Statistical Meetings, volume 7, page 11.
  * Wang and Xie, (2022) Wang, H. and Xie, Y. (2022).  Sequential change-point detection: Computation versus statistical performance.  arXiv preprint arXiv:2210.05181.
  * Ward et al., (2023) Ward, K., Romano, G., Eckley, I., and Fearnhead, P. (2023).  A constant-per-iteration likelihood ratio test for online changepoint detection for exponential family models.  arXiv preprint arXiv:2302.04743.
  * Wei and Xie, (2022) Wei, S. and Xie, Y. (2022).  Online kernel cusum for change-point detection.  arXiv preprint arXiv:2211.15070.
  * Xie et al., (2023) Xie, L., Moustakides, G. V., and Xie, Y. (2023).  Window-limited cusum for sequential change detection.  IEEE Transactions on Information Theory, pages 1–1.
  * Yu et al., (2020) Yu, Y., Padilla, O. H. M., Wang, D., and Rinaldo, A. (2020).  A note on online change point detection.  arXiv preprint arXiv:2006.03283.
  * Zou et al., (2014) Zou, C., Yin, G., Feng, L., Wang, Z., et al. (2014).  Nonparametric maximum likelihood approach to multiple change-point problems.  The Annals of Statistics, 42(3):970–1002.

###  .1 Proofs

### Proof of Proposition 3

For $n=0$ the recursion follows trivially. For $n>0$, at each iteration, we
can either have

$\displaystyle\mathcal{Q}_{n}$
$\displaystyle=\max_{\theta}\sum_{t=1}^{n}h(x_{t},\theta)-\max_{\theta}Q_{n}(\theta)=$
$\displaystyle=\max_{\theta}\sum_{t=1}^{n}h(x_{t},\theta)-\max_{\theta}\sum_{t=1}^{n}h(x_{t},\theta),$

if there is no change up to time $n$, or for a change at time $\tau\in
1,\dots,n-1$:

$\displaystyle\mathcal{Q}_{n}$
$\displaystyle=\max_{\theta}\sum_{t=1}^{n}h(x_{t},\theta)-\max_{\theta}Q_{n}(\theta)=$
$\displaystyle=\max_{\theta}\sum_{t=1}^{n}h(x_{t},\theta)-\max_{\theta}\left[Q_{n-1}(\theta)+h(x_{n},\theta)\right]=$
$\displaystyle=\max_{\theta}\sum_{t=1}^{n}h(x_{t},\theta)-$
$\displaystyle\max_{\theta}\left[\max_{\theta_{0}}\sum_{t=1}^{\tau}h(x_{t},\theta_{0})+\sum_{t=\tau+1}^{n-1}h(x_{t},\theta)+h(x_{n},\theta)\right].$

$\Box$

### Proof of Theorem 4

The proof follows almost exactly the one derived for Theorem 3 in Romano et
al., (2023). While they assume a piecewise constant process plus i.i.d. noise
from a continuous distribution, we denote how this assumption can be relaxed
to include any exchangeable real-valued random process with an additive cost
function, such as our Bernoulli process. Their proof relies on the following
lemmas.

###### Lemma 5

For $i\leq j\leq k$

$\mathcal{I}_{i:k}\quad\subseteq\quad\mathcal{I}_{i:j}\cup\mathcal{I}_{j+1:k}$
(12)

###### Lemma 6

The set of $\tau$ in $\mathcal{I}_{i:j}$ are the extreme points of the largest
convex minorant of the sequence
$S_{i:j}=\\{\sum_{t=i}^{i}x_{t},\dots,\sum_{t=i}^{j}x_{t}\\}$.

We notice that if we write the cost of a segment with a change at
$\tau\in\mathbb{N}$ as

$\begin{split}q_{i:j,\tau}(\theta_{0},\theta_{1})=\sum_{t=i}^{\tau}\left[x_{t}\log\theta_{0}+(1-x_{t})\log(1-\theta_{0})\right]+\\\
\sum_{t=\tau+1}^{j}\left[x_{t}\log\theta_{1}+(1-x_{t})\log(1-\theta_{1})\right],\end{split}$

then for any $\tau,\tau^{\prime}$:

$\begin{split}&q_{i:j,\tau}(\theta_{0},\theta_{1})-q_{i:j,\tau^{\prime}}(\theta_{0},\theta_{1})=\\\
&=\sum_{t=\tau}^{\tau^{\prime}}\left[x_{t}\log\left(\frac{\theta_{1}}{\theta_{0}}\right)+(1-x_{t})\log\left(\frac{1-\theta_{1}}{1-\theta_{0}}\right)\right]\end{split}$

this does not depend on $i$ and $j$. This allows us to extend Lemmas 5 and 6
to the Bernoulli case.

And as both Andersen, (1955) and Abramson, (2012) assume an exchangeable real-
valued random variable, we can adapt Lemma 7 of Romano et al., (2023) to the
Bernoulli case:

###### Lemma 7

Assuming the $x_{t}$ follows an exchangeable real-valued process $\forall t\in
i:j$, then $E(\\#\mathcal{I}_{i:j})=\sum_{1}^{j-i-1}1/(t+1)$.

Having shown that the three lemmas apply to the Bernoulli case we can simply
follow the proof from Romano et al., (2023). $\Box$

###  .2 Detailed Description of Simulation Scenarios

In this section we describe the simulation scenarios, detailing the sampling
for respectively the pre-change and post-change distributions. The Cauchy
change-in-scale was obtained by sampling $y_{1:\tau}$ from a Cauchy with
location 0 and scale 1, i.e., $y_{1:\tau}\sim\textit{Cauchy}(0,1)$ and for
post change $y_{\tau+1:n}\sim Cauchy(0,5)$. Following similar notation, the
Gaussian change-in-mean was obtained by $y_{1:\tau}\sim N(0,1)$ and
$y_{\tau+1:n}\sim N(1,1)$. The multimodal scenario was obtained by sampling
$y_{t}\sim N(0,1)$ with probability $\alpha$ and $y_{t}\sim N(10,1)$ with
probability $(1-\alpha)$, for $\alpha=2/3$, $t=1,\dots,\tau$ and for
$\alpha=1/3$, $t=\tau+1,\dots,n$ after the change.

In the OU scenario, the data is generated by simulating a discrete Ornstein-
Uhlenbeck process (also known as mean-reverting random walk) plus additional
Gaussian noise. This is $y_{t}=\nu_{t}+\epsilon_{t}$, where $\epsilon_{t}\sim
N(0,\sigma_{\epsilon}^{2})$ is a white noise process and $\nu_{t}$ is an
autoregressive process defined by the recursion $\nu_{t}=\nu_{t-1}-\theta
f_{t-1}-\theta\nu_{t-1}+\sigma_{\nu}w_{t-1}$, where $\theta=0.1,\ w_{t}\sim
N(0,\sigma_{\eta}^{2})$. The sequence $f_{t}$ is added to encode a shift in
the mean of the Ornstein-Uhlenbeck process: for $t=1,\dots,\tau$, $f_{t}=0$,
and for $t=\tau+1,\dots,n$, $f_{t}=-10$.

The sinusoidal scenario is generated by sampling $y_{t}\sim N(\mu_{t},1)$ for
$t=1,\dots,n$. The mean $\mu_{t}$ is given by the function $\mu_{t}=\sin(\pi
ft)$, where $f=0.2$ is the frequency of the sinusoidal function. From the
changepoint $\tau$, the mean starts to decrease exponentially with rate
$\lambda=0.005$, i.e. for $t=\tau+1,\dots,n$, $\mu_{t}=A\sin(\pi
ft)\exp(-\lambda t)$.

Finally, the change in tails scenario is generated by sampling $y_{t}\sim
t_{\nu}$ for $t=1,\dots,n$, where $t_{\nu}$ is a Student’s t-distribution with
$\nu$ degrees of freedom. After the changepoint, for $t=\tau+1,\dots,n$, we
randomly selected a $5^{th}$ of the observations and increase by a random
amount drawn from a Poisson distribution with mean 10.

###  .3 Additional Simulation Results

In Section 5 we have commented on the choice of $M$, the number of quantiles.
As for the quantile points, when those cannot be provided directly by the
practitioners, those needs to be estimated through a probation period. We
therefore performed a simulation study to see the effect of the probation
period (amount of training observations) for the quantile estimation on the
performances of NP-FOCuS. As in the main study, we present detection delay and
false positive rate in Tables 2(a) and 2(b). We found that, in general, even
as little as 100 observations are enough for estimating the quantiles. In many
cases the length of the probation period has marginal impact on the
performances. The major differences are seen on the sinusoidal scenario, where
25 and 50 observations are not enough to capture the full domain of the
distribution, resulting in significantly more false positives or slower
detections.

scenario | 25 | 50 | 100 | 150 | 200
---|---|---|---|---|---
OUmean | 102.97 | 94.29 | 87.99 | 78.46 | 80.71
cauchy | 46.19 | 46.45 | 33.98 | 34.47 | 43.82
gauss | 29.57 | 28.3 | 22.26 | 27.4 | 26.58
multim | 71.6 | 61.24 | 44.86 | 55.9 | 60.88
sinusoidal | 154.68 | 252.3 | 165.8 | 140.92 | 130.61
tails | 36.55 | 55.85 | 46.97 | 49.33 | 44.4
(a) Detection delay

scenario | 25 | 50 | 100 | 150 | 200
---|---|---|---|---|---
OUmean | 0 | 0 | 0 | 0.01 | 0.01
cauchy | 0 | 0 | 0.01 | 0.01 | 0
gauss | 0 | 0 | 0.01 | 0 | 0
multim | 0.01 | 0.01 | 0.03 | 0.02 | 0.01
sinusoidal | 0.08 | 0.03 | 0.03 | 0.04 | 0.05
tails | 0.02 | 0 | 0 | 0 | 0
(b) False positive rate

Table 2: Average detection delay and false positive rate for all scenarios of
our simulation study over a range of different probation periods.

To see how the methods deal with an earlier change, we replicate the study of
Section 4.1, but we shift the change at time $\tau=800$. Results are
summarised in Figure 9.

Figure 9: Proportions of changes detected within $t$ observations following
the change in six different change scenarios. The change is denoted by the
vertical dotted line at $t=800$.

We denote how there is a significant drop in performances in NEWMA across all
simulations scenarios. This is particularly evident in the Cauchy and
multimodal scenario. The reduction in detection power is attributable to a
smaller number of random features, in order to stabilise the statistics faster
and achieve shorter burn-in periods. This result is in line with what observed
in Section VI-B of Keriven et al., (2020).