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.
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.056375 sec
#> Completed: 400/2000 - in 0.108558 sec
#> Completed: 600/2000 - in 0.160803 sec
#> Completed: 800/2000 - in 0.216577 sec
#> Completed: 1000/2000 - in 0.267866 sec
#> Completed: 1200/2000 - in 0.317435 sec
#> Completed: 1400/2000 - in 0.367638 sec
#> Completed: 1600/2000 - in 0.416679 sec
#> Completed: 1800/2000 - in 0.469027 sec
#> Completed: 2000/2000 - in 0.519497 secWith 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: 8.9
#> - Computational time: 0.52 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
#> 49 64 130 148 301 307 319 325The package also provides a method for plotting the change points.
We can assess convergence of the latent order posterior chain, for
example, by inspecting the traceplot of its log-likelihood with
coda::traceplot.
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.03948 sec
#> Completed: 400/2000 - in 0.078843 sec
#> Completed: 600/2000 - in 0.119545 sec
#> Completed: 800/2000 - in 0.160387 sec
#> Completed: 1000/2000 - in 0.202748 sec
#> Completed: 1200/2000 - in 0.245142 sec
#> Completed: 1400/2000 - in 0.287592 sec
#> Completed: 1600/2000 - in 0.33003 sec
#> Completed: 1800/2000 - in 0.372545 sec
#> Completed: 2000/2000 - in 0.415026 sec
table(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 2 3 4 5
#> 42 134 102 24 34plot(out, loss = "binder", plot_freq = TRUE)
#> Warning in salso::salso(mcmc_chain, loss = "binder", maxNClusters =
#> maxNClusters, : The number of clusters equals the default maximum possible
#> number of clusters.
Function detect_cp can also be used to detect change
points on survival functions. We consider the synthetic dataset
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 1.90273 sec
#> Completed: 400/2000 - in 3.70433 sec
#> Completed: 600/2000 - in 5.50799 sec
#> Completed: 800/2000 - in 7.29868 sec
#> Completed: 1000/2000 - in 9.09697 sec
#> Completed: 1200/2000 - in 10.8842 sec
#> Completed: 1400/2000 - in 12.7115 sec
#> Completed: 1600/2000 - in 14.5795 sec
#> Completed: 1800/2000 - in 16.3786 sec
#> Completed: 2000/2000 - in 18.163 sec
print(out)
#> DetectCpObj object
#> Type: change points detection on an epidemic diffusionAlso here, with function plot we can plot the survival
function and the position of the change points.
