Likelihood Scores for Sparse Signal and Change-Point Detection

Consider a large number $N$ of data streams containing change-points. We consider the situation in which all data up to a given time is available for analysis, so each data stream is an observed sequence of length $T$. At each change-point one or more of the sequences undergo distribution change. The objective is to identify these change-points and the sequences undergoing distribution change. Of interest here is the identification of these change-points when there is sparsity, that is when the number of sequences undergoing change is small compared to $N$. More specifically we want to know the minimum magnitude of change for which the distribution change can be detected under sparsity. And secondly we want to be have an algorithm that is able to detect, with high probability, change-points under the minimum detectable change. See Niu, Hao and Zhang [23] and Wang and Samworth [30] for applications to engineering, genomics and finance.

A typical strategy to deal with sparsity is to subject local scores to thresholding or penalization before summing them up across sequence. Algorithms employing this strategy include the Sparsified Binary Segmentation (SBS) (Cho and Fryzlewicz [8]), the double CUSUM (DC) (Cho [7]), the Informative Sparse Projection (INSPECT) (Wang and Samworth [30]) and the scan algorithm of Enikeeva and Harachaoui [12]. The strategy was also employed by Mei [22], Xie and Siegmund [34] and Wang and Mei [31] in sequential change-point detection on multiple sequences, and Zhang, Siegmund, Ji and Li [37] to detect distribution deviations from known baselines on multiple sequences. Thresholding and penalization suppress noise by removing small and moderate scores, mostly from the majority of sequences without change, thus enhancing the signals from the sparse sequences with changes. It is however unlikely that we are able to specify a threshold or penalization parameter that is optimal at all levels of sparsity.

The higher-criticism (HC) test statistic, proposed by Tukey [29] to check for significantly large number of small p-values, uses multiple thresholds for sparse mixture detection. The number of p-values below a threshold is transformed to a higher-criticism score and this score is maximized over all thresholds. The Berk and Jones [3] test statistic uses multiple thresholds as well but it applies a different p-value transformation. The HC test statistic was shown by Donoho and Jin [9] to be optimal in the detection of a sparse normal mixture. Cai, Jeng and Jin [4] extended it to detect intervals in multiple sequences where the means of a sparse fraction of the sequences deviate from a known baseline and showed that the HC test statistic is optimal. Chan and Walther [5] considered sequence length much larger than number of sequences with detection boundaries that are more complex. They showed that the HC test statistic achieves detection at these boundaries and is optimal in more general settings. They also showed that the Berk-Jones test statistic achieves the same optimality. The HC and Berk-Jones test statistics have the advantage of not requiring a threshold to be specified in advance. However at each threshold we are unable to differentiate a p-value just below the threshold and a p-value that is below the threshold by a large margin. This loss of information results sub-optimality when identifying the exact location where change has occurred.

Our approach here is to convert the p-values into likelihood scores for testing sparse sequences. It can be considered to be a soft form of thresholding in which p-values that are close to zero are penalized less than p-values that are barely significant. Unlike the HC and Berk-Jones test statistic only one transformation of the p-values is needed. It retains the advantage of the HC and Berk-Jones test statistics in not specifying a threshold in advance, but unlike the HC and Berk-Jones test statistics a smaller p-value always has a larger likelihood score.

Since the likelihood scores are transformations of p-values, the proposed method can be applied to any type of distribution changes and it can handle data types that vary across sequences. Our theory however requires a specific distribution family for neat asymptotics and we consider here in particular either normal or Poisson data. We show optimality up to the correct asymptotic constants. For sparse normal change-points these constants are two-dimensional extensions of those in Ingster [17] and Donoho and Jin [9] for sparse normal mixture detection. These constants have been discussed in the context of sparse normal change-point detection assuming a known baseline in Chan and Walther [5] and Chan [6]. For sparse Poisson change-points the constants are new and different from sparse normal constants.

The optimality of multiple sequence identification of change-points up to the correct constant is new. Previous works on optimality for normal data are up to the correct order of magnitude though they go beyond the i.i.d. model, for example Pilliat, Carpentier and Verzelen [26] considered sparse change-point detection in time-series with normal errors. Liu, Guo and Samworth [21] showed optimality up to the best order for normal errors. However they considered a different problem that assumes at most one change-point.

The algorithm we propose here has two steps in the identification of two change-points. The first detection screening step applies the Screening and Ranking (SaRa) idea of Niu and Zhang [24]. The second estimation step for more precise location of change-points uses the CUSUM-like procedure of Wild Binary Segmentation (WBS), cf. Fryzlewicz [14].. This two-step approach saves computation time because the fast screening step evaluates a large number of segments whereas the computationally intensive estimation step is only applied when a change-point has been detected during screening. In contrast for WBS the estimation step is applied on a large number of randomly generated segments. Unlike in Niu and Zhang [24] we do not apply the BIC criterion of Zhang and Siegmund [36] to determine the number of change-points. Instead critical values are specified in advance and binary segmentation, cf. Olshen et al. [25], is applied to detect the change-points sequentially.

An alternative to binary segmentation is estimating the full set of change-points at one go by applying global optimization and making use of dynamic programming to manage the computational complexity. This was employed by the HMM algorithms of Yao [35] and Lai and Xing [20], the multi-scale SMUCE algorithm of Frick, Munk and Sieling [13] and the Bayesian Likelihood algorithm of Du, Kao and Kou [10]. These methods are however designed for single sequence segmentation. Niu, Hao and Zhang [23] provides an excellent background of the historical developments.

The outline of this paper is as follows. In Section II we introduce the sparse likelihood (SL) scores and show that they are optimal in the detection of sparse normal mixtures. In Section III we extend SL scores to detect change-points in multiple sequences. In Section IV we show that SL scores are optimal for change-point detection when the observations are normal or Poisson. In Section V we perform simulation studies on the SL scores. In the appendices we prove the optimality of SL scores.

We start with the simpler problem of detecting a sparse mixture, with the objective of motivating the sparse likelihood score.

Let ${\bf p}=(p^{1},\ldots,p^{N})$ be independent p-values of $N$ null hypotheses and let $p^{(1)}\leq\cdots\leq p^{(N)}$ be the sorted p-values. Tukey (1976) proposed the higher-criticism test statistic

with ${\rm HC}({\bf p})=0$ if $Np^{(n)}>n$ for all $n$, for the overall test that all null hypotheses are true.

Donoho and Jin [9] showed that the HC test statistic is optimal for detecting a sparse fraction of false null hypotheses. Consider test scores $Z^{n}\sim{\rm N}(0,1)$ when the $n$th null hypothesis is true and $Z^{n}\sim{\rm N}(\mu_{N},1)$ for some $\mu_{N}>0$ when the $n$th null hypothesis is false. Define

Donoho and Jin [9] showed that on the sparse mixture $(1-\epsilon){\rm N}(0,1)+\epsilon{\rm N}(\mu_{N},1)$ no algorithm is able to achieve, as $N~{}\rightarrow~{}\infty$,

for testing $H_{0}$: $\epsilon=0$ versus $H_{1}$: $\epsilon=N^{-\beta}$, if $\mu_{N}=\sqrt{2\nu\log N}$ for $\nu<\rho_{Z}(\beta)$. They also showed that the HC test statistic achieves (3) when $\nu>\rho_{Z}(\beta)$ and is thus optimal. Type I error refers to the conclusion of $H_{1}$ when $H_{0}$ is true whereas Type II error refers to the conclusion of $H_{0}$ when $H_{1}$ is true. Ingster (1997, 1998) established the detection lower bound showing that (3) cannot be achieved when $\nu<\rho_{Z}(\beta)$.

Like the HC test statistic, the Berk and Jones [3] test statistic

achieves (3) when $\nu>\rho_{Z}(\beta)$.

We introduce the sparse likelihood scores in Section II.A and show that they achieve (3) in the detection of sparse mixtures, when $\nu>\rho_{Z}(\beta)$, in Section II.B.

Let $X^{n}_{t}$ denote the $t$th observation of the $n$th sequence for $1\leq t\leq T$ and $1\leq n\leq N$. Consider first the model

We are interested in the detection and estimation of

For $s<t$, let $\bar{X}_{st}^{n}=(t-s)^{-1}\sum_{u=s+1}^{t}X_{u}^{n}$. To check for a change of mean on the $n$th sequence at location $t$, select $s<t<u$ and let p-value

In the sparse likelihood algorithm we combine these p-values using $\ell_{N}({\bf p}_{stu})$, where ${\bf p}_{stu}=(p_{stu}^{1},\ldots,p_{stu}^{N})$. When the data follow some other distributions, the corresponding likelihood ratio statistic and p-value can be computed accordingly.

Sparse likelihood scores detects well when only a small fraction of the sequences undergo change of mean. For $T$ large computing the sparse likelihood score for all $(s,t,u)$ is expensive. Instead we combine the approximating set idea of Walther (2010) to first space out the $(s,t,u)$ that are evaluated, and to apply the CUSUM-type scores used in WBS to estimate the change-point location accurately only when the first step indicates a change-point.

In addition to computational savings, through this two-step approach we are able to incorporate multi-scale penalization terms similar to those used in Dümbgen and Spokoiny (2001) and the SMUCE algorithm of Frick, Munk and Sieling (2014), to ensure optimality not only at all levels of sparse change-points, but also at all orders of change magnitudes.

Let $1\leq h_{1}<h_{2}<\cdots$ and $1\leq d_{1}<d_{2}<\cdots$ be integer-valued sequences with $h_{i}\geq d_{i}$ for all $i$. Let $K_{i}(g)=\lfloor\tfrac{g-1}{d_{i}}\rfloor$. Define

The elements of ${\cal A}_{i}(g)$ are the indices where sparse likelihood scores for windows of length $h_{i}$ are computed. Initially we have the full dataset ${\bf X}_{1:T}=(X^{n}_{t}:1\leq t\leq T,1\leq n\leq N)$ and after one or more change-points have been estimated, it is split into sub-datasets ${\bf X}_{b:e}=(X^{n}_{t}:b\leq t\leq e,1\leq n\leq N)$, with length $g=e-b+1$. We check for change-points in ${\bf X}_{b:e}$ using windows specified by ${\cal A}_{i}(g)$.

Let the penalized sparse likelihood scores

For $g\geq 1$, let $i_{g}=\max\{i:h_{i}+d_{i}\leq g\}$. The detection of change-points within ${\bf X}_{b:e}$, with window lengths at least $h_{i_{0}}$, is as follows.

There are two steps in SL-estimate in the estimation of a change-point, when the largest penalized score exceeds the critical value $c$. The first is the identification of an interval $(s(ij),u(ij))$, associated with the largest penalized score, within which a change-point lies. The second is the estimation of the change-point within this interval. In the approximating set ${\cal A}_{i}(g)$, neighboring windows are located $d_{i}$ apart, hence we are unable to estimate the change-points accurately in the first step. Accurate estimation is carried out, with more intensive computations within $(s(ij),u(ij))$, in the second step. Since the second step is performed only after an interval has been identified as containing a change-point, performing this two-step procedure saves computations in regions where scores are generally small and the likelihood of change-points is low.

After a change-point has been identified, we split the dataset into two and execute the same algorithm on each split dataset. To avoid repetitive computations, we start from the window length $h_{i_{0}}$ used in the evaluation of the change-point spltting the dataset, instead of starting from the smallest window length $h_{1}$, on the split datasets. The use of a set of representative set of window-lengths $h_{i}$ for computational savings in change-point detection have been proposed in Willsky and Jones [33]. The recursive segmentation algorithm for the computation of the estimated change-point set $\widehat{\boldsymbol{\tau}}$ is given below, with initialization at $(c,1,1,T,\emptyset)$.

In Figure 2 we show that the critical values of the sparse likelihood algorithm, for a specified Type I error probability, is stable over $N$. Contributing factors include $\ell_{N}(p)$ having a mean that is close to zero and $\ell_{N}({\bf p})$ having exponential tail probabilities not depending on $N$, see (6), when $p$ and $p^{n}$ are uniformly distributed.

Let $\boldsymbol{\mu}=(\mu_{t}^{n}:1\leq t\leq T,1\leq n\leq N)$ and let $J=(\#\boldsymbol{\tau})$ be the number of change-points. We show that the sparse likelihood algorithm is optimal for normal observations in Section IV.A and for Poisson observations in Section IV.B. Consider $T\rightarrow\infty$ as $N\rightarrow\infty$ such that

In Theorems 2 and 4 we specify the detection boundary for asymptotically zero Type I and II error probabilities. Analogous detection boundaries for a single sequence is given in Arias-Castro, Donoho and Huo [1, 2].

In Theorems 3 and 5 we show that Type I and II error probabilities of the sparse likelihood algorithm go to zero at the detection boundary.

Recall that $i_{T}=\max\{i:h_{i}+d_{i}\leq T\}$. Consider the sparse likelihood algorithm with $d_{i}$ and $h_{i}$ satisfying

and critical values $c_{T}$ satisfying

For the sparse likelihood algorithm select parameters $\lambda_{1}>0$ and

We satisfy (11) when $h_{i}\sim\exp(\tfrac{i}{\log i})$ and $d_{i}\sim\tfrac{h_{i}}{i}$ as $i~{}\rightarrow~{}\infty$. Moreover (12) holds because

Condition (11) ensures that the set of $(h_{i},d_{i})$ is sufficiently dense to detect change-points optimally. Condition (12) is required for (13) to hold. The first half of condition (13) ensures Type II error probability goes to 0. The second half ensures Type I error probability goes to 0.

We follow here the simulation set-up in Sections 5.1 and 5.3 of Wang and Samworth [30]. Assume that the random variables are normal with variances that are unknown but equal within sequence. These variances are estimated using median absolute differences of adjacent observations and after normalization, the random variables are treated like unit variance normal.

In the first study there is exactly one change-point $\tau_{1}$. Consider $\mu_{t}^{n}=0$ for $t\leq\tau_{1}$ and all $n$. For $t>\tau_{1}$, let

The objective is to estimate $\tau_{1}$ assuming we know there is exactly one change-point. We estimate $\tau_{1}$ here by

where $\ell^{\rm pen}_{N}$ is the penalized sparse score with $\lambda_{1}=1$ and $\lambda_{2}=\sqrt{\tfrac{\log T}{\log\log T}}$.

We simulate the probabilities that $|\widehat{\tau}_{1}-\tau_{1}|\leq k$ for $k=3$ and 10, and compare against the INSPECT algorithm and the scan algorithm of Enikeeva and Harchaoui [12]. These two algorithms have the best numerical performances in Wang and Samworth [30]. The comparisons in Table I show that the sparse likelihood algorithm performs well.

In the second study there are three change-points within $N=200$ sequences of length $T=2000$, at $\tau_{1}=500$, $\tau_{2}=1000$ and $\tau_{3}=1500$. At each change-point exactly $40$ sequences undergo mean changes. Six scenarios are considered, corresponding to

$1\leq j\leq 3$ and $1\leq n\leq 40$, for $r=0.4,0.6$ and $k=0,20,40$. For $k=0$, the mean changes are within the same 40 sequences at all three change-points, whereas for $k=40$ the mean changes at all three change-points are on distinct sequences. For $k=20$, there is partial overlap of the sequences having mean changes at adjacent change-points. The number of estimated change-points over 100 simulated datasets on each sequence is recorded, as well as the adjusted Rand index (ARI), see Rand [27] and Hubert and Arabie [16], to measure the quality of the change-point estimation.

In the application of the sparse likelihood algorithm, we select $h_{1}=1$ and $h_{i+1}=\lceil 1.1h_{i}\rceil$ for $i\geq 1$, and $d_{i}=\lfloor h_{i}/i\rfloor$, for a total of $i_{T}=61$ window lengths. We select critical value $c_{T}=5$ and parameters $\lambda_{1}=1$, $\lambda_{2}=\sqrt{\tfrac{\log T}{\log\log T}}\doteq 1.94$.

Wang and Samworth [30] showed that INSPECT achieves average ARI of 0.90 when $r=0.6$ and either 0.73 (for $k=20$) or 0.74 (for $k=0$ and 40) when $r=0.4$, comparable to sparse likelihood, see Table II.

In addition to INSPECT, Wang and Samworth [30] considered DC, SBS and scan, as well as the CUSUM aggregration algorithms of Jirak [19] and Horváth and Hušková [15], with average ARI in the range 0.77–0.87 when $r=0.6$ and 0.68–0.72 when $r=0.4$.