# Bayesian change point detection

Really naive question. I have a time series. I know how to perform segmentation (like binary segmentation algorithm). The goal is to find intervals generated from different probabilistic models.

But I have all information about possible models (distribution shape, variance, mean). So for each time point I have likelihood of each model and its prior. => I can calculate posterior for each time point, any model and any interval.

Problem: if I just segment the time series using maximum posterior probability, I will have too many change points. HMM can be a solution, but it also takes only one point into account and does not "look" at the whole interval. Also it is difficult to apply HMM for non-normal data.

It can be solved with sliding window, but it is not clear how to choose size of sliding window.

Is there an algorithm for this type of Bayesian change point detection (when you know possible models)? Like HMM, but takes the interval into account and can work with any parametric distribution? Heuristic algorithm is good too.

How can I apply maximum likelihood clustering for this problem?

UPD: Simulation of the problem:

variances <- runif(1000,0.01,0.5)

coverages <- c()

for (i in seq(1:100)) {
coverages <- c(coverages, rnorm(1, mean=0, sd=variances[i]))
}

for (i in seq(101:200)) {
coverages <- c(coverages, rnorm(1, mean=-log(2), sd=variances[i] / 0.75))
}

for (i in seq(201:300)) {
coverages <- c(coverages, rnorm(1, mean=log(3/2), sd=variances[i] * 0.75))
}

for (i in seq(301:1000)) {
coverages <- c(coverages, rnorm(1, mean=0, sd=variances[i]))
}

plot(coverages)


In real life, I know possible variances and means for each time point. I need to infer prevalence of one of the models inside the segment.

Briefly, the package mcp does Bayesian change point regression. As of v0.2, it takes Gaussian, Binomial, Bernoulli, and Poisson. Modeling your data as four intercept-only segments:

model = list(
y ~ 1,  # Intercept
~ 1,  # etc...
~ 1,
~ 1
)

library(mcp)
df = data.frame(x = seq_along(coverages), y = coverages)
fit = mcp(model, df, par_x = "x")


Let's plot it with a prediction interval, just for fun (green dashed lines). The blue curves are posterior densities for the change point locations. The gray lines are random draws from the posterior.

plot(fit, q_predict = T)


You can use plot_pars() to plot individual parameter estimates. Here are the summaries. where cp_* are the change point estimates:

summary(fit))

Iterations: 9000 from 3 chains.
Segments:
1: y ~ 1
2: y ~ 1 ~ 1
3: y ~ 1 ~ 1
4: y ~ 1 ~ 1

Population-level parameters:
name    mean  lower    upper Rhat n.eff
cp_1 101.280  99.38 103.0000    1  5627
cp_2 199.562 199.00 200.4314    1  5038
cp_3 299.365 296.85 301.7760    1  2340
int_1  -0.047  -0.11   0.0104    1  5614
int_2  -0.620  -0.68  -0.5592    1  5792
int_3   0.423   0.37   0.4838    1  6463
int_4  -0.018  -0.04   0.0036    1  5382
sigma_1   0.295   0.28   0.3082    1  5963


Read more on the mcp website. Disclaimer: I am the developer of mcp.

• That's a great answer! Just a bit too late, I already wrote a paper using different method biorxiv.org/content/10.1101/837971v1 - so I hope next generation of researchers will use your package, I am sure it will work better than mine simple algorithm... Jan 10 '20 at 17:43
• Great that you found (and published!) a solution. Consider accepting an answer (not necessarily mine) if it solves the problem you originally posed. It makes it easier for future users to find solutions when they come here. Jan 10 '20 at 23:08

The two good papers on this subject are below:

These do not apply a clustering algorithm but take the interval (since the last change point) into account as you have asked for. And they work with parametric distributions. The paper by Adams and Mackay (the first one) also has the algorithm implemented in MatLab and Python.

• Thank you very much! It seems exactly like what I wanted. They look really nice! I will read them carefully, probably, my search is finished =) Feb 23 '16 at 13:56
• =( unfortunately, the thing that I did not notice in the first paper: We further assume that for each partition ρ, the data within it are i.i.d. from some probability distribution. For my situation I know distribution parameters, so I can answer the question: "what are parameters for time point $t$ and model $H_i$", but the data within is not i.i.d. so distributions for time points $t$ and $t+1$ can be different even for the same model $H_i$. Feb 23 '16 at 16:18

This is more of a comment than an answer but it's too long to be a comment:

The "bible" for sequential analysis is probably 2014's book Sequential Analysis: Hypothesis Testing and Changepoint Detection by Alexander Tartakovsky. It is seemingly exhaustive in its coverage of the topic.

http://www.amazon.com/Sequential-Analysis-Hypothesis-Changepoint-Probability-ebook/dp/B00MMOIWTS/ref=sr_1_1?ie=UTF8&qid=1445511005&sr=8-1&keywords=sequential+analysis+tartakovsky

That said, in June 2014 Columbia sponsored The Fifth International Workshop in Sequential Methodologies which brought together the latest and greatest practitioners in the field. Tartakovsky was on the organizing committee.

See the "Detailed Program" link on the conference website for abstracts and papers. There's probably something there targeted specifically to your question.

Response lifted from this CV thread as written by me: sequential estimators for a proportion

• I would like to buy this book, but can not. I have a stipend, not a salary...70 euro is quite a lot for me =( I'll try to find smth in the Detailed Program. Feb 23 '16 at 16:27
• Can you get it through interlibrary loan? Feb 25 '16 at 13:09

Numerous packages are available in R for changepoint or breakpoint detection but the majority of them are non-Bayesian. Many such packages are touched in @Jonas Lindeløv's blog post: https://lindeloev.github.io/mcp/articles/packages.html. Glad that he also pointed to his site in the answer above.

In addition to mcp, bcp is also a popular Bayesian changepoint detection model and it has been actually used to analyze copy-number alteration sequence data--the same use case as your publication. For completeness, I put some quick results here with your sample data:

set.seed(1234)
variances = runif(1000, 0.01, 0.5)
coverages = c()
for (i in seq(1:100)) {
coverages = c(coverages, rnorm(1, mean=0, sd=variances[i]))
}
for (i in seq(101:200)) {
coverages <- c(coverages, rnorm(1, mean=-log(2), sd=variances[i] / 0.75))
}
for (i in seq(201:300)) {
coverages <- c(coverages, rnorm(1, mean=log(3/2), sd=variances[i] * 0.75))
}
for (i in seq(301:1000)) {
coverages <- c(coverages, rnorm(1, mean=0, sd=variances[i]))
}

library(bcp)
out = bcp(coverages)
plot(out)


Here is the bcp output:

Another package for Bayesian changepoint detection is Rbeast, which handles only time-series or sequence-like data and therefore is not as versatile as mcp or bcp. But it is useful for your use scenario. Rbeast is aimed also to decompose time series into periodic and trend components. Since your sequence contains no periodic/seasonal component, season='none' is set in the code below and only the trend component is fitted:

library(Rbeast)
out = beast(coverages, season='none')
plot(out)


Below is the Rbeast output:

For each segment,Rbeast fits either a linear (1st order polynomial) or constant (0-th order poly) model; the average order of the polynomials needed to adequately fit the trend is estimated over time and depicted as the Order_t curve in the figure: all close to 0, suggesting an overall flat curve. For your sample time series, since we know the segments are constant beforehand, we can enforce this strong prior by setting the min and max order of the segments (polynomial) to both 0: tcp.minmax=c(0,0) so that only constant lines are fitted.

out = beast(coverages, season='none', tcp.minmax=c(0,0) )
plot(out)


New contributor
zhaokg is a new contributor to this site. Take care in asking for clarification, commenting, and answering. Check out our Code of Conduct.