Change point package in R using Reversible-Jump MCMC Bayesian approach

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)

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')

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