4

I am looking for a way (or at least an R package) to perform Bayesian changepoint analysis with Reversible-jump MCMC approach.

I will apply this for detecting changepoints in Typhoon time series.

This is my reference paper: https://journals.ametsoc.org/doi/pdf/10.1175/JCLI-D-13-00744.1

Figure from the above paper

I would also like to plot the posterior probability mass function for each change point.

Here is the sample data for this:

structure(list(V1 = c(7L, 6L, 4L, 4L, 4L, 2L, 5L, 4L, 4L, 4L, 
6L, 7L, 8L, 6L, 10L, 7L, 9L, 5L, 1L, 4L, 5L, 5L, 2L, 2L, 5L, 
1L, 2L, 4L, 0L, 3L, 6L, 3L, 6L, 1L, 5L, 3L, 4L, 0L, 2L, 4L)), class = 
"data.frame", row.names = c(NA, 
-40L))

I found this R-package, but it is not applied for changepoint analysis: https://cran.r-project.org/web/packages/rjmcmc/rjmcmc.pdf

Can anyone point me in the right package or at least help me on how to do this in R? I'll appreciate any help.

Lyndz
  • 347
  • 1
  • 13
  • 30
  • I might be misunderstanding your question, but it looks like all you need to create is an histogram using [hist](https://www.rdocumentation.org/packages/graphics/versions/3.6.2/topics/hist) as long as you already have the data from the mcmc. Otherwise, [this package](https://cran.r-project.org/web/packages/coda/index.html) might be helpful. – Trusky Jun 04 '20 at 15:35
  • @Trusky--i checked that one, but there is no reversible-jump MCMC. – Lyndz Jun 04 '20 at 15:38

1 Answers1

0

Just in case that extra thoughts are still needed for this old question, the majority of Bayesian changepoint or breakpoint detection packages in R are implemented via Gibbs sampling not Reversible-Jump MCMC sampling. One exception is a package Rbeast (https://cran.r-project.org/web/packages/Rbeast/index.html) I implemented that uses a hybrid Reversible-Jump MCMC sampler to estimate changepoint probabilities. As a caveat, Rbeast is formulated for Gaussian-like data not Poisson data (i.e., count data in your reference paper). In any case, you can still test it on count data, for example, by log-transforming. In SouthWood's book Ecological Methods (fifth edition), an example was given on Page 475 to transform count data via square root to make it more Gaussian-like before applying Rbeast.

Here are some sample code and quick results using the sample data you provided. Note that Rbeast does both time series decomposition (if a periodic component is present) and changepoint detection simultaneously, which makes it different from stl or changepoint functions that either just decompose time series or only detect breaks. Because your sample data has no seasonal/periodic component, in the example below, season='none' is specified in the beast function:

Y= c(7L, 6L, 4L, 4L, 4L, 2L, 5L, 4L, 4L, 4L, 6L, 7L, 8L, 6L, 10L, 7L, 9L, 5L, 1L, 
     4L, 5L, 5L, 2L, 2L, 5L, 1L , 2L, 4L, 0L, 3L, 6L, 3L, 6L, 1L, 5L,3L, 4L, 0L, 2L, 4L)
        
library(Rbeast)

# Rbeast is also a tool for decomposing a time series into seasonal and trend components, 
# but here your data is trend-only, so season='none' is used.
out=beast(Y,season='none')
plot(out)

Sample output

The vertical dashed lines pinpoint the most probable locations of changepoints. The green Pr(tcp) graph depicts the probability of changepoint occurring at individual point of time, which corresponds roughly to Figure d and e of your referenced paper but without differentiation of the orders of changepoints (e.g. first changepoint, second changepoint, ...).

The posterior probability mass function you look for is stored in the output out$trend$ncpPr. Here is a bar plot of it.

ncpPr=out$trend$ncpPr[,1]
ncp  = (1:length(ncpPr)) -1
barplot( ncpPr ~ ncp, data=data.frame(ncp, ncpPr), xlab='num of changepoints', ylab='posterior prob')

enter image description here

Another great package for this kind of analysis is bfast, which is not MCMC-based.

zhaokg
  • 282
  • 2
  • 14