# Nonparametric changepoint detection for series with variable number of measurements across time

I have been looking at a lot of recent changepoint detection algorithms ( *-PELT, NEWMA, ...) but it seems they all work on a single (or multiple) variable(s) that are composed of each a single value for each date.

My problem is a bit different as I have a variable amount of values per "date" (could even be represented using CDF or KDE) and I'd like to detect changes in behavior of those values. (For example changes in mean, standard deviation, shape, etc).

So instead of having series of single values, for example:

x0 = 0.1
x1 = 0.5
x2 = 0.3
x3 = 0.4
x4 = 2.5
x5 = 2.1
x6 = 2.3


I instead have series of multiple values (count per "date" can change), for example:

x0 = (0.1,0.11,0.45,0.26,...)
x1 = (0.5,0.3,0.4,0.43,...)
x2 = (0.3,0.2)
x3 = (0.4,0.21,0.32,0.54)
x4 = (2.5,2.1,2.65,2.57,...)
x5 = (2.1,2.15,2.6,2.33, 2.41)
x6 = (2.3, 2.12, 2.39, 2.54, 2.16)


I had a few ideas but that I don't like very much:

• Computing a descriptive statistic (mean, median, stddev) for each date, and apply changepoint detections to those
• This can get quite expensive
• This doesn't seem reliable
• assign each value of a "date" to multiple fake "dates"
• Can and will skew the results
• There is a big loss of information

Is there some algorithm that could could work with such data?

Edit: http://www.jmlr.org/papers/v20/16-155.html Could be answering this question, still have to read it.

• What does a "changepoint" look like in this context? In a univariate or regular multivariate problem what is changing is clear. Are the smaller dimensional dates simply missing values in some of the dimensions, i.e. x2=(NA,0.3,NA,NA,0.2,NA,...)? – adunaic Jul 29 '20 at 10:14
• No the values are measures of the same system, I just don't have the same amount of those per date. What I'm trying to detect is if there is a change in distribution of those values. – Lectem Jul 30 '20 at 9:50
• So is it that there are several measurements per day but taken at different times in the day and you have just stacked them by day? If you have the individual times within the day then you should create a univariate time series using that. If you don't have the individual times but know that the order is the same order they were collected then you can analyze them in that order without the dates, without loss of power. The only time this becomes a problem is if you are wanting to fit seasonal or auto-correlation in your changepoint model. – adunaic Jul 31 '20 at 12:42
• Well it's kind of a weird spot I think as my date isn't really a date but a version of the system being measured, and I could have multiple measurements of the same version happening at the same time or overlapping, but each measurement has multiple values. It's not easy to describe but basically my X axis is not the real time but the version/configuration of what is being measured. In theory I could even have a measurement done on day 1 for version A, another on day2 for version B, and then another one on day3 for version A. (each measurement yielding multiple values) – Lectem Jul 31 '20 at 21:31
• Take a look at the approach of Adams & MacKay: arxiv.org/abs/0710.3742 it is easily generalizable to data in any kind of space. And also the approach of Scargle: ui.adsabs.harvard.edu/abs/1998ApJ...504..405S/abstract – pglpm Aug 5 '20 at 8:10

If you are open to using R, here is a solution using mcp. mcp can infer the location of changes in means (worked examples), variances (worked example), autocorrelation (worked example), and any combination of these.

Set up data:

x0 = c(0.1,0.11,0.45,0.26)
x1 = c(0.5,0.3,0.4,0.43)
x2 = c(0.3,0.2)
x3 = c(0.4,0.21,0.32,0.54)
x4 = c(2.5,2.1,2.65,2.57)
x5 = c(2.1,2.15,2.6,2.33, 2.41)
x6 = c(2.3, 2.12, 2.39, 2.54, 2.16)

df = data.frame(
x = c(rep(0, length(x0)), rep(1, length(x1)), rep(2, length(x2)), rep(3, length(x3)), rep(4, length(x4)), rep(5, length(x5)), rep(6, length(x6))),
y = c(x0, x1, x2, x3, x4, x5, x6)
)


We model this as a single change in intercept and a change in variance:

model = list(
y ~ 1 + sigma(1),  # intercept and variance
~ 1 + sigma(1)  # new intercept and new variance
)
fit = mcp(model, df, par_x = "x")


Here are some plots of the fitted means (left) and variances (right):

plot(fit) + plot(fit, which_y = "sigma")


Here are the parameter estimates:

> summary(fit)

# Population-level parameters:
#     name mean lower upper Rhat n.eff
#     cp_1 3.50 3.053  4.00    1  6433
#    int_1 0.32 0.239  0.41    1  5015
#    int_2 2.35 2.229  2.47    1  5276
#  sigma_1 0.15 0.095  0.22    1  3067
#  sigma_2 0.22 0.139  0.32    1  3171


cp_1 is the change point, int_k is the intercept in segment $$k$$ and sigma_k is the variance in segment $$k$$. In this case, the change in intercept clearly is very informative while the change in variance is less pronounced.

You can also model continuously increasing variances (y ~ sigma(1 + x), etc.