Gaetano Romano, Edward Austin, Idris Eckley, and Paul Fernhead
2021-06-25
Introduction to NUNC
NUNC is an acronym for Non-Parametric Unbounded Changepoint detection and is designed to perform the detection of a changes in the distribution of data that has been drawn from an unknown distribution. Furthermore, the algorithm runs in a sliding window. This allows for sequential changepoint detection to take place, along with affording NUNC the ability to deal with non-stationarity in the stream of data. The primary purpose of NUNC is to detect changes in non-stationary or complex distributions of data, where existing parametric detection techniques may fail. The second purpose of NUNC is to perform the analysis online. This means that the data observed by time \(t\) is processed before the next data point arrives, and requires that the memory requirements of the algorithm do not overflow as the algorithm is left running and \(t \rightarrow \infty\).
The NUNC package contains functions that implement both offline sequential changepoint detection from a pre-recorded dataset, and online changepoint detection where the input is a data generating process rather than a fixed set of data. Furthermore, the package contains implementations of all three variants of the NUNC algorithm. The NUNC algorithms are R wrappers for C++ code, providing fast changepoint detection performance. Furthermore, the online implementation of NUNC utilises the Boost library, and a circular buffer structure, to allow for a stream of data observed in real time to be monitored for a change in distribution.
So that NUNC is able to identify changes in distribution for data that has been drawn from an unknown distribution, a non-parametric approach is taken. This means that NUNC specifies the distribution based only on the observed data, and this distribution changes as different data is observed, rather than relying on classical assumptions such as Gaussianity. As will be discussed in the mathematical background Section, NUNC compares the empirical CDF for the pre and post change data and, if the difference is significant enough, then a changepoint is identified in the data. The comparison takes place using a cost function, and this cost function is formulated by aggregating the difference between the two CDFs at K different quantiles.
An example implementation of NUNC for an online data generating process is presented below. Further examples for both NUNC online, and NUNC offline, are presented later.
set.seed(10)onlinedata <-c(rnorm(100, 1, 10), rcauchy(100, 12, .1)) f <-function() { out <-onlinedata[i] i <<-i +1out }i <-1NUNC_run <-NUNC(f, w =25, beta =8, K =6, method ="local")#f is the data generating process, w is the window size, beta is the threshold and K the number of quantilesNUNC_run$changepoint #changepoint location#> [1] 98NUNC_run$t #time detected (edge of sliding window)#> [1] 114
This vignette is divided into three parts. The first explores the mathematical background for NUNC, the second presents the NUNC package and the functions it contains, and the final section uses NUNC to analyse a dataset drawn from an accelerometer.
Mathematical Background for NUNC
NUNC is based upon the comparison of the empirical CDF for independent data points \(x_1, ..., x_t\) drawn from a distribution \(F_{1:t}\). The empirical CDF, \(\hat{F}_{1:t}\), provides an estimate for the distribution of the data, and is given by \[
\hat{F}_{1:t}(q) = \frac{1}{t} \left\{ \sum_{j=1}^{t} \mathbb{I}(x_j < q) + 0.5 \times \mathbb{I}(x_j = q) \right\},
\] for some quantile \(q\). Using the independence assumption, it can be shown that \(t\hat{F}(q) \sim \mathrm{Binom}(t, F(q))\), and a likelihood for a set of data points can then be computed. We write the log-likelihood of the segment \(x_{\tau_1+1}, \ldots, x_{\tau_2}\) as \[
\mathcal{L}(x_{\tau_1+1:\tau_2}; q) = (\tau_2 - \tau_1) \left[\hat{F}_{\tau_1+1:\tau_2}(q) \log(\hat{F}_{\tau_1+1:\tau_2}(q)) - (1-\hat{F}_{\tau_1+1:\tau_2}(q))\log(1-\hat{F}_{\tau_1+1:\tau_2}(q)) \right].
\]
This likelihood can then be used to define a statistic for the detection of a change at a single quantile of the distribution based on the likelihood ratio test. This statistic is defined as \[
\max_{1\leq \tau \leq t} \mathcal{L}(x_{1:\tau}; q) + \mathcal{L}(x_{\tau+1:t}; q) - \mathcal{L}(x_{1:t}; q),
\] and then in order to increase power we aggregate this statistic across \(K\) quantiles. This gives the following test statistic \[
C_K(x_{1:t}) = \max_{1 \leq \tau \leq t} \frac{1}{K} \sum_{k=1}^K 2 \left[ \mathcal{L}(x_{1:\tau}; q_k) + \mathcal{L}(x_{\tau+1:t}; q_k) - \mathcal{L}(x_{1:t}; q_k) \right],
\] and it is this test statistic that is used by NUNC to determine if a change occurs.
NUNC Local is based on this test statistic, although restricts the search to a point inside the sliding window of size W. Furthermore, a point is only identified as a change if the test statistic takes a significant enough value, and this significance is measured by the exceedance of a threshold denoted by \(\beta\). The stopping rule for for a change at time \(t - W + 1 \leq \tau \leq t\) is \[
\max_{t-W+1 \leq \tau \leq t} \sum_{k=1}^K 2 \left[ \mathcal{L}(x_{t-W+1:\tau}; q_k) + \mathcal{L}(x_{\tau+1:t}; q_k) - \mathcal{L}(x_{t-W+1:t}; q_k) \right] \geq K\beta.
\] Using this test statistic, NUNC Local can be seen as identifying changes in the distribution of the data inside the sliding window.
One drawback of NUNC Local is that \(W\) checks are performed inside each sliding window, despite the test statistic being correlated between nearby points and the power of the test being reduced at the edges of the window. To handle this, a subset of the points in the window can be tested instead. This is implemented in the package using the argument in the NUNC functions.
In contrast to NUNC Local, NUNC Global does not search for changes inside the sliding window but instead compares the data in the window to the historic data that has been seen so far. In order to store the historic information in a memory efficient manner we again fix \(K\) quantiles and update the longrun CDF, denoted by \(z^{(t)}_W(\cdot)\), each time a point leaves the sliding window. This update step is denoted by the recursive equation \[
z^{(t+1)}_W(q) = \frac{1}{t-W+1} \bigg [ (t-W) z^{(t)}_W(q) + F_{t-W+1:t-W+1}(q).
\bigg ].
\]
This can then be used to compute the likelihood functions for the data, and a change detection test performed in a similar manner to as defined above.
One question that remains is how many quantiles to use, and what value they should take. In line with Haynes (2017) we propose that \(K = \lceil a \log(W) \rceil\), and select the value of the quantiles from the sample of data using the formula \[
q_k = \hat{F}^{-1} \left(1 + \left(2W+1\right) \exp \left[ \frac{c}{K}(2k-1) \right] \right)^{-1}.
\] The advantage of selecting the quantiles using this formula is that is gives more weight to the tails of the distribution than would be offered by using quantiles that are evenly spaced.
By selecting the \(K\) quantiles using the sample of data, the need to understand the underlying distribution for the data is avoided. Despite this, in some circumstances NUNC may be used to monitor a process with a baseline from a known distribution. In this case, it may be advantageous to specify the quantiles from this known distribution, rather than estimating them from a sample of \(W\) points. The reason for this is that the estimates of the tail quantiles may be poor in situations where \(W\) is not large enough. It is this issue that NUNC Semi-Parametric addresses, by allowing the user to specify the number of quantiles to use through use of the argument.
Having discussed the mathematics underpinning NUNC, we can now explore how the NUNC package can be used to detect changes in the distribution of data. We first present a series of examples that detail the different functions in the NUNC package. After this, we analyse an accelerometer dataset that is contained in the package.
Detecting Changepoints Using NUNC
Three different forms of NUNC are provided in the NUNC package: NUNC Local, NUNC Global, and NUNC Semi-Parametric. NUNC Local tests for changes in the distribution of the data inside the sliding window, whereas NUNC Global tests for changes in the distribution of the data by asking whether the observations in the window are drawn from the same distribution as the historic observations. In both these methods, the test uses quantiles of the empirical CDF that have been selected using the observed data. In some cases, however, it is argued that there will be a strong belief as to what the distribution of the data is, and so NUNC Semi-Parametric allows for the user to specify the quantiles that are used.
In this section we will explore the different variants, and present code for implementing the algorithms on several different examples. We begin by examining the offline variants of NUNC. For NUNC Local we consider a change in a set of mixture distributions:
set.seed(1)x_pre <-c(rnorm(500), rnorm(200, 1, 4), rnorm(300, 2, 2))x_pre <-sample(x_pre, 1000)x_post <-c(rnorm(100), rnorm(200, -3, 1), rnorm(200, 3, 2), rnorm(500, -1, 1))x_post <-sample(x_pre, 1000)x <-c(x_pre, x_post)#we use [1:3] to supress the data outputNUNCoffline(data = x, w =200, beta =12, K =4*log(200), method ="local")[1:3]#> $t#> [1] 1181#> #> $changepoint#> [1] 1064#> #> $method#> [1] "local"
We see that NUNCoffline has returned three outputs: , the time the change is detected; , the location of the change in the window; and , which is the method used. We see that a similar analysis can be obtained in a shorter period of time using the argument. This is specified as an additional argument using the input: this is a containing the argument and can be written as follows:
A <-Sys.time()NUNCoffline(data = x, w =200, beta =12, K =4*log(200), method ="local")$changepoint#> [1] 1064B <-Sys.time()NUNCoffline(data = x, w =200, beta =12, K =4*log(200), method ="local", params =list(max_checks =30))$changepoint#> [1] 1069C <-Sys.time()paste("NUNC Local has a", B-A)#> [1] "NUNC Local has a 2.38081789016724"paste("NUNC Local on a subset of points in the window has a ", C-B)#> [1] "NUNC Local on a subset of points in the window has a 0.420536994934082"
We next consider analysing a set of data using NUNC Global. We illustrate this using an example of data that is heavy tailed:
We see that the input structure is similar, however the output structure also contains an argument called ; this is the value of the longrun CDF at the \(K\) quantiles for the data.
The third form of NUNC is NUNC Semi-Parametric. We explore how this can provide better results than NUNC Local in cases where the quantiles for the data are pre-specified correctly by examining a scenario where there is a high variance in the distribution for the data under the null hypothesis. This scenario makes estimation of the quantiles from a small sample of data difficult.
As in the case, the additional argument is specified using the argument. We see in the following example we obtain a correct changepoint when using NUNC Semi-Parametric, whereas an incorrect changepoint is returned using NUNC Local.
Online analysis of a data stream can also be performed using NUNC, and this is illustrated by the following examples using the function. This takes a similar set of inputs to the function, however rather than take a argument, we instead use {f =}, where is a data generating function that returns a single number at each time point.
set.seed(4) onlinedata <-c(rnorm(300, 1, 10), rcauchy(1200, 12, .1)) f <-function() { out <-onlinedata[i] i <<-i +1#uncomment and run with x11() or windows() to see real time plotting#plot(tail(onlinedata[1:i], 80), ylab = "test data", type = "l")out }i <-1NUNC(f, 50, beta =10, K =15, params =list(max_checks =3), method ="local")#> $t#> [1] 324#> #> $changepoint#> [1] 299#> #> $method#> [1] "local"#> #> attr(,"class")#> [1] "NUNCout"i <-1NUNC(f, 50, beta =10, K =15)#> $t#> [1] 345#> #> $changepoint#> [1] 294#> #> $zVals#> [1] 0.97457627 0.03050847 0.03728814 0.05423729 0.11186441 0.20677966#> [7] 0.30847458 0.42711864 0.57288136 0.74576271 0.84745763 0.90847458#> [13] 0.92881356 0.95254237 0.96949153#> #> $method#> [1] "global"#> #> attr(,"class")#> [1] "NUNCout"
In the above example, the data is fed into the function one point at a time and analysis performed. If the code is run with the plotting aspect of the function uncommented, a user will also be able to see the stream of data plotted whilst it is being observed.
Analysing Accelerometer Data Using NUNC
Now that a user is familiar with the different functions contained in the NUNC package, we can explore a set of data created using an accelerometer. This dataset is contained within the package, and contains 100 instances where the movement of a Dualshock4 controller has been recorded. The movement is captured by a change in the y axis value recorded by the motion sensor in the controller.
data(accelerometer)
Each entry in the list is a data frame containing three variables: The four actions that can take place are: picking up the controller, sitting with the controller and moving it, sliding the controller along a surface, and shaking the controller. In every example, the movement of the controller changes the distribution of y axis values.
We can explore an example of this change in distribution:
We can now use NUNC to detect changes in the distribution of the data. For example
accelerometer[[52]]$changepoint#> [1] 497NUNCoffline(data = accelerometer[[52]]$y, w =100, beta =17, K =10, method ="local",params =list(max_checks =30))$changepoint#> [1] 573NUNCoffline(data = accelerometer[[52]]$y, w =100, beta =17, K =10, method="global")$changepoint#> [1] 535
In both instances of NUNC, there is some delay between the emergence and the detection of the change. This is to be expected given the sequential nature of the algorithm. We also observe the better performance of NUNC Global, which again is to be expected given the fact that the pre-change distribution is stationary.
We are then able to plot the output of this analysis using the function, with the changepoint location indicated by the red line:
We can build upon this and perform a more detailed analysis of the accelerometer dataset using NUNC. In particular, we can assess both the detection power, and the detection delay, for both NUNC Global and NUNC Local. Furthermore, the detection delay is given by for NUNC Local and NUNC Global respectively, along with the standard deviation of this delay. For both variants of NUNC, we use the threshold given by a theoretical result presented in Austin (2021). This is set to control the probability of false alarm under the null hypothesis at the 10% level.
The results for this experiment are given in the table below.
Table showing results from NUNC analysis of accelerometer dataset
Power
Delay
Delay Sd
Local
0.72
69.65
136.38
Local Grid
0.72
74.53
142.86
Global
0.74
83.31
94.64
The code to implement this example is given below:
NUNC_local_detected <-NUNC_grid_detected <-NUNC_global_detected <-0NUNC_local_delay <-NUNC_grid_delay <-NUNC_global_delay <-rep(NA, 100)#compute theoretical thresholds for the test to control false alarms at 10%local_threshold <-find_beta(alpha =0.1, t =2000, w =100, K =15, method ="local")global_threshold <-find_beta(alpha =0.1, t =2000, w =100, K =15, method ="global")for(i in1:100){ local_run <-NUNCoffline(accelerometer[[i]]$y, 100, local_threshold, 15, "local") grid_run <-NUNCoffline(accelerometer[[i]]$y, 100, local_threshold, 15, "local",params =list(max_checks =30)) global_run <-NUNCoffline(accelerometer[[i]]$y, 100, local_threshold, 15, "global")#if the changepoint is detected, and is not before the accelerometer is moved,#record a successful detectionif(local_run$changepoint !=-1&local_run$t >=accelerometer[[i]]$changepoint){ NUNC_local_detected <-NUNC_local_detected +1 NUNC_local_delay[i] <-local_run$t -accelerometer[[i]]$changepoint } if(grid_run$changepoint !=-1&grid_run$t >=accelerometer[[i]]$changepoint){ NUNC_grid_detected <-NUNC_grid_detected +1 NUNC_grid_delay[i] <-grid_run$t -accelerometer[[i]]$changepoint } if(global_run$changepoint !=-1&global_run$t >=accelerometer[[i]]$changepoint){ NUNC_global_detected <-NUNC_global_detected +1 NUNC_global_delay[i] <-global_run$t -accelerometer[[i]]$changepoint }print(i)}#we can now compute the power and average detection delaylocal_power <-NUNC_local_detected /100grid_power <-NUNC_grid_detected /100global_power <-NUNC_global_detected /100local_delay <-mean(NUNC_local_delay, na.rm =TRUE)local_delay_sd <-sd(NUNC_local_delay, na.rm =TRUE)grid_delay <-mean(NUNC_grid_delay, na.rm =TRUE)grid_delay_sd <-sd(NUNC_grid_delay, na.rm=TRUE)global_delay <-mean(NUNC_global_delay, na.rm =TRUE)global_delay_sd <-sd(NUNC_global_delay, na.rm =TRUE)
