Notice of Update

I have made some changes to the Bayesian Blocks code linked from here, on 24th November 2021. Also I note the coming and going of a “BayesianBlocks” package on CRAN which contained an optinterval function also based upon PELT.  It has been withdrawn from CRAN. As in the case of most such withdrawals, it probably happened because no one is willing to maintain it.

The Post

The Bayesian blocks algorithm of Scargle, Jackson, Norris, and Chiang has an enthusiastic user community in astrostatistics, in data mining, and among some in machine learning. It is a dynamic programming algorithm (see VanderPlas referenced below) and, so, exhibits optimality when used without performance-related shortcuts. There is an implementation by Scargle and others in MATLAB, and a popular version exists in the astropy and astroML modules for Python. Scargle introduced the idea about 2001, and has repeatedly improved it. I searched for a comparable implementation for R, the programming language I use most heavily and could not find one.

In time, however, I discovered the work of Killick, Fearnhead, and Eckley, led to it by an oceanographer colleague. They build upon the work of Jackson, Scargle, Barnes, Arabhi, Alt, Gioumousis, Gwin, Sangtrakulcharoen, Tan, and Tsai, but find a pruning step which lets the algorithm achieve the same result but in linear time. Their algorithm, known as PELT, has been incorporated into two R packages on CRAN, changepoint and changepoint.np, and has found use in oceanography and other areas. Still, a one-for-one substitute for the Scargle-Jackson-Norris-Chiang Bayesian blocks was missing, and it is the purpose of this post to remedy that.

• R. Killick, P. Fearnhead, And I. A. Eckley, “Optimal detection of changepoints with a linear computational cost”, Journal of the American Statistical Association, December 2012, 107(500), DOI: 10.1080/01621459.2012.737745
• R. Killick, I. A. Eckley, K. Ewans, P. Johnson, “Detection of changes in the characteristics of oceanographic time-series using changepoint analysis”, Ocean Engineering, 2010, 37(13), 10.1016/j.oceaneng.2010.04.009
• R. Killick, I. A. Eckley, “changepoint: An R package for changepoint analysis”, Journal of Statistical Software, June 2014, 58(3), 10.18637/jss.v058.i03
• B. Jackson, J. D. Scargle, D. Barnes, S. Arabhi, A. Alt, P. Gioumousis, E. Gwin, P. Sangtrakulcharoen, L. Tan, T. T. Tsai, “An algorithm for optimal partitioning of data on an interval”, IEEE Signal Processing Letters, February 2005, 12(2), 10.1109/LSP.2001.838216
• J. D. Scargle, J. P. Norris, B. Jackson, J. Chiang, “Studies in astronomical time series analysis. VI. Bayesian Block representations”, The Astrophysical Journal, 2013 February 20, 764(167), 10.1088/0004-637X/764/2/167
• J. VanderPlas, “Dynamic programming in Python: Bayesian Blocks”, Pythonic Perambulations, 12 September 2012, a blog entry
• bayesian_blocks, a component of astropy.stats, the astropy module for Python
• 4.1 Bayesian Blocks: Histograms the right way, a component of 4. Unsupervised Learning: Density Estimation, the astroML module for Python
• J. D. Scargle, “Bayesian Blocks in two or more dimensions: Image segmentation and cluster analysis”, arXiv:math/0111128v1

Jake VanderPlas is one of the astroML authors and, specifically that of the Bayesian blocks implementation therein. He offers a description of the Scargle algorithm, explaining it as an instance of a dynamic programming code pattern. He also offers a clever mixture distribution consisting of five Cauchy distributions for illustrating his code, a mixture distribution I have adopted for illustrating this, my implementation of Bayesian blocks in R. (R code is at link.)

Comparing its performance to that of a standard histogramming facility, truehist in the MASS package, we see the following:

(Click on figure to see a larger image, then use browser Back Button to return to blog.)

The result from VanderPlas’ implementation can be seen at his blog, just above the Conclusion at its bottom.

But VanderPlas’ mixture also can serve as a series for another test. Bayesian blocks has been proposed as a signals representation as well, one which can be used for compression as well as specifying features. The implementation here, based upon the Killick-Fearnhead-Eckley PELT, can do that as well, as can be seen by the figure below:

(Click on figure to see a larger image, then use browser Back Button to return to blog.)

Important Note
Alas, the PELT implementation used for this is based upon the PELT method of cpt.meanvar in the changepoint package of R, available from CRAN. That method apparently works for nothing but the Gamma model. In particular, it breaks badly for the Normal model. Thanks to Ben Groebe for pointing this out. He tried it on these data:

Note how Gamma-like the rendering is.

The changepoint package does handle Gaussian data via cpt.meanvar, but the PELT method is not used in that case (and, so, no Bayesian blocks) but, rather, the BinSeg method. (See the documentation for cpt.meanvar.)

The code is available from my Google space as a gzip‘d R file. Please feel free to use for anything you like and, not only would I be happy to answer questions about it, I’d like to know what kinds of things you are using it for, or receiving suggestions for improvements. Please comment below!

I am also intending to generalize the code so it can be used for multidimensional Bayesian blocks, as Scargle has suggested (in J. D. Scargle, “Bayesian Blocks in two or more dimensions: Image segmentation and cluster analysis”, arXiv:math/0111128v1), but, of course, based upon the Killick-Fearnhead-Eckley PELT. My particular interest is cluster-finding.

Considered Remark

While I have used this Bayesian blocks implementation via PELT for a couple of useful things, I have re-thought the entire changepoint detection philosophy which underlies its operation. In short, it is attempting to do a changepoint detection in a non-parametric manner, without bringing explicit knowledge into the definition of a changepoint. I s’ppose that this is possible in some general sense, if the data are like the series segmentation above, which was, after all, the kind of application which Scargle and, later, Killick, Fearnhead, and Eckley of PELT had in mind.

My preferred way of doing change detection now is via a longitudinal regression model tailored to the problem or series.