BayesChange Package: Bayesian Change Point Analysis

• BayesChange is a comprehensive framework implementing fully Bayesian change point detection using product-partition models and exact split-merge MCMC.
• It employs a Pitman–Yor process prior on ordered partitions to effectively cluster time-dependent curves based solely on change point locations.
• The package leverages C++-backed routines and optimized algorithms to efficiently analyze univariate, multivariate, and epidemic SIR count data.

BayesChange is a computationally efficient R package—backed by C++—providing fully Bayesian, product-partition-based methods for change point detection and clustering in univariate and multivariate time series, as well as in daily epidemic SIR count data. Uniquely, BayesChange can cluster multiple time-dependent curves or survival functions by partitioning them solely according to the locations of their change points. The package centers on exact split-merge Markov chain Monte Carlo (MCMC) with analytic marginal likelihoods, a Pitman–Yor process prior on ordered partitions, and an R interface with S3 methods for post-processing, visualization, and point estimation of change point and clustering structures.

1. Modeling Framework

1.1 Change Point Model for Single Series

Given a sequence $y = (y_1,\dots,y_T)$, univariate ($d=1$) or multivariate ($d>1$), the data are assumed to admit an unobserved ordered partition $\rho = \{A_1, ..., A_m\}$ of the time index into contiguous blocks, with each block modeled by a common parameter $\theta_j^*$.

• Partition Prior: The prior on partitions is based on the exchangeable partition probability function (eppf) of a Pitman–Yor process $\mathrm{PY}(\sigma,\delta)$, restricted to ordered partitions:

$$\Pr(\rho) = \frac{T!}{m!\,\prod_{j=1}^m |A_j|!}\;\frac{\prod_{j=1}^{m-1}(\delta + j\,\sigma)}{(\delta+1)_{T-1}\;\prod_{j=1}^m (1-\sigma)_{|A_j|-1}}$$

with discount parameter $\sigma\in(0,1)$ and strength $\delta>-\sigma$.

• Block-wise Likelihood: Within each block, data are generated under an AR(1)/Ornstein–Uhlenbeck model:

$$y_t\mid y_{t-1},\theta_j^* \sim \mathcal N\left(\mu_j+\phi(y_{t-1}-\mu_j),\,(1-\phi^2)\Sigma_j\right)$$

with analytically marginalizable conjugate Normal–Gamma (univariate) or Normal–Inverse-Wishart (multivariate) priors.

• Posterior: The joint posterior is

$d=1$0

1.2 Clustering by Change Point Structure

With $d=1$1 curves $d=1$2, clustering is performed by grouping curves that share identical change points. The associated hierarchical model places a Dirichlet mixture prior over all possible ordered partitions (of size $d=1$3):

• Each cluster $d=1$4 shares enforced change point configuration $d=1$5.
• Cluster allocation prior: $d=1$6.

This structure supports clustering solely by shared change-point location patterns—a capability not present in prior R packages.

2. Posterior Simulation Methodology

2.1 Split-Merge MCMC for Single Series

The sampler performs MCMC over ordered partitions using three steps at each iteration:

1. Split: With probability $d=1$7, selects a block and splits at a random interior location.
2. Merge: Otherwise, merges two consecutive blocks.
3. Shuffle: When $d=1$8, shuffles assignments across two neighboring blocks, keeping block sizes fixed.

Acceptance probability uses the Metropolis–Hastings ratio:

$d=1$9

Hyperparameters $d>1$0 and $d>1$1 are updated via Metropolis–Hastings or Gibbs steps.

2.2 Clustering Sampler

The clustering algorithm jointly samples the allocation partition ($d>1$2) and the set of cluster atoms ($d>1$3), as follows:

1. Selects two curves at random; if in the same cluster, split their cluster; otherwise, merge.
2. For each new cluster, samples a new $d>1$4 from an instrumental mixture of single-series posteriors.
3. Accepts or rejects via the full joint posterior.
4. Periodically resamples each $d>1$5 conditional on its cluster's curves, using internal split-merge MCMC.

Both algorithms are detailed in the package appendix.

3. Implementation Strategies and Computational Aspects

• Backend: All core routines are implemented in C++ via Rcpp, RcppArmadillo, and RcppGSL for high-efficiency linear algebra and random variate generation.
• Key Functions: Univariate, multivariate, and epidemic SIR variants—detect_cp_uni, detect_cp_multi, and detect_cp_epi—execute MCMC in explicit for-loops with blockwise MH steps. Helper functions (e.g., AlphaSplit_UniTS) permit $d>1$6(block-size) per update.
• Clustering functions: (clust_cp_uni, clust_cp_multi, clust_cp_epi) include a second-level split-merge MCMC and precompute Dirichlet mixture normalization constants.
• Data Structures: Partitions are maintained as length-$d>1$7 integer label vectors with corresponding break-point indexes; only local data is recomputed at each step, minimizing unnecessary recalculation.

This architecture yields linear-in-$d>1$8 scaling and enables practical MCMC sampling for hundreds or thousands of time points and series.

4. R-Level User Interface and Workflow

4.1 Wrappers

• detect_cp() for single-series detection:
  • Inputs: data, n_iterations, params (hyperparameters), kernel type, etc.
• clust_cp() for curve clustering:
  • Inputs: data (matrix/array), n_iterations, alpha_SM (Dirichlet weight), etc.

4.2 Hyperparameters

• Univariate detection: params = list(a, b, c, prior_var_phi, prior_delta_c, prior_delta_d)
• Multivariate: params = list(m_0, k_0, nu_0, S_0, prior_var_phi, prior_delta_c, prior_delta_d)
• Epidemic SIR: params = list(M, xi, a0, b0, I0_var)
• Clustering: additionally set alpha_SM, normalization control parameters B, L.

4.3 Output and S3 Methods

Method	Description	Typical Class
detect_cp()	Returns MCMC block-labels and parameter chains	"DetectCpObj"
clust_cp()	Returns cluster labels and partition trajectories	"ClustCpObj"
posterior_estimate(obj, ...)	SALSO search for “best” partition/cluster
plot()	ggplot2 overlay of estimated change points/clusters

Further S3 methods (print, summary) document runtime diagnostics and sampling details. Plotting functions provide raw CP frequency as well as loss-minimizing point estimates.

5. Empirical Example in R

A typical analysis for a univariate time series ($d>1$9) with known change points is as follows:

$\theta_j^*$2 Convergence diagnostics include traceplots (e.g., out$$phi_MCMC</code>), summarizing acceptance rates, and inspecting the frequency of time steps being flagged as change points.</p> <h2 class='paper-heading' id='performance-and-comparative-assessment'>6. Performance and Comparative Assessment</h2> <p>Benchmarking against other packages:</p> <div class='overflow-x-auto max-w-full my-4'><table class='table border-collapse w-full' style='table-layout: fixed'><thead><tr> <th>Package</th> <th>Dimension</th> <th>Features Offered</th> <th>Relative Speed</th> </tr> </thead><tbody><tr> <td><strong>bcp</strong></td> <td>Univariate</td> <td>Change point detection</td> <td>5× slower</td> </tr> <tr> <td><strong>ppmSuite</strong></td> <td>Univariate/multi</td> <td>Some PPM-based detection</td> <td>Lacks clustering, less efficient</td> </tr> <tr> <td><strong>HDcpDetect</strong></td> <td>Univariate</td> <td>Frequentist, no <a href="https://www.emergentmind.com/topics/probabilistic-projective-measurement-ppm" title="" rel="nofollow" data-turbo="false" class="assistant-link" x-data x-tooltip.raw="">PPM</a> prior</td> <td>—</td> </tr> </tbody></table></div> <p>On$$\rho = \{A_1, ..., A_m\}$0 (real) and $\rho = \{A_1, ..., A_m\}$1 (sim) iterations, BayesChange is 5× faster than bcp and 8× faster for $\rho = \{A_1, ..., A_m\}$2 multivariate data. Clustering across $\rho = \{A_1, ..., A_m\}$3 series of length $\rho = \{A_1, ..., A_m\}$4 typically takes $\rho = \{A_1, ..., A_m\}$53 seconds for $\rho = \{A_1, ..., A_m\}$6 MCMC iterations. Editor's term: These numbers represent order-of-magnitude improvements in wall-clock run time over prior R packages for Bayesian change point tasks of comparable complexity.

BayesChange is indicated when: (i) full Bayesian posterior uncertainty (over number and position of change points) is required, (ii) the data are multivariate or have simultaneous change points, (iii) clustering by change-point location pattern is needed, or (iv) there is a preference for C++-backed sampling at scale.

7. Tuning Guidelines, Interpretation, and Extensibility

• Prior Controls: $\rho = \{A_1, ..., A_m\}$7 modulate the expected number of blocks. $\rho = \{A_1, ..., A_m\}$8 yields a Dirichlet-process-like regime favoring few blocks; larger $\rho = \{A_1, ..., A_m\}$9 intensifies the prior weight on many, smaller blocks.
• Hyperparameter Choices: $\theta_j^*$0 yields diffuse Normal–Gamma priors when uncertain; stronger priors should reflect external information.
• Clustering Dirichlet Weight: Low $\theta_j^*$1 (<1) leans toward fewer cluster-specific CP patterns; high values allow for more differentiation.
• Extensions: To introduce novel block likelihoods (e.g., for Poisson counts with seasonality), one writes the marginal block likelihood in C++, exports it, and creates an R wrapper with a new kernel.

These features make BayesChange a comprehensive and computationally efficient solution for unsupervised change point detection and clustering in time series and count/survival models, with robust uncertainty quantification and flexible architecture for methodological extension.