BayesChange Tutorial

Here we provide a brief tutorial of the BayesChange package. The BayesChange package contains two main functions: one that performs change points detection on time series and epidemic diffusions and one that perform clustering of time series and epidemic diffusions with common change points. Here we briefly show how to implement these.

library(BayesChange)

Detecting change points

The function detect_cp provide a method for detecting change points, it is based on the work Martínez and Mena (2014) and on Corradin et al. (2024)

Depending on the structure of the data, detect_cp might perform change points detection on univariate time series or multivariate time series. We import dataset eu_inflation that contains the standardized monthly inflation rates from the Harmonized Index of Consumer Prices (HICP) for the 12 COICOP expenditure categories across European Union countries. The data span the period from February~1997 to December~2024 resulting in a matrix of \(12\) rows and \(355\) columns.

data("eu_inflation")

Now we can run the function detect_cp, as arguments of the function we need to specify the number of iterations, the number of burn-in steps and a list with the the autoregressive coefficient phi for the likelihood of the data, the parameters a, b, c for the priors and the probability q of performing a split at each step. Since we deal with time series we need also to specify kernel = "ts".

out <- detect_cp(data = eu_inflation[1,],
                 n_iterations = 2000, n_burnin = 500, q = 0.5,
                 params = list(prior_var_phi = 0.1, prior_delta_c = 1, prior_delta_d = 1), kernel = "ts")
#> Completed:   200/2000 - in 0.051367 sec
#> Completed:   400/2000 - in 0.104224 sec
#> Completed:   600/2000 - in 0.151622 sec
#> Completed:   800/2000 - in 0.197464 sec
#> Completed:   1000/2000 - in 0.248899 sec
#> Completed:   1200/2000 - in 0.299252 sec
#> Completed:   1400/2000 - in 0.352791 sec
#> Completed:   1600/2000 - in 0.413137 sec
#> Completed:   1800/2000 - in 0.470574 sec
#> Completed:   2000/2000 - in 0.526228 sec

With the methods print and summary we can get information about the algorithm.

print(out)
#> DetectCpObj object
#> Type: change points detection on univariate time series
summary(out)
#> DetectCpObj object
#> Change point detection summary:
#>  - Data: univariate time series
#>  - Burn-in iterations: 500 
#>  - MCMC iterations: 1500 
#>  - Average number of detected change points: 7.63 
#>  - Computational time: 0.53 seconds
#>  
#> Use plot() for a detailed visualization or posterior_estimate() to analyze the detected change points.

In order to get a point estimate of the change points we can use the method posterior_estimate that uses the method salso by David B. Dahl and Müller (2022) to get the final latent order and then detect the change points.

cp_est <- posterior_estimate(out, loss = "binder")
cumsum(table(cp_est))[-length(table(cp_est))] + 1
#>   1   2   3   4   5   6   7   8 
#>  45 203 241 301 305 322 323 326

The package also provides a method for plotting the change points.

plot(out, loss = "binder")

We can assess convergence of the latent order posterior chain, for example, by inspecting the traceplot of its log-likelihood with coda::traceplot.

coda::traceplot(out$lkl_MCMC, ylab = "Log-Likelihood")

If we instead consider a matrix of data, detect_cp automatically performs a multivariate change points detection method. We define the parameters.

params_multi <- list(m_0 = rep(0,3),
                     k_0 = 1,
                     nu_0 = 10,
                     S_0 = diag(0.1,3,3),
                     prior_var_phi = 0.1,
                     prior_delta_c = 1,
                     prior_delta_d = 1)

Arguments k_0, nu_0, phi_0, m_0, prior_delta_c, prior_delta_d and prior_var_phi correspond to the parameters of the prior distributions for the multivariate likelihood.

out <- detect_cp(data = eu_inflation[1:3,], n_iterations = 2000,
          n_burnin = 500, q = 0.5, params = params_multi, kernel = "ts")
#> Completed:   200/2000 - in 0.091843 sec
#> Completed:   400/2000 - in 0.185125 sec
#> Completed:   600/2000 - in 0.278698 sec
#> Completed:   800/2000 - in 0.372207 sec
#> Completed:   1000/2000 - in 0.46859 sec
#> Completed:   1200/2000 - in 0.565817 sec
#> Completed:   1400/2000 - in 0.663724 sec
#> Completed:   1600/2000 - in 0.761897 sec
#> Completed:   1800/2000 - in 0.859877 sec
#> Completed:   2000/2000 - in 0.961587 sec

table(posterior_estimate(out, loss = "binder"))
#> 
#>   1   2   3   4   5   6 
#>  40 124  36  78  24  34
plot(out, loss = "binder", plot_freq = TRUE)

Function detect_cp can also be used to detect change points on survival functions. We consider the synthetic dataset epi_synthetic

data("epi_synthetic")

To run detect_cp on epidemiological data we need to set kernel = "epi". Moreover, besides the usual parameters, we need to set the number of Monte Carlo replications M for the approximation of the integrated likelihood and the recovery rate xi. a0 and b0 are optional and correspond to the parameters of the gamma distribution for the integration of the likelihood.

params_epi <- list(M = 250, xi = 1/8, a0 = 4, b0 = 10, I0_var = 0.1)

out <- detect_cp(data = epi_synthetic, n_iterations = 2000, n_burnin = 500,
                 q = 0.25, params = params_epi, kernel = "epi")
#> Completed:   200/2000 - in 3.47599 sec
#> Completed:   400/2000 - in 7.11926 sec
#> Completed:   600/2000 - in 10.838 sec
#> Completed:   800/2000 - in 14.322 sec
#> Completed:   1000/2000 - in 17.9212 sec
#> Completed:   1200/2000 - in 21.4089 sec
#> Completed:   1400/2000 - in 24.8651 sec
#> Completed:   1600/2000 - in 28.3214 sec
#> Completed:   1800/2000 - in 31.7704 sec
#> Completed:   2000/2000 - in 35.312 sec

print(out)
#> DetectCpObj object
#> Type: change points detection on an epidemic diffusion

Also here, with function plot we can plot the survival function and the position of the change points.

plot(out)

Clustering time dependent data with common change points

BayesChange contains another function, clust_cp, that cluster respectively univariate and multivariate time series and survival functions with common change points. Details about this methods can be found in Corradin et al. (2026)

In clust_cp the argument kernel must be specified, if data are time series then kernel = "ts" must be set. Then the algorithm automatically detects if data are univariate or multivariate.

We consider for this example dataset stock_uni that contains the daily mean stock prices for the 50 largest companies (by market capitalization) in the Standard&Poor’s 500 Index from January 1, 2020 to January 1, 2022.

data("stock_uni")

Arguments that need to be specified in clust_cp are the number of iterations n_iterations, the number of elements in the normalisation constant B, the split-and-merge step L performed when a new partition is proposed and a list with the parameters of the algorithm, the likelihood and the priors..

params_uni <- list(a = 1,
                   b = 1,
                   c = 1,
                   phi = 0.1)

out <- clust_cp(data = stock_uni[1:5,], n_iterations = 2000, n_burnin = 500,
                L = 1, q = 0.5, B = 1000, params = params_uni, kernel = "ts")
#> Normalization constant - completed:  100/1000 - in 0.018069 sec
#> Normalization constant - completed:  200/1000 - in 0.036561 sec
#> Normalization constant - completed:  300/1000 - in 0.056982 sec
#> Normalization constant - completed:  400/1000 - in 0.074768 sec
#> Normalization constant - completed:  500/1000 - in 0.093225 sec
#> Normalization constant - completed:  600/1000 - in 0.112537 sec
#> Normalization constant - completed:  700/1000 - in 0.131584 sec
#> Normalization constant - completed:  800/1000 - in 0.151159 sec
#> Normalization constant - completed:  900/1000 - in 0.171933 sec
#> Normalization constant - completed:  1000/1000 - in 0.191917 sec
#> 
#> ------ MAIN LOOP ------
#> 
#> Completed:   200/2000 - in 0.448281 sec
#> Completed:   400/2000 - in 0.803392 sec
#> Completed:   600/2000 - in 1.16513 sec
#> Completed:   800/2000 - in 1.50551 sec
#> Completed:   1000/2000 - in 1.84308 sec
#> Completed:   1200/2000 - in 2.17524 sec
#> Completed:   1400/2000 - in 2.53146 sec
#> Completed:   1600/2000 - in 2.87629 sec
#> Completed:   1800/2000 - in 3.22044 sec
#> Completed:   2000/2000 - in 3.55272 sec

posterior_estimate(out, loss = "binder")
#> Warning in salso::salso(mcmc_chain, loss = "binder", maxNClusters =
#> maxNClusters, : The number of clusters equals the default maximum possible
#> number of clusters.
#> [1] 1 2 3 1 1

Method plot for clustering univariate time series represents the data colored according to the assigned cluster.

plot(out, loss = "binder")
#> Warning in salso::salso(mcmc_chain, loss = "binder", maxNClusters =
#> maxNClusters, : The number of clusters equals the default maximum possible
#> number of clusters.

Method plot_psm shows the posterior similarity matrix of the clustering. Selecting reorder = TRUE we can choose to order the matrix depending on the clustering obtained.

plot_psm(out, reorder = TRUE)

If time series are multivariate, data must be an array, where each element is a multivariate time series represented by a matrix. Each row of the matrix is a component of the time series. Here we use dataset stock_multi that contains for each company the daily opening and closing stock prices.

data("stock_multi")
params_multi <- list(m_0 = rep(0,2),
                     k_0 = 1,
                     nu_0 = 10,
                     S_0 = diag(1,2,2),
                     phi = 0.1)

out <- clust_cp(data = stock_multi[,,1:5], n_iterations = 2500, n_burnin = 500,
                L = 1, B = 1000, params = params_multi, kernel = "ts")
#> Normalization constant - completed:  100/1000 - in 0.035732 sec
#> Normalization constant - completed:  200/1000 - in 0.071291 sec
#> Normalization constant - completed:  300/1000 - in 0.106814 sec
#> Normalization constant - completed:  400/1000 - in 0.142554 sec
#> Normalization constant - completed:  500/1000 - in 0.178019 sec
#> Normalization constant - completed:  600/1000 - in 0.213486 sec
#> Normalization constant - completed:  700/1000 - in 0.249102 sec
#> Normalization constant - completed:  800/1000 - in 0.284532 sec
#> Normalization constant - completed:  900/1000 - in 0.319914 sec
#> Normalization constant - completed:  1000/1000 - in 0.355925 sec
#> 
#> ------ MAIN LOOP ------
#> 
#> Completed:   250/2500 - in 0.812958 sec
#> Completed:   500/2500 - in 1.60206 sec
#> Completed:   750/2500 - in 2.42219 sec
#> Completed:   1000/2500 - in 3.2348 sec
#> Completed:   1250/2500 - in 4.03843 sec
#> Completed:   1500/2500 - in 4.83728 sec
#> Completed:   1750/2500 - in 5.6336 sec
#> Completed:   2000/2500 - in 6.43706 sec
#> Completed:   2250/2500 - in 7.23635 sec
#> Completed:   2500/2500 - in 8.04338 sec

posterior_estimate(out, loss = "binder")
#> [1] 1 2 3 4 1
plot(out, loss = "binder")

Finally, if we set kernel = "epi", clust_cp cluster survival functions with common change points. Also here details can be found in Corradin et al. (2026)

Data are a matrix where each row is the number of infected at each time. Inside this package is included the dataset epi_synthetic_multi with multivariate synthetic epidemic diffusions.

data("epi_synthetic_multi")

params_epi <- list(M = 100, xi = 1/8,
                   alpha_SM = 1,
                   a0 = 4,
                   b0 = 10,
                   I0_var = 0.1,
                   avg_blk = 2)

out <- clust_cp(epi_synthetic_multi[,10:150], n_iterations = 2000, n_burnin = 500,
                L = 1, B = 1000, params = params_epi, kernel = "epi")

posterior_estimate(out, loss = "binder")
plot(out, loss = "binder")
Corradin, Riccardo, Luca Danese, Wasiur R. KhudaBukhsh, and Andrea Ongaro. 2024. Model-Based Clustering of Time-Dependent Observations with Common Structural Changes. https://arxiv.org/abs/2410.09552.
Corradin, Riccardo, Luca Danese, Wasiur R. KhudaBukhsh, and Andrea Ongaro. 2026. “Model-Based Clustering of Time-Dependent Observations with Common Structural Changes.” Statistics and Computing 36 (1): 7. https://doi.org/10.1007/s11222-025-10756-x.
David B. Dahl, Devin J. Johnson, and Peter Müller. 2022. “Search Algorithms and Loss Functions for Bayesian Clustering.” Journal of Computational and Graphical Statistics 31 (4): 1189–201. https://doi.org/10.1080/10618600.2022.2069779.
Martínez, Asael Fabian, and Ramsés H. Mena. 2014. “On a Nonparametric Change Point Detection Model in Markovian Regimes.” Bayesian Analysis 9 (4): 823–58. https://doi.org/10.1214/14-BA878.