Abstract
Time series, as frequently the case in neuroscience, are rarely stationary, but often exhibit abrupt changes due to attractor transitions or bifurcations in the dynamical systems producing them. A plethora of methods for detecting such change points in time series statistics have been developed over the years, in addition to test criteria to evaluate their significance. Issues to consider when developing change point analysis methods include computational demands, difficulties arising from either limited amount of data or a large number of covariates, and arriving at statistical tests with sufficient power to detect as many changes as contained in potentially highdimensional time series. Here, a general method called Paired Adaptive Regressors for Cumulative Sum is developed for detecting multiple change points in the mean of multivariate time series. The method’s advantages over alternative approaches are demonstrated through a series of simulation experiments. This is followed by a real data application to neural recordings from rat medial prefrontal cortex during learning. Finally, the method’s flexibility to incorporate useful features from stateoftheart change point detection techniques is discussed, along with potential drawbacks and suggestions to remedy them.
Keywords: change point, cumulative sum, adaptive regression splines, nonstationary, bootstrap test, blockpermutation, behaviour, spike counts
1 Introduction
Stationary data are the exception rather than the rule in many areas of science (Aston Kirch, 2012; Elsner ., 2004; Z. Fan ., 2015; Gärtner ., 2017; Latimer ., 2015; Paillard, 1998; Shah ., 2007; Stock Watson, 2014). Time series statistics often change, sometimes abruptly, due to transitions in the underlying system dynamics, adaptive processes or external factors. In neuroscience, both behavioural time series (Durstewitz ., 2010; Powell Redish, 2016; AC. Smith ., 2004) and their neural correlates (Durstewitz ., 2010; Gärtner ., 2017; Latimer ., 2015; Powell Redish, 2016; Roitman Shadlen, 2002) exhibit strongly nonstationary features which relate to important cognitive processes such as learning (Durstewitz ., 2010; Powell Redish, 2016; AC. Smith ., 2004) and perceptual decision making (Hanks Summerfield, 2017; Latimer ., 2015; Roitman Shadlen, 2002). As such, identifying nonstationary features in behavioural and neural time series becomes necessary, both for interpreting the data in relation to the potential influences generating those features, and for removing those features from the data in order to perform statistical analyses that assume stationary observations (Hamilton, 1994; Shumway Stoffer, 2010). Abrupt jumps in time series statistics form one important class of nonstationary events. These are often caused by bifurcations, which, in turn, may occur with gradual changes in parameters of the underlying system (Strogatz, 2001). Consequently, they are of wide interest to both statistical data analysis and the study of dynamical systems, and are commonly referred to as change points (CP; Chen Gupta, 2012).
Detecting CPs has a long and varied history in statistics, and we will not attempt to exhaustively survey the different approaches, including regression models (Brown ., 1975; Quandt, 1958), Bayesian techniques (Chernoff Zacks, 1964) and cumulative sum (CUSUM) statistics (Basseville, 1988; Page, 1954), to name but a few, within the limited scope of this article. Instead, we refer the reader to the excellent reviews on the topic (Aminikhanghahi Cook, 2017; Bhattacharya, 1994; Chen Gupta, 2012) and focus on the offline CUSUM class of methods (Hinkley, 19711) to which PARCS belongs (as opposed to sequential CUSUM methods, Page, 1954, that locate a CP online, while the time series is evolving), specifically methods that aim at detecting CPs in the mean of the time series. CUSUMbased methods are powerful, easy to implement, and are backed up by an extensive literature, theoretical results and various extensions to multiple CPs and multivariate scenarios, making them an ideal starting point. These methods assume that the time series is piecewise stationary in the statistic under consideration (e.g., piecewise constant mean) and rely on a cumulative sum transformation of the time series. Commonly, atmostonechange (AMOC) is identified by maximumtype statistics (Kirch, 2007) at the extremum of the curve resulting from that transformation (Antoch ., 1995; Basseville, 1988).
Extending the CUSUM method to multiple CPs usually involves repetitive partitioning of the time series upon each detection (binary segmentation methods; Bai, 1997; Cho Fryzlewicz, 2015; Fryzlewicz, 2014; Olshen ., 2004; Scott Knott, 1974). This segmentation procedure, however, may hamper detection in later iterations as the reduction in number of observations depletes statistical power exponentially fast as more CPs are to be retrieved. In this article, we develop the PARCS (Paired Adaptive Regressors for Cumulative Sum) method which offers a straightforward extension that leverages the full time series in order to detect multiple CPs, thus providing a new solution to this issue. PARCS rests on the fact that a CUSUM transformation of the data relates to computing an integral transformation of the piecewise constant mean time series model, resulting in a piecewise linear mean function that bends at potential CPs and could be approximated by adaptive regression spline methods (Friedman, 1991; Friedman Silverman, 1989; Stone ., 1997). Namely, rather than attempting to approximate the discontinuous time series mean directly (Efron ., 2004; Vert Bleakley, 2010)
, the PARCS model is an approximation to the continuous CUSUMtransformed time series by a piecewise linear function. The bending points of the PARCS model are each defined by a pair of nonoverlapping piecewise linear regression splines that are first selected by a twostage iterative procedure.
The PARCS model is further refined by a nonparametric CP significance test based on bootstraps (Antoch Hušková, 2001; Dumbgen, 1991; Hušková, 2004; Kirch, 2007; Matteson James, 2014). While analytically derived parametric tests may usually be preferable over bootstrapbased tests due to better convergence and coverage of the tails, in the current CP setting closed form expressions for parametric tests are hard to come by and are usually replaced by approximations (Gombay Horváth, 1996; Horváth, 1997)
. In this case, tests based on bootstraps are preferable since they are known to converge faster to the limit distribution of the test statistic
(often they are also not as conservative as parametric approximations for datasets of a relatively small size; Antoch Hušková, 2001; Csörgö Horváth, 1997; Kirch, 2007). In order to accommodate the possibility of temporally dependent noise in the data (Antoch ., 1997; Horváth, 1997; Picard, 1985), model selection is carried out by a nonparametric blockpermutation bootstrap procedure (Davison Hinkley, 1997; Hušková Slabỳ, 2001; Kirch, 2007)developed specifically for PARCS, which relies on a test statistic that quantifies the amount of bending at each candidate CP. Since model estimation is based on linear regression, PARCS is also effortlessly extended to spatially independent,
multivariate time series.The article is structured as follows. Section 2.1 introduces the CUSUM method for AMOC detection. We then develop the PARCS method, presenting in Section 2.2 the procedure for inferring a nested model that allows for significance testing of multiple CPs, followed in Section 2.3 by an outline of the nonparametric permutation test procedure for refining the PARCS model further. Results in Section 3 illustrate that PARCS improves on several issues inherent in classical methods for change point analysis. In Section 3.1, we compare the PARCS approach to the CUSUM method in detecting a single CP, followed in Section 3.2 by a comparison with standard binary segmentation in detecting multiple CPs. We also demonstrate in Section3.3 that PARCS is successful in detecting CPs in spatially independent, multivariate time series. We then present in Section 3.4 an example from the neurosciences, in which neural and behavioural CPs are compared during operant ruleswitching learning (Durstewitz ., 2010). Finally, we discuss in Section 4 the PARCS approach in relation to other stateoftheart CP detection methods, along with drawbacks and potential extensions.
2 Methods
This section outlines the CUSUM method and the PARCS extension to multiple CPs, in addition to a nonparametric permutation technique to test for the statistical significance of CPs as identified by PARCS. For generality, the formulation assumes temporally dependent observations in the time series, independent observations being a special case.
2.1 CUSUM: Cumulative Sum of Differences to the Mean
A class of methods for identifying a single CP in the mean relies on computing a CUSUM transformation of the time series . A useful formulation that allows for dependent observations in the time series is given by the moving average (MA) step model (Antoch ., 1997; Horváth, 1997; Lombard Hart, 1994; Kirch, 2007),
(1) 
where a jump in the time series mean from baseline to occurs after time step , the change point. The step parameter or weight is positive (negative) when the time series mean increases (decreases) following . The largest integer such that noise coefficient defines a finite order of the MA process, which is 0 for temporally independent observations. We will assume that the MA process is stationary, which will always be the case if it is finite, with
independent and identically distributed (i.i.d.) random variables
(for an infinite process, points for may be considered unobserved, and coefficients have to fulfil certain conditions to make the process stationary, as given, for instance, in Shumway Stoffer, 2010). The Gaussian noise assumption in the MA process can be relaxed, as long as the noise process has zero mean and finite, constant variance
(see Antoch ., 1997; Horváth, 1997; Lombard Hart, 1994; Kirch, 2007, for theoretical results on the more general form of dependent noise). The discrete Heaviside step function, , is defined by,Identifying the presence of a CP requires testing the null hypothesis,
, against the alternative, (Antoch ., 1995; Lombard Hart, 1994). This begins by inferring the time of the step according to a CP locator statistic. A typical offline CP locator statistic is the maximum point of the weighted absolute cumulative sum of differences to the mean (Antoch Hušková, 2001; Horváth, 1997),(2) 
where is the arithmetic mean of the time series (see Figure 1A). The first term on the righthand side corrects for bias toward the centre, where more centrallylocated points are downweighed by an amount controlled by parameter . Other CUSUMbased locator statistics exist with different biascorrecting terms and cumulative sum transformations (Antoch ., 1997; Bhattacharya, 1994; Jirak, 2012; Kirch, 2007). As outlined in the Discussion, PARCS may be modified to include such biascorrecting terms as well. However, as we will demonstrate, PARCS can significantly reduce centre bias even without recourse to such a term. To show this, we will mostly deal with the generic case, , when comparing PARCS to the CUSUM transformation as defined in Eq. 2. This has the added advantage of avoiding having to select an optimal power or an optimal weight factor, a choice that usually depends on prior assumptions on the CP’s potential location (Bhattacharya, 1994). As such, and unless stated otherwise, the term CUSUM transformation will refer, thereof, to the cumulative sum of differences to the mean,
(3) 
where the maximum value,
(4) 
defines a test statistic by which it is decided whether to reject the null hypothesis.
Given potentially dependent observations, , as defined by the model in Eq. 1, nonparametric bootstrap testing proceeds by blockpermutation (Davison Hinkley, 1997; Hušková, 2004; Kirch, 2007), such that temporal dependence in the data is preserved (see Section 3.2). The candidate CP is identified according to Eq. 2 and its associated test statistic is computed by Eq. 4. Estimates and are retrieved from the arithmetic means of before and after using the model in Eq. 1. By subtracting from the time series we arrive at a time series that provides an estimate of the null distribution. The stationary time series is split into blocks of size , chosen such that temporal dependencies are mostly preserved in the permuted time series (Davison Hinkley, 1997). One way to do so is to select the block size to be larger than the order of the underlying MA process, (since the autocorrelation function of an MA process cuts off at order ; Davison Hinkley, 1997). This requires identifying the order which can be determined from the conform time series by inspecting its autocorrelation function (J. Fan Yao, 2003) for different time lags . The autocorrelation function’s asymptotic distribution (Kendall Stuart, 1983),
(5) 
provides a test statistic for deciding the largest time lag at which to reject the null hypothesis , given some preset significance level . The resulting blocks are randomly permuted and each permutation is CUSUMtransformed according to Eq. 3 to compute an conform sample of the test statistic in Eq. 4 (note that we do not know the true step parameter or the true CP , of course, such that this procedure will only yield an estimate of the distribution). A sufficiently large number of permutations results in samples of an conform empirical distribution function (EDF) that weighs every sample equally. The candidate CP is detected when the test statistic as computed from the original time series satisfies , where is a preset significance level and the inverse of the EDF, defined as the largest value out of permutations (Davison Hinkley, 1997; Durstewitz, 2017).
2.2 PARCS: Paired Adaptive Regressors for Cumulative Sum
The PARCS method for estimating multiple CPs rests on the fact that the integral of a piecewise constant function is piecewise linear. The AMOC model as defined in Eq. 1 assumes a piecewise stationary MA process, consisting of two segments with constant mean. A process consisting of segments generalises Eq. 1 to data containing atmostchange,
(6) 
The CUSUM transformation of this process as given by Eq. 3 corresponds to the numerical integration of a piecewise stationary process . That is, is approximately (due to the noise) piecewise linear (exactly piecewise linear in the mean; see Figure 1B). If points at which bends were known, the latter can be fitted by a weighted sum of local piecewise linear basis functions or splines, centred at the knots ,
This fit corresponds to modelling the expected value of , conditioned on spline pair set , resulting in model inference,
which is a simple regression problem that can be solved by estimating the intercept and coefficients that minimise the meansquareerror,
(7) 
However, in the multiple CP detection setting (assuming is known), optimal knot placement is not known a priori, but can be inferred by adaptively adjusting knot locations (Friedman, 1991; Friedman Silverman, 1989; Stone ., 1997) to maximally satisfy the goodnessoffit criterion in Eq. 7. In other words, and as shown in Figures 1C,D, the problem of identifying multiple CPs is replaced by the equivalent problem of inferring the order PARCS model (or PARCS model),
(8) 
with associated tuple that best fits the CUSUM transformation of the time series. Regression coefficients in model 8 are real numbers, while knots in the present time series context are positive integers, excluding the first and last time steps, .
Fitting the PARCS model is based on a forward/backward spline selection strategy (PL. Smith, 1982) with added CP ranking stage and proceeds as outlined in Algorithm 1. Starting with an empty PARCS model, containing only the intercept , a forward sweep increases model complexity to a forward upper bound order by adding at each iteration the spline pair
, not yet contained in the model, that decreases residual meansquareerror the most. A reasonable heuristic for setting
is 2 to 3 times (assuming is known or given some liberal guess). This is followed by a backward pruning iteration, in which the spline pair whose removal increases residual meansquareerror the least is dropped from the model. Pruning removes those knots that were added at the beginning of the forward phase which became redundant as the model was refined by later additions (Friedman Silverman, 1989). This stage continues until the number of knots reaches the preset final upper bound of model complexity , i.e., knots are pruned. Knots are then sorted in descending order according to the amount of explained variance. The ranking iteration returns a nested model by pruning the PARCS model further, down to the PARCS model. The first knot to be pruned, reducing the number of knots to , explains the least variance and is placed last as in the tuple . The last knot to be pruned explains the most variance and is placed first as . Note that regression coefficients are reestimated every time a knot is added to or removed from the PARCS model.The model can be effortlessly extended to the multiple response setting in the case of spatially independent time series (extension to a nondiagonal MA covariance matrix, Stone ., 1997, will be considered elsewhere). Given independent, piecewise stationary MA processes with common CPs ,
(9) 
where , the corresponding multivariate CUSUM transformation is fitted by the multiple response, PARCS model, conditioned on common spline pairs,
using Algorithm 1. Returning CPs that are common to all variables is done by using the goodnessoffit criterion in Eq. 7, averaged over all responses .
2.3 PARCS Model Selection by BlockPermutation Bootstrap
The piecewise linear PARCS formulation, Eq. 8, of the CUSUM transformation in Eq. 3 bends at the CPs. Due to the presence of noise in the original time series , some noise realisations may appear as slight bends in the CUSUMtransformed time series, leading PARCS to return false CPs. As such, the amount of bending at knot can be used as a test statistic for bootstrap significance testing that can refine the PARCS model further. No bending indicates either a constant fit, , or a smooth linear fit, (see also Figure 1C). Thus, a suitable test statistic that quantifies the amount of bending at is given by,
(10) 
where for multivariate time series, the test statistic is the average over all time series.
Before describing the blockpermutation bootstrap method for PARCS, we outline a procedure for identifying the order of the MA noise process, provided as pseudocode in Algorithm 2. First, an conform time series is computed by regressing out the PARCS model of Eq. 8 from the CUSUMtransformed time series and then inverting the CUSUM transformation. This is followed by inspecting the autocorrelation function of for different time lags . The largest time lag at which the null hypothesis is rejected, given some preset significance level is then returned as the order , given some predefined upper bound of MA order, .
Given the tuple CP set returned by Algorithm 1 and an estimate of the dependent normal noise order by Algorithm 2, a blockpermutation bootstrap test returns the subset of significant CPs, as outlined in Algorithm 3. First, an conform time series is computed. For each CP , starting with the one ranked highest, all CPsplines already deemed significant by the bootstrap test are regressed out of . A PARCS model with the remaining knots, including , is estimated and the test statistic , evaluated at according to Eq. 10, is computed. Knot is tested for significance against an conform EDF, estimated through blockpermutation bootstrapping: A total of bootstrap samples is generated from the conform series by randomly permuting blocks of size . For each of these bootstrap samples test statistic is evaluated at knot location , yielding an EDF
which assigns equal probability
to each bootstrapped . A significant is then added to , or rejected as false discovery otherwise. The procedure repeats for the CP next in the rank order. Similar to Algorithm 1, regression coefficients are reestimated every time a knot is added to or removed from the PARCS model.3 Results
We first evaluate the PARCS method on synthetic data in single and multiple CP detection settings, followed by a real data example on detecting behavioural and neural change points during rule learning.
3.1 Alleviating CUSUM Bias in AMOC Detection
We first compare the CUSUM method for detecting a single CP to the PARCS approach in order to evaluate the effect of each method on the centre bias in CP detection. Both white and MA Gaussian noise are considered. We also compare PARCS to the CUSUM locator statistic of Eq. 2 with (the maximum likelihood estimator of CP location under the assumption of i.i.d. Gaussian noise) and identify conditions under which one method is preferable over the other.
Univariate time series of length are simulated according to the step model in Eq. 1 with different levels of white Gaussian noise, , and different ground truth CP locations, . Baseline is set to and step parameter to . Note that in the step model with white Gaussian noise, increasing is equivalent to reducing . A single CP was identified by using the CUSUM method and estimating the PARCS model, both followed by bootstrap significance testing with permutations, nominal significance level , and blocks of size (since noise is independent in this example). Each parameter configuration was repeated for 1000 noise realisations.
We compare bias in CP detection toward the centre of the time series in both the CUSUM and PARCS methods. We measure this centre bias by , which is positive when estimate falls onto the side located toward the centre from , and is negative otherwise. As expected given the choice in Eq. 2, the CUSUM method shows a strong centre bias which increases for lower signaltonoise ratio and more peripheral CPs (see Figure 2A). The CUSUM method’s power decreases for harder parameter settings (higher and more peripheral ) in that true CPs are missed in more of the realisations. The PARCS method results in a significant reduction in centre bias but does not eliminate it completely, and yields more misses relative to CUSUM if both are run at the same nominal level (see Figure 2B). Summary comparison between the two methods is shown in Figure 2C for two exemplary CP locations, .
To fully appreciate the source of CUSUM centre bias and its reduction by PARCS, time series realisations with the two hardest parameter settings ( and ) are considered in Figure 2D, which compares the distribution of in the
of realisations in which both CUSUM and PARCS returned a CP. The histograms show a strongly skewed, heavytailed distribution for CUSUM, compared to a more symmetric distribution around 0 for PARCS, indicating only little bias. Most of the centre bias in PARCS is accounted for by outliers. This is illustrated by excluding outliers in the boxplots, which show a median of 1 time step centre bias in PARCS against median centre bias of 4 time steps in the case of CUSUM. Note that measuring centre bias as defined above does not differentiate between biased detections and false discoveries where, in extreme cases, a CP may be detected beyond the middle point
of the time series, corresponding to centre bias greater than . However, this scenario rarely occurred in the simulation results reported here.While PARCS reduces centre bias, the simulation results above indicate that it behaves more conservatively than CUSUM at the same nominal level. In principle, false rejection rates (type II errors) may be reduced by adjusting the level, at the same time producing more false discoveries. In order to assess how well the nominal significance level agrees with the empirical type I error rate (false discoveries), 1000 white Gaussian noise realisations of length are simulated with . Conclusions drawn from this analysis are largely the same for larger signaltonoise ratio (results not shown). A single CP was extracted using the CUSUM method and estimating the PARCS model on these time series conforming to the null hypothesis . Type I error rates at different nominal levels are shown in Figure 2E as probabilityprobability (PP) plots, depicting the nominal, , against the empirical, , probabilities of accepting the null hypothesis when the null hypothesis is true. While the empirical type I error rate of CUSUM perfectly agrees with the nominal significance level, for PARCS, in contrast, the empirical rate of false discoveries tends toward 0% for and remains smaller than 1% for as large as . This entails that PARCS behaves highly conservatively, and that the nominal level may be adjusted considerably upward without strongly influencing the false discovery rate. On the other hand, despite being more conservative, Figure 2I shows that the receiver operating characteristic (ROC) curve for PARCS, depicting the method’s false discovery rate against its power for different nominal levels, consistently lies above that of CUSUM. For estimating the statistical power of each method, 1000 white Gaussian noise realisations with one CP at a random location in the range (and ) are simulated and type II error rates at different nominal levels are computed. This ROC analysis indicates that for every nominal level for CUSUM, there exists at least one nominal for PARCS such that PARCS has both higher power (fewer type II errors) and lower false detection rate (fewer type I errors), making it the preferable method. This point is explored in more detail later in the context of multiple CP detection in Section 3.2. For shorter time series, PARCS behaves similarly as for longer time series in our simulations, while for CUSUM type I error rates now start to fall below the nominal level as well (Figures 2FH). The area under the ROC curves become smaller for shorter time series for both methods, but the ROC curve of PARCS consistently lies above that of CUSUM in those cases as well (Figures 2JL).
Next, we examine the behaviour of the tests with dependent noise. In the case of temporally dependent noise, an appropriate block size for the bootstrap procedure could be determined by inspecting the autocorrelation function of (see Eq. 5 and Algorithm 2). 1000 noise realisations of length are drawn from an order2 MA process with coefficients and . Since increasing noise variance in the temporally dependent case is not equivalent to decreasing the step parameter, we repeat the analysis with the same parameters from the temporally independent case above but varying the step parameter, and considering two levels of Gaussian noise, .
Figure 3 shows results of the comparison for (top row) and
(bottom row). Similar to the white noise case, the CUSUM method’s centre bias increases for smaller signaltonoise ratio (smaller
), and PARCS, in comparison, consistently reduces centre bias. For the same nominal level, PARCS misses more of the true CPs than CUSUM for peripheral CPs when (top panels of Figures 3A,B), but the type II error rate of the two methods is more comparable in the high noise case, (bottom panels of Figures 3A,B), despite the PARCS method’s more conservative behaviour (far lower type I error rates) in this setting as well. A summary comparison between the two methods is shown in Figure 3C for two exemplary CP locations, . Despite a significant reduction when using PARCS, both centre bias distributions in the 62% of realisations with a CP identified by the two methods, and with the two hardest parameter settings (, and ) remain strongly skewed (bottom panel of Figure 3D).So far, we compared PARCS to the CUSUM statistic with . It is intuitive when developing PARCS to choose for the CUSUM transformation in Eq. 2, since this directly corresponds to the numerical integral of the time series upon which the PARCS approach is based (but see Section 4). Besides, under certain conditions, the CUSUM method using the test statistic with is more sensitive than that with (Antoch ., 1995). However, the CUSUM statistic with returns the maximum likelihood estimator of CP location in an AMOC scenario when noise in the step model of Eq. 1
is i.i.d. and normally distributed, leading theoretically to the strongest centre bias reduction under those conditions
(Antoch Hušková, 2001). We therefore also compare PARCS to this maximum likelihood CUSUM estimator here, henceforth referred to as CUSUM.Univerariate time series of length are simulated according to the step model Eq. 1 with different ground truth CP locations, . We consider only the scenario with largest white Gaussian noise variance, , in this analysis, for which PARCS showed the largest centre bias. A single CP was identified by using the CUSUM method and estimating the PARCS model, both followed by bootstrap significance testing. Other parameters are as in the previous analyses above. CUSUM results in a significant reduction in centre bias compared to PARCS in three of the most peripheral ground truth CPs, , but does not eliminate it completely (see Figure 4A). For the same nominal level, CUSUM also shows lower type II error rates compared to PARCS for these CPs as reported in Table 1, but recall that PARCS has a far lower type I error rate than CUSUM for the same choice of nominal (cf. Figures 2EH and Kirch, 2007). The two methods are comparable in the quality of their detections for all other ground truth CP locations, , with PARCS having a slight advantage.
In order to assess how well the two methods fare in the small sample size limit, and to characterise the convergence behaviour of the bootstrap procedure in each method, we repeat the same analysis for shorter series lengths, . Ground truth CPs are set to the same relative location within the time series as in the simulations. As summarised in Table 1, detection rates deteriorate as series length decreases, as does the bias relative to series length (where the relative location within the series with respect to the periphery is more relevant than the absolute CP location; see Figures 4B,C). Especially for , PARCS performs mostly better than CUSUM, giving higher detection rates (see Table 1) and smaller centre bias (Figure 4C) in the majority of ground truth CPs, although it is still more conservative with near 0% type I error rate (given the bootstrap resolution; see Figure 2G). As we show next, this is a particularly important advantage of PARCS over the CUSUMbased methods when detecting multiple CPs, since CUSUMbased techniques rely on dissecting the time series into smaller segments in this case, reducing sample size at each iteration.
method  type II error rate  

CUSUM  07  03  02  01  01  03  12  
100  
PARCS  17  03  01  01  01  03  16  
CUSUM  30  20  16  13  16  26  37  
50  
PARCS  44  22  12  08  10  19  41  
CUSUM  53  38  35  32  38  44  59  
26  
PARCS  68  41  32  24  29  37  58 
3.2 Detecting Multiple CPs in Univariate Data
For the scenario with multiple CPs, we assess the performance of PARCS in comparison to the CUSUM method with standard binary segmentation (Bai, 1997; Scott Knott, 1974) for univariate data with white Gaussian noise. Standard binary segmentation is known to mislocate CPs in some scenarios, but modifications to the segmentation procedure have been proposed for solving this problem (Fryzlewicz, 2014). We show that PARCS provides an alternative approach. We then discuss a fundamental practical problem in statistical testing when using segmentation methods in general that is avoided by PARCS. Through comparison with standard binary segmentation, we illustrate conditions under which using such methods becomes infeasible. We then consider temporally dependent noise in univariate time series.
The binary segmentation method (Bai, 1997; Scott Knott, 1974) for detecting multiple CPs proceeds as follows (pseudocode can be found in Fryzlewicz, 2014): If, according to a CUSUM test criterion, a CP is detected over the full time series, the series is partitioned at . The procedure is repeated on the resulting left and right segments, potentially returning two additional CPs, and , respectively. The procedure is terminated when no more CPs are detected after subsequent partitioning. Similar to the single CP case, we use bootstrap testing in deciding the significance of a CP at each stage. In the present context, we will refer to this CUSUMbased binary segmentation method simply by ‘CUSUM’.
Three different processes with white Gaussian noise, , of the form given by Eq. 6 and of length are simulated for 1000 realisations each. A baseline and two CPs at time steps and are set in all three scenarios. Weights are set to , and (see insets in top row of Figure 5). For the present comparison, binary segmentation is terminated after at most one partitioning, as this completely suffices to compare the methods (note that this allows CUSUM to detect up to three potential CPs, with only two present in the series). Similarly, the PARCS model is estimated, and both methods use the corresponding permutation bootstrap test with and . In order to compare type I and type II error rates between the two methods, we set nominal levels to 0.05 and 0.30 for CUSUM and PARCS, respectively, which is expected to return about the same factual type I error rates of 5% for both methods in series of length , according to Figures 2EG.
A first look at Figure 5 suggests that both CUSUM and PARCS detect CPs that are close to the ground truth. A more detailed comparison with respect to type I and type II error rates and the quality of CP detections is provided in Table 2. The quality of detections using accuracy scores is defined as the correct detection rate within a range from the ground truth CP location, adjusted for type I errors by an additive term of , where is the factual (empirical) level. This way, the accuracy score is an overall performance measure that takes into account both type I and type II errors, and how far off the detected CP is from the true one.
The first scenario (left panels in Figure 5) has the hardest parameter setting, since is both more peripheral and smaller in magnitude than . CUSUM first detects the easier CP, followed by detecting at the left hand segment as , but with a lower accuracy score as confirmed in Table 2. Similarly, PARCS returns and as the first and second rank CPs, respectively, with accuracy scores higher than those of CUSUM. The relatively low accuracy scores at detecting in both methods are due to its peripheral location and small magnitude. Both type I and type II error rates are markedly lower in PARCS.
method  error rate  accuracy score  

type I  type II  
CUSUM  10 / 14 / 41  08 / 02 / 01  72 / 92 / 77  91 / 78 / 45  
100  
PARCS  02 / 03 / 02  04 / 00 / 01  80 / 96 / 95  96 / 74 / 76  
CUSUM  12 / 18 / 19  27 / 21 / 11  34 / 56 / 61  74 / 38 / 16  
50  
PARCS  04 / 04 / 03  13 / 02 / 05  51 / 82 / 82  85 / 52 / 52  
CUSUM  12 / 22 / 18  41 / 44 / 16  28 / 44 / 71  77 / 26 / 25  
26  
PARCS  06 / 07 / 06  24 / 09 / 10  37 / 75 / 76  79 / 47 / 47 
The second and third scenarios (centre and right panels in Figure 5, respectively) are easier in terms of ground truth parameter settings, since the more peripheral CP has the larger magnitude. In the second scenario (centre panels in Figure 5), the two methods are comparable with regards to their overall accuracy scores as defined above (many detections lie outside the accuracy score range, especially for ), but PARCS has the lower type I and type II error rates. CUSUM has a higher rate of false discoveries than in the first scenario. Its first detection is the higher magnitude CP, which comes with a lower accuracy than PARCS, followed by a detection at the right hand segment with a higher accuracy than PARCS. The relatively low accuracy rates for detecting in both methods are due to its small magnitude. The third scenario (right panels in Figure 5) is an example of a setting in which standard binary segmentation may fail in correctly allocating CPs (see Fryzlewicz, 2014, for a binary segmentation approach that solves this problem). While the performance of PARCS remains about the same as in the second scenario, CUSUM’s first detection diverges from either of the two ground truth CPs in a large number of realisations (see topright panel in Figure 5). The large type I error rate (more than three times that in the first and second scenarios) markedly reduces the accuracy scores, which are substantially lower than for PARCS for both ground truth CPs (see Table 2). Another factor behind the high type I error rates of any iterative procedure including binary segmentation is that the same CP could be detected again at later iterations when an earlier detection is slightly biased.
We now compare the impact of shorter series on the performance of binary segmentation methods and PARCS. We consider 1000 noise realisations from each of the three scenarios with . Two ground truth CPs are set to the same relative location within the time series as for . As seen in Table 2, both type I and type II error rates increase in both methods with the decrease in series length, with the exception of type I error rates for CUSUM in the third scenario. PARCS is consistently the superior method, having both higher statistical power and less false discoveries than CUSUM. While accuracy scores predictably decrease with shorter series length, comparison between the two methods remains qualitatively similar to the case in the first and third scenario, while there is a marked change in the second scenario: While CUSUM is slightly superior in accurately detecting for , PARCS progressively surpasses CUSUM in accuracy as series length decreases. This behaviour is a result of deterioration in the power of the bootstrap test statistic for shorter series. Not only does the overall sample size decrease, but CUSUM in the second scenario with is tasked after a potential first detection of with bootstrapping the CUSUM test statistic on segments as short as 4 or 5 time steps only, which is statistically infeasible, either when using bootstraps or approximate parametric tests (Cho Fryzlewicz, 2015; Fryzlewicz, 2014; Olshen ., 2004). We stress that this is a fundamental drawback to any method that relies on partitioning, and is not specific to standard binary segmentation.
The limitations of binary segmentation methods become more obvious when noise is temporally dependent. For instance, given a time series of length with parameters as in the second scenario and an order2 MA noise process, blocks of size are required for proper blockpermutation. If CUSUM first detected accurately, , the left hand segment would be only 20 time steps long. This allows for only 7 blocks, yielding 5040 possible permutations. This number drops to 720 permutations had been detected only 2 time steps further to the left, which makes it hard to approximate the EDF of the CUSUM statistic reliably. In addition, specifying the block size first requires estimating the MA process order by approximating the conform time series (see Algorithm 2), but potential CPs are not known a priori, due to the recursive nature of binary segmentation methods.
Given these considerations, we focus on PARCS only as we now move over to the case of detecting multiple CPs in series with temporally dependent noise. 1000 noise realisations of length are drawn from an order2 MA process with and coefficients and . Other parameters are as in the previous analysis. In Figure 6, we illustrate distributions of correct detections for the second scenario only (with the other two scenarios qualitatively comparable to their counterparts in Figure 5), but error and accuracy rates on these scenarios are reported as well. After removing the three step changes following PARCS CP detection, the majority of residual time series (70% as shown in bottom panel in Figure 6A) had an autocorrelation that cuts off at the correct order of the ground truth MA noise process, i.e. with being the last coefficient that lies outside the 95% confidence bounds (dashed lines; top panel in Figure 6A; results for the first and third scenarios are comparable). Blockpermutation bootstrap testing with nominal and is carried out on these series with blocks of size 3 and, for other time series, according to the estimated order in Figure 6A (with an upper bound of 10 on block size). Exactly two CPs are detected in more than 99.5% of realisations in all scenarios. Figure 6B shows that the distribution of correct detections for the second scenario is largely concentrated around the ground truth CPs. Accuracy scores in each scenario are, respectively, 96%, 99%, and 99% for and 99%, 89%, and 89% for . Note also the oscillation in the example realisation in Figure 6C, which results from dependent noise with a negative MA coefficient .
3.3 Detecting Multiple CPs in Multivariate Data
The PARCS method’s ability to detect multiple CPs in spatially independent, multivariate time series is demonstrated in Figure 7 on 1000 realisations of length with covariates and white Gaussian noise, . Parameters are set to , , , and . The scaling parameter controls signaltonoise ratio in the time series and is initially set to 1.0. Given these parameter values, the two CPs are not represented in all covariates of the time series, as exemplified in Figure 7A, rendering CPs harder to detect from the averaged univariate time series (with steps differing in sign across the covariates partially cancelling each other, resulting in small weights, and , for the resulting univariate time series).
Following PARCS CP detection, augmented by a bootstrap test with nominal , and , exactly two CPs are detected in 99.9% of realisations. Accuracy scores are 99.8% and 98% for and , respectively. The lower variance in detections, as seen in Figure 7B, is due to the higher average absolute weight compared to .
We then test the method’s performance for smaller signaltonoise ratios with . Figure 7C shows that correct detection rates within a range from the ground truth CPs for different values of remain above 50%, even for magnitudes as small as . These rates are a result of PARCS leveraging CP information from multiple covariates simultaneously, rather than depleting the signal through averaging.
3.4 Detecting Neural Events that Reflect Learning
A previous study by one of the authors and colleagues exemplifies the practical value of change point detection in neuroscience (Durstewitz ., 2010). These authors demonstrated that acquiring a new behavioural rule in rats is accompanied by sudden jumps in behavioural performance, which in turn is reflected in the activity of neural units recorded simultaneously in the medial prefrontal cortex (mPFC). In the current section, we revisit part of these data to showcase PARCS in a real data scenario.
Before moving to the demonstration, it is important to note that the data in question are not normally distributed and potentially include linear trends (Durstewitz ., 2010) not accounted for by the step models in Eqs. 1, 6 and 9. As such, some preprocessing may be necessary for a statistical analysis that is more consistent with the step model assumptions (this may include detrending and potentially some mild smoothing with Gaussian kernels; see Durstewitz ., 2010)
. However, to keep the present demonstration simple, PARCS was applied directly to the data with minimal preprocessing, which only involves squareroottransforming the neural count data for bringing them closer to a Gaussian distribution and stabilising the variance
(Kihlberg ., 1972).In order to show that PARCS can still return reasonable CP estimates under these nonGaussian conditions, we first test its performance on simulated spike count data, before applying it to the empirical data. We simulate 1000 realisations of length with covariates according to a Poisson process. Parameters are set to , , , and . This choice of parameters results in average firing rates that are comparable in their means to the white Gaussian noise case (cf. Figures 7A and 8
A) and to the low firing rates often observed in mPFC neurons. One obvious diversion from Gaussian assumptions in the case of Poisson noise is that the variance is not constant anymore, but is equal to the means within each of the segments separated by true CPs. Following squareroottransforming the data and PARCS
CP detection, augmented by a bootstrap test with nominal , and , exactly two CPs are detected in 92% of realisations. Accuracy scores are 98% and 70% for and , respectively. The lower accuracy scores compared to the white Gaussian noise scenario are due to the lower signaltonoise ratios, resulting from the increase in noise variance with firing rates (cf. Figures 7B and 8B). Nevertheless, these results sufficiently justify the use of PARCS in the present context.We now turn to the experimentally obtained dataset. Six animals were trained on a twochoice deterministic operant rule switching task which proceeds as follows: At the beginning of the session, the animal follows a previously acquired behavioural rule whereby it responds to a visual cue by a lever press for attaining a reward (visual rule). Unknown to the animal, reward contingencies are switched after 20 trials to a novel spatial rule, in which attaining the reward requires pressing a certain baited lever (right or left), regardless of the visual cue. The session is terminated when the animal reaches a preset criterion that indicates that the new rule behaviour has been learnt. In addition to the binary behavioural data of lever presses over trials, spike counts emitted by mPFC units during the 3 seconds following cue onset were collected through single unit recording techniques. Neural and behavioural data from one animal are shown in Figure 9A and 9B, respectively. Trials corresponding to the steady state visual and spatial rule (first and last 20 trials, respectively) are not considered in the analysis.
Time series with one or two significant CPs were described in the original study byDurstewitz . (2010), so PARCS models are estimated for each animal for both the multivariate neural data (multiple response PARCS model; Figure 9A) and the univariate behavioural data (Figure 9B), in addition to PARCS models for the behavioural data as summarised in Figure 9C. As shown in Figure 9B, one neural CP in that exemplary animal matches its behavioural counterpart. The second neural CP, while not as close to its behavioural counterpart, is only 7 trials apart, and the two are highly correlated across animals, as shown in Figure 9C, concurring with the original findings of Durstewitz . (2010). Besides the significant correlation, the corresponding black linear regression line lies very close to the diagonal, indicating that neural and behavioural CPs are not only correlated, but are almost equal. Moreover, those authors report that data from many animals contain at most a single CP (also note the low weight of one of the CPs as estimated from one of the animals using PARCS). Comparing the PARCS behavioural CP to its neural counterpart (blue circles in Figure 9C) shows that correlation remains high and significant. A sample size of is usually recommended for evaluating the significance of a correlation. However, the corresponding linear regression line is also close to the diagonal, in further support for the reliability of this result, and in agreement with the original results despite different procedures: For the neural data in the original study, CUSUMbased detection was performed on a multivariate discrimination statistic defined across the whole neural population, while here, the model was determined directly from the multiple spike count data.
4 Discussion
In the current article, we introduced PARCS, a method for detecting multiple step changes, or CPs, in potentially multivariate, temporally dependent data, supported by a bootstrapbased nonparametric test. We also showed that PARCS substantially reduces centre bias in estimating CPs compared to the most basic specification of the CUSUM method, and presented conditions under which it compares to or outperforms the maximumlikelihood CUSUM statistic. Furthermore, we demonstrated that PARCS may achieve higher sensitivity (statistical power) than CUSUMbased methods while at the same time having lower type I errors in multiple CP scenarios, mainly because PARCS can make use of the full time series while CUSUMbased methods rely on segmenting the time series for detecting multiple CPs. We finally confirmed previous results pertaining to the acquisition of a new behavioural rule and the role of the medial prefrontal cortex in this process.
As already apparent from some of our simulation studies, the basic PARCS method as introduced here leaves room for improvement. In the presence of a single CP, we showed that PARCS strongly reduces the amount of bias toward the centre that results from the direct application of the most basic form of the CUSUM locator statistic. Theoreticallygrounded modifications to the CUSUM transformation that reduce this amount of bias rely on downweighing more centrallylocated points (Kirch, 2007). As shown with PARCS, this problem is not quite as severe. Nevertheless, since PARCS approximates the CUSUM transformation using a regression model, similar downweighing could be incorporated into the PARCS procedure as well by using weighted least squares instead of regular least squares (Hastie ., 2009), which is a straightforward amendment. Furthermore, the PARCS method currently requires a liberal guess of the number of CPs in advance, followed by refinements through nonparametric bootstrap testing. It is desirable, however, especially when no prior information on is available, to have statistical tests as termination criteria for the forward and backward stages. In adaptive regression spline methods (Friedman, 1991; Friedman Silverman, 1989; Stone ., 1997), there is strong empirical evidence (Hinkley, 1969, 19712) backed by theoretical results (Feder, 1975) that the difference in residual meansquareerror between two nested models that differ in one additional knot is well approximated, albeit conservatively, by a scaled
statistic on 4 degrees of freedom
(Friedman, 1991). This led to one nonparametric termination recipe that is based on generalised crossvalidation (Craven Wahba, 1979). Another approach is to infer the piecewise linear regression model with the aid of a parametric test for specifying the number and location of knots, without recourse to iterative procedures (Liu ., 1997). Unfortunately, neither approach is directly applicable to PARCS, since they both require assumptions that are not met in the CUSUMtransformed time series. The CUSUM transformation of the time series is a nonstationary ARMA process. Deriving reasonable generalised crossvalidation (Craven Wahba, 1979; Friedman, 1991; Friedman Silverman, 1989), Fratio (Durstewitz, 2017; Hastie ., 2009) or parametric (Liu ., 1997) test statistics require currently unknown corrections to those tests which account for nonstationarity and the particular form of the ARMA model underlying the CUSUMtransformed data.When multiple CPs are present in the data, PARCS can outperform standard binary segmentation (Bai, 1997; Scott Knott, 1974). Other segmentation methods also solve the problem of mislocating CPs inherent in the standard procedure (Cho Fryzlewicz, 2015; Fryzlewicz, 2014; Olshen ., 2004). Wild binary segmentation (WBS; Fryzlewicz, 2014), for instance, relies on sampling local CUSUM transformations of randomly chosen segments of the time series. The candidate CP with the largest value among sampled CUSUM curves is returned to be tested against a criterion, followed by binary segmentation. WBS is preferable to PARCS in that its test statistic and termination criterion when noise is independent are backed up by rigorous theory, and may be the favourable method when segments are large enough for the test statistic to converge. If series of only limited length are available, however, WBS may run into similar problems as standard binary segmentation for CUSUM, since each detection is still followed by partitioning the data further. WBS also, to the best of our knowledge, currently lacks a thorough analysis on the behaviour of its test statistic for dependent data. It is tempting to speculate on the potential for a hybrid method that capitalises on the desirable features of both methods. Computational demands arise in WBS from the need to choose segment range parameters by sampling few thousand CUSUM curves to which PARCS may offer an easy and efficient workaround: Fryzlewicz (2014) demonstrated that the optimal WBS segment choice is the segment bounded by the two CPs closest to the target CP from each side. PARCS could thus provide an informed selection of boundaries by returning candidate CPs in the data and use these to demarcate segments, rather than random sampling as in WBS.
In dealing with multivariate data, recent methods tackled the computational demands of having a large number of covariates and sparse CP representations (Cho Fryzlewicz, 2015; Wang Samworth, 2018). These methods rely on low dimensional projections of the multivariate CUSUM curve that preserve the CPs and follow this projection by a binary segmentation method. Since PARCS for multivariate time series is also based on the CUSUM transformation, it is straightforward to leverage the computational savings provided by such projection methods in reducing the dimensionality of the PARCS input, while avoiding the drawbacks of binary segmentation methods. This may offer a route for extending PARCS to the important case of multivariate CP detection in mutually dependent time series with spatial dependence, a configuration which these projection methods also consider (Cho Fryzlewicz, 2015; Wang Samworth, 2018). Alternatively, nondiagonal covariance structure in multivariate series may be accounted for by extending the PARCS formulation to the multivariate regression spline realm (Friedman, 1991; Stone ., 1997).
Finally, when analysing the neural and behavioural data during the rule switching task, we mentioned that data may also contain trends that are not accounted for by step change time series models (Durstewitz ., 2010). Caution must be made when analysing real data using CP detection methods in that these methods, PARCS included, assume a step change model underlying the generation of the data and hence may attempt to approximate trends and other nonstationary features by a series of step changes, a point made more explicit by Fryzlewicz (2014) (Durstewitz ., 2010, therefore removed trends around candidate CPs first). Hence, to avoid wrong conclusions with respect to the source and type of nonstationarity in experimental time series, it may be necessary to either augment change point detection by adequate preprocessing (Durstewitz ., 2010) or to generalise time series models for CP detection to include other forms of nonstationarity.
Author Contributions
HT, DD conceived study; HT developed and implemented methods; HT carried out simulations; HT analysed data; HT, DD interpreted results; HT prepared figures; HT wrote manuscript; HT, DD revised and finalised manuscript.
Funding
This research was funded by grants to DD by the German Research Foundation (DFG) (SPP1665, DU 354/82) and through the German Ministry for Education and Research (BMBF) via the e:Med framework (01ZX1311A & 01ZX1314E).
Acknowledgements
The authors thank Dr. Georgia Koppe and Dr. Eleonora Russo for discussions and Dr. Jeremy Seamans for providing the neural and behavioural data.
Data Availability Statement
Method implementation will be freely available on an online repository upon publication.
References
 Aminikhanghahi Cook (2017) Aminikhanghahi2017Aminikhanghahi, S. Cook, DJ. 2017. A survey of methods for time series change point detection A survey of methods for time series change point detection. Knowl. Inf. Syst.512339–367. 10.1007/s101150160987z
 Antoch . (1997) Antoch1997Antoch, J., Hušková, M. Prášková, Z. 1997. Effect of dependence on statistics for determination of change Effect of dependence on statistics for determination of change. J. Stat. Plan. Inference602291–310. 10.1016/S03783758(96)001383
 Antoch . (1995) Antoch1995Antoch, J., Hušková, M. Veraverbeke, N. 1995. Changepoint problem and bootstrap Changepoint problem and bootstrap. J. Nonparametr. Statist.52123–144. 10.1080/10485259508832639
 Antoch Hušková (2001) Antoch2001Antoch, J. Hušková, M. 2001. Permutation tests in change point analysis Permutation tests in change point analysis. Stat. Probab. Lett.53137–46. 10.1016/S01677152(01)000098
 Aston Kirch (2012) Aston2012Aston, JAD. Kirch, C. 2012. Evaluating stationarity via changepoint alternatives with applications to fMRI data Evaluating stationarity via changepoint alternatives with applications to fMRI data. Ann. Appl. Stat.641906–1948. 10.1214/12AOAS565
 Bai (1997) Bai1997Bai, J. 1997. Estimating multiple breaks one at a time Estimating multiple breaks one at a time. Econ. Theory133315–352. 10.1017/S0266466600005831
 Basseville (1988) Basseville1988Basseville, M. 1988. Detecting changes in signals and systems–a survey Detecting changes in signals and systems–a survey. Automatica243309–326. 10.1016/00051098(88)900738
 Bhattacharya (1994) Bhattacharya1994Bhattacharya, PK. 1994. Some aspects of changepoint analysis Some aspects of changepoint analysis. Lect. Notes Monogr. Ser.2328–56. 10.1214/lnms/1215463112
 Brown . (1975) Brown1975Brown, RL., Durbin, J. Evans, JM. 1975. Techniques for testing the constancy of regression relationships over time Techniques for testing the constancy of regression relationships over time. J. R. Stat. Soc. Series B Stat. Methodol.149–192.
 Chen Gupta (2012) Chen2012Chen, J. Gupta, AK. 2012. Parametric statistical change point analysis: With applications to genetics, medicine, and finance Parametric statistical change point analysis: With applications to genetics, medicine, and finance. Basel, SwitzerlandSpringer Science+Business Media, LLC.
 Chernoff Zacks (1964) Chernoff1964Chernoff, H. Zacks, S. 1964. Estimating the current mean of a normal distribution which is subjected to changes in time Estimating the current mean of a normal distribution which is subjected to changes in time. Ann. Math. Stat.353999–1018. 10.1214/aoms/1177700517
 Cho Fryzlewicz (2015) Cho2015Cho, H. Fryzlewicz, P. 2015. Multiplechangepoint detection for high dimensional time series via sparsified binary segmentation Multiplechangepoint detection for high dimensional time series via sparsified binary segmentation. J. R. Stat. Soc. Series B Stat. Methodol.772475–507. 10.1111/rssb.12079
 Craven Wahba (1979) Craven1978Craven, P. Wahba, G. 1979. Smoothing Noisy Data with Spline Functions: Estimating the Correct Degree of Smoothing by the Method of Generalized CrossValidation Smoothing noisy data with spline functions: Estimating the correct degree of smoothing by the method of generalized crossvalidation. Numer. Math.314317–403. 10.1007/BF01404567
 Csörgö Horváth (1997) Csorgo1997Csörgö, M. Horváth, L. 1997. Limit theorems in changepoint analysis Limit theorems in changepoint analysis ( 18). New York, NYJohn Wiley & Sons Inc.
 Davison Hinkley (1997) Davison1997Davison, AC. Hinkley, DV. 1997. Bootstrap methods and their application Bootstrap methods and their application. New York, NYCambridge University Press.
 Dumbgen (1991) Dumbgen1991Dumbgen, L. 1991. The asymptotic behavior of some nonparametric changepoint estimators The asymptotic behavior of some nonparametric changepoint estimators. Ann. Stat.1471–1495. 10.1214/aos/1176348257
 Durstewitz (2017) Durstewitz2017Durstewitz, D. 2017. Advanced Data Analysis in Neuroscience: Integrating Statistical and Computational Models Advanced data analysis in neuroscience: Integrating statistical and computational models. Springer International Publishing.
 Durstewitz . (2010) Durstewitz2010Durstewitz, D., Vittoz, NM., Floresco, SB. Seamans, JK. 2010. Abrupt Transitions between Prefrontal Neural Ensemble States Accompany Behavioral Transitions during Rule Learning Abrupt transitions between prefrontal neural ensemble states accompany behavioral transitions during rule learning. Neuron663438–448. 10.1016/j.neuron.2010.03.029
 Efron . (2004) Efron2004Efron, B., Hastie, T., Johnstone, I. Tibshirani, R. 2004. Least angle regression Least angle regression. Ann. Stat.322407–499. 10.1214/009053604000000067

Elsner . (2004)
Elsner2004Elsner, JB., Niu, X. Jagger, TH.
2004.
Detecting shifts in hurricane rates using a Markov chain Monte Carlo approach Detecting shifts in hurricane rates using a Markov chain Monte Carlo approach.
J. Clim.17132652–2666. 10.1175/15200442(2004)017¡2652:DSIHRU¿2.0.CO;2  J. Fan Yao (2003) Fan2003Fan, J. Yao, Q. 2003. Nonlinear time series: nonparametric and parametric methods Nonlinear time series: nonparametric and parametric methods. New York, NYSpringer.
 Z. Fan . (2015) Fan2015Fan, Z., Dror, RO., Mildorf, TJ., Piana, S. Shaw, DE. 2015. Identifying localized changes in large systems: Changepoint detection for biomolecular simulations Identifying localized changes in large systems: Changepoint detection for biomolecular simulations. Proc. Natl. Acad. Sci. U.S.A.112247454–9. 10.1073/pnas.1415846112
 Feder (1975) Feder1975Feder, PI. 1975. The log likelihood ratio in segmented regression The log likelihood ratio in segmented regression. Ann. Stat.84–97.
 Friedman (1991) Friedman1991Friedman, JH. 1991. Multivariate adaptive regression splines Multivariate adaptive regression splines. Ann. Stat.1911–67. 10.1214/aos/1176347963
 Friedman Silverman (1989) Friedman1989Friedman, JH. Silverman, BW. 1989. Flexible Parsimonious Smoothing and Additive Modeling Flexible parsimonious smoothing and additive modeling. Technometrics3113–21. 10.1080/00401706.1989.10488470
 Fryzlewicz (2014) Fryzlewicz2014Fryzlewicz, P. 2014. Wild binary segmentation for multiple changepoint detection Wild binary segmentation for multiple changepoint detection. Ann. Stat.4262243–2281. 10.1214/14AOS1245
 Gärtner . (2017) Gartner2017Gärtner, M., Duvarci, S., Roeper, J. Schneider, G. 2017. Detecting joint pausiness in parallel spike trains Detecting joint pausiness in parallel spike trains. J. Neurosci. Methods28569–81. 10.1016/j.jneumeth.2017.05.008
 Gombay Horváth (1996) Gombay1996Gombay, E. Horváth, L. 1996. On the rate of approximations for maximum likelihood tests in changepoint models On the rate of approximations for maximum likelihood tests in changepoint models. J. Multivar. Anal.561120–152. 10.1006/jmva.1996.0007
 Hamilton (1994) Hamilton1994Hamilton, JD. 1994. Time series analysis Time series analysis. Princeton, NJPrinceton University Press.
 Hanks Summerfield (2017) Hanks2017Hanks, TD. Summerfield, C. 2017. Perceptual Decision Making in Rodents, Monkeys, and Humans Perceptual decision making in rodents, monkeys, and humans. Neuron93115–31. 10.1016/j.neuron.2016.12.003
 Hastie . (2009) Hastie2009Hastie, T., Tibshirani, R. Friedman, J. 2009. The Elements of Statistical Learning: Data Mining, Inference, and Prediction The elements of statistical learning: Data mining, inference, and prediction (second ). New York, NYSpringer.
 Hinkley (1969) Hinkley1969Hinkley, DV. 1969. Inference about the intersection in twophase regression Inference about the intersection in twophase regression. Biometrika563495–504. 10.1093/biomet/56.3.495
 Hinkley (19711) Hinkley1971aHinkley, DV. 19711. Inference about the changepoint from cumulative sum tests Inference about the changepoint from cumulative sum tests. Biometrika583509–523. 10.1093/biomet/58.3.509
 Hinkley (19712) Hinkley1971bHinkley, DV. 19712. Inference in twophase regression Inference in twophase regression. J. Am. Stat. Ass.66336736–743. 10.1080/01621459.1971.10482337
 Horváth (1997) Horvath1997Horváth, L. 1997. Detection of changes in linear sequences Detection of changes in linear sequences. Ann. Inst. Stat. Math.492271–283. 10.1023/A:1003110912735
 Hušková (2004) Huskova2004Hušková, M. 2004. Permutation principle and bootstrap in change point analysis Permutation principle and bootstrap in change point analysis. L. Horváth B. Szyszkowicz (), Asymptotic Methods in Stochastics: Festschrift for Miklós Csörgő Asymptotic methods in stochastics: Festschrift for Miklós Csörgő ( 44, 273–291). American Mathematical Soc.
 Hušková Slabỳ (2001) Huvskova2001Hušková, M. Slabỳ, A. 2001. Permutation tests for multiple changes Permutation tests for multiple changes. Kybernetika375605–622.
 Jirak (2012) Jirak2012Jirak, M. 2012. Changepoint analysis in increasing dimension Changepoint analysis in increasing dimension. J. Multivar. Anal.111136–159. 10.1016/j.jmva.2012.05.007
 Kendall Stuart (1983) Kendall1983Kendall, MG. Stuart, A. 1983. The Advanced Theory of Statistics The advanced theory of statistics ( 3). London, UKGriffin.
 Kihlberg . (1972) Kihlberg1972Kihlberg, J., Herson, J. Schotz, W. 1972. Square root transformation revisited Square root transformation revisited. Appl. Statist.2076–81. 10.1214/aos/1176343000
 Kirch (2007) Kirch2007Kirch, C. 2007. Block permutation principles for the change analysis of dependent data Block permutation principles for the change analysis of dependent data. J. Stat. Plan. Inference13772453–2474. 10.1016/j.jspi.2006.09.026
 Latimer . (2015) Latimer2015Latimer, KW., Yates, JL., Meister, MLR., Huk, AC. Pillow, JW. 2015. Singletrial spike trains in parietal cortex reveal discrete steps during decisionmaking Singletrial spike trains in parietal cortex reveal discrete steps during decisionmaking. Science3496244184–187. 10.1126/science.aaa4056
 Liu . (1997) Liu1997Liu, J., Wu, S. Zidek, JV. 1997. On segmented multivariate regression On segmented multivariate regression. Stat. Sin.497–525.
 Lombard Hart (1994) Lombard1994Lombard, F. Hart, J. 1994. The analysis of changepoint data with dependent errors The analysis of changepoint data with dependent errors. Lect. Notes Monogr. Ser.23194–209. doi:10.1214/lnms/1215463125
 Matteson James (2014) Matteson2014Matteson, DS. James, NA. 2014. A nonparametric approach for multiple change point analysis of multivariate data A nonparametric approach for multiple change point analysis of multivariate data. J. Am. Stat. Assoc.109505334–345. 10.1080/01621459.2013.849605
 Olshen . (2004) Olshen2004Olshen, AB., Venkatraman, E., Lucito, R. Wigler, M. 2004. Circular binary segmentation for the analysis of arraybased DNA copy number data Circular binary segmentation for the analysis of arraybased DNA copy number data. Biostatistics54557–572. 10.1093/biostatistics/kxh008
 Page (1954) Page1954Page, E. 1954. Continuous inspection schemes Continuous inspection schemes. Biometrika411/2100–115. 10.2307/2333009
 Paillard (1998) Paillard1998Paillard, D. 1998. The timing of Pleistocene glaciations from a simple multiplestate climate model The timing of pleistocene glaciations from a simple multiplestate climate model. Nature3916665378. 10.1038/34891
 Picard (1985) Picard1985Picard, D. 1985. Testing and estimating changepoints in time series Testing and estimating changepoints in time series. Adv. Appl. Probab.174841–867. 10.1017/S0001867800015433
 Powell Redish (2016) Powell2016Powell, NJ. Redish, AD. 2016. Representational changes of latent strategies in rat medial prefrontal cortex precede changes in behaviour Representational changes of latent strategies in rat medial prefrontal cortex precede changes in behaviour. Nat. Commun.712830. 10.1038/ncomms12830
 Quandt (1958) Quandt1958Quandt, RE. 1958. The estimation of the parameters of a linear regression system obeying two separate regimes The estimation of the parameters of a linear regression system obeying two separate regimes. J. Am. Stat. Assoc.53284873–880. 10.1080/01621459.1958.10501484
 Roitman Shadlen (2002) Roitman2002Roitman, JD. Shadlen, MN. 2002. Response of neurons in the lateral intraparietal area during a combined visual discrimination reaction time task Response of neurons in the lateral intraparietal area during a combined visual discrimination reaction time task. J. Neurosci.22219475–9489.

Scott Knott (1974)
Scott1974Scott, AJ. Knott, M.
1974.
A cluster analysis method for grouping means in the analysis of variance A cluster analysis method for grouping means in the analysis of variance.
Biometrics303507–512. 10.2307/2529204  Shah . (2007) Shah2007Shah, SP., Lam, WL., Ng, RT. Murphy, KP. 2007. Modeling recurrent DNA copy number alterations in array CGH data Modeling recurrent DNA copy number alterations in array CGH data. Bioinformatics2313i450–i458. 10.1093/bioinformatics/btm221
 Shumway Stoffer (2010) Shumway2010Shumway, RH. Stoffer, DS. 2010. Time series analysis and its applications: With R examples Time series analysis and its applications: With R examples (3rd ). New York, NYSpringer.
 AC. Smith . (2004) Smith2004Smith, AC., Frank, LM., Wirth, S., Yanike, M., Hu, D., Kubota, Y.Brown, EN. 2004. Dynamic Analysis of Learning in Behavioral Experiments Dynamic analysis of learning in behavioral experiments. J. Neurosci.242447–461.
 PL. Smith (1982) Smith1982Smith, PL. 1982. Curve fitting and modeling with splines using statistical variable selection techniques Curve fitting and modeling with splines using statistical variable selection techniques NASA Report 166034. Langley Research Center, Hampton, VA.
 Stock Watson (2014) Stock2014Stock, JH. Watson, MW. 2014. Estimating turning points using large data sets Estimating turning points using large data sets. J. Econom.178368–381. 10.1016/j.jeconom.2013.08.034
 Stone . (1997) Stone1997Stone, CJ., Hansen, MH., Kooperberg, C. Truong, YK. 1997. Polynomial splines and their tensor products in extended linear modeling: 1994 Wald memorial lecture Polynomial splines and their tensor products in extended linear modeling: 1994 Wald memorial lecture. Ann. Stat.2541371–1470. 10.1214/aos/1031594728
 Strogatz (2001) Strogatz2001Strogatz, SH. 2001. Nonlinear dynamics and chaos: With applications to physics, biology, chemistry, and engineering Nonlinear dynamics and chaos: With applications to physics, biology, chemistry, and engineering. Perseus Books Publishing.
 Vert Bleakley (2010) Vert2010Vert, JP. Bleakley, K. 2010. Fast detection of multiple changepoints shared by many signals using group LARS Fast detection of multiple changepoints shared by many signals using group LARS. JD. Lafferty, CKI. Williams, J. ShaweTaylor, RS. Zemel A. Culotta (), Advances in Neural Information Processing Systems 23 Advances in neural information processing systems 23 ( 2343–2351). Curran Associates, Inc.
 Wang Samworth (2018) Wang2018Wang, T. Samworth, RJ. 2018. High dimensional change point estimation via sparse projection High dimensional change point estimation via sparse projection. J. R. Stat. Soc. Series B Stat. Methodol.80157–83. 10.1111/rssb.12243
Comments
There are no comments yet.