---
title: 'BayesChange Package: Bayesian Change Point Analysis'
url: https://www.emergentmind.com/topics/bayeschange-package
type: topic
---

# BayesChange Package: Bayesian Change Point Analysis

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
  $$
  p(\rho\mid y)\;\propto\;\Pr(\rho)\prod_{j=1}^m p\left(y_{A_j}\mid \phi, \text{(priors)}\right)
  $$

### 1.2 Clustering by Change Point Structure

With \(n\) curves \(\{y^{(i)}\}\), 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 \(2^{T-1}\)):
- Each cluster \(B_k\) shares enforced change point configuration \(\rho^{*(k)}\).
- Cluster allocation prior: \(\pi\sim\Dirichlet(\alpha,\dots,\alpha)\).

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 \(q\), selects a block and splits at a random interior location.
2. **Merge**: Otherwise, merges two consecutive blocks.
3. **Shuffle**: When \(m>1\), shuffles assignments across two neighboring blocks, keeping block sizes fixed.

Acceptance probability uses the Metropolis–Hastings ratio:
$$
\alpha = \min\left\{1, \frac{\Pr(\rho^\text{new})\,p(y\mid\rho^\text{new})}{\Pr(\rho^\text{old})\,p(y\mid\rho^\text{old})} \cdot \frac{q_{\rm back}}{q_{\rm fwd}} \right\}
$$
Hyperparameters $(\sigma, \delta)$ and $\phi$ are updated via Metropolis–Hastings or Gibbs steps.

### 2.2 Clustering Sampler

The clustering algorithm jointly samples the allocation partition (\(\lambda\)) and the set of cluster atoms (\(\{\rho^{*(k)}\}\)), 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 $\rho^*$ from an instrumental mixture of single-series posteriors.
3. Accepts or rejects via the full joint posterior.
4. Periodically resamples each $\rho^{*(k)}$ 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 $O$(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-$T$ integer label vectors with corresponding break-point indexes; only local data is recomputed at each step, minimizing unnecessary recalculation.

This architecture yields linear-in-$T$ 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 (\(T=200\)) with known change points is as follows:

```r
library(BayesChange)
set.seed(123)
T<-200; phi<-0.1
cp <- c(1,51,151,201)
data <- numeric(T)
for (s in seq_len(length(cp)-1)) {
  mu <- c(0,1.5,0)[s]
  sd_ <- c(0.13,0.15,0.12)[s]
  for (t in cp[s]:(cp[s+1]-1)) {
    if (t==cp[s]) data[t]<-rnorm(1,mu,sd_)
    else data[t] <- phi*data[t-1] + (1-phi)*mu + rnorm(1,0,sd_*sqrt(1-phi^2))
  }
}
params_uni <- list(a=1,b=1,c=1, prior_var_phi=0.05, prior_delta_c=1, prior_delta_d=1)
out <- detect_cp(data, n_iterations=8000, n_burnin=3000, q=0.25, params=params_uni, kernel="ts", print_progress=FALSE)
cp_est_labels <- posterior_estimate(out, loss="binder")
unique_blocks <- split(seq_along(cp_est_labels), cp_est_labels)
est_cps <- cumsum(sapply(unique_blocks,length))[-length(unique_blocks)] + 1
print(est_cps) # Should be ≈ c(51,151)
plot(out, loss="binder", plot_freq=TRUE)
```
Convergence diagnostics include traceplots (e.g., `out$phi_MCMC`), summarizing acceptance rates, and inspecting the frequency of time steps being flagged as change points.

## 6. Performance and Comparative Assessment

Benchmarking against other packages:

| Package         | Dimension       | Features Offered                         | Relative Speed    |
|-----------------|----------------|------------------------------------------|-------------------|
| **bcp**         | Univariate      | Change point detection                   | 5× slower         |
| **ppmSuite**    | Univariate/multi| Some PPM-based detection                 | Lacks clustering, less efficient |
| **HDcpDetect**  | Univariate      | Frequentist, no PPM prior                | —                 |

On $T=500$ (real) and $1000$ (sim) iterations, BayesChange is 5× faster than bcp and 8× faster for $d=3$ multivariate data. Clustering across $n=20$ series of length $T=200$ typically takes $\approx$3 seconds for $5000$ 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**: $\sigma, \delta$ modulate the expected number of blocks. $\sigma \approx 0, \delta > 0$ yields a Dirichlet-process-like regime favoring few blocks; larger $\sigma$ intensifies the prior weight on many, smaller blocks.
- **Hyperparameter Choices**: $(a, b, c) = 1$ yields diffuse Normal–Gamma priors when uncertain; stronger priors should reflect external information.
- **Clustering Dirichlet Weight**: Low $\alpha_{\rm SM}$ (<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.

Source: https://www.emergentmind.com/topics/bayeschange-package