library(changepoints)

This is a simple guide for offline multiple change points detection on univariate mean.

There are 3 methods implemented for this purpose:

DP.univar: performs dynamic programming for univariate mean change points detection.
- CV.search.DP.univar: perform grid search to select the tuning parameter (associated with the l_0 penalty) through Cross-Validation.
BS.univar: performs standard binary segmentation for univariate mean change points detection.
- BS.univar.CPD: Perform univariate mean change points detection based on binary segmentation with the threshold parameter automatically selected based on the sSIC score.
WBS.univar: performs wild binary segmentation for univariate mean change points detection.
- WBS.univar.CPD: Perform univariate mean change points detection based on wild binary segmentation with the threshold parameter automatically selected based on the sSIC score.

In addition, function local.refine.univar performs local refinement for an initial change point estimation.

Simulate data

delta = 5 # 2*delta represents the minimum gap between boundaries
sigma2 = 1 # error variance

set.seed(1234)
y = c(rep(0, 50), rep(2, 50), rep(0, 50), rep(-2, 50)) + rnorm(200, mean = 0, sd = sqrt(sigma2)) # univariate time series
cpt_true = c(50, 100, 150)
n = length(y) # sample size

Perform dynamic programming

gamma_set = c(0.01, 0.5, 1, 5, 10, 50) # a set of tuning parameters for DP
DP_result = CV.search.DP.univar(y, gamma_set, delta) # grid search through cross-validation
min_idx = which.min(DP_result$test_error) # select gamma achieves the minimum validation error
cpt_DP_hat = unlist(DP_result$cpt_hat[[min_idx]]) # estimated change points by DP
cpt_DP_hat
#> [1]  53  99 149
Hausdorff.dist(cpt_DP_hat, cpt_true)
#> [1] 3
cpt_DPlr_hat = local.refine.univar(cpt_DP_hat, y) # perform local refinement
cpt_DPlr_hat
#> [1]  51 100 150
Hausdorff.dist(cpt_DPlr_hat, cpt_true)
#> [1] 1

Perform standard binary segmentation with user specified threshold parameter

tau_BS = 4 # user specified threshold parameter for BS
BS_object = BS.univar(y, 1, n, delta)
BS_result = thresholdBS(BS_object, tau_BS)
BS_result
#> $BS_tree_trimmed
#>                                       levelName
#> 1 Binary Segmentation Tree Trimmed with tau = 4
#> 2  °--N150                                     
#> 3      °--N51                                  
#> 4          °--N100                             
#> 
#> $cpt_hat
#>   location     value level
#> 1      150 15.408017     1
#> 2       51  8.753484     2
#> 3      100 10.720866     3
cpt_BS_hat = sort(BS_result$cpt_hat[,1]) # estimated change points by BS
cpt_BS_hat
#> [1]  51 100 150
Hausdorff.dist(cpt_BS_hat, cpt_true)
#> [1] 1
cpt_BSlr_hat = local.refine.univar(cpt_BS_hat, y) # perform local refinement
cpt_BSlr_hat
#> [1]  51 100 150
Hausdorff.dist(cpt_BSlr_hat, cpt_true)
#> [1] 1

Perform standard binary segmentation with automatically tuned threshold parameter

BS_CPD_result = tuneBSunivar(BS_object, y)
BS_CPD_result$cpt
#> [1]  51 100 150
Hausdorff.dist(BS_CPD_result$cpt, cpt_true)
#> [1] 1
BS_CPD_LR = local.refine.univar(BS_CPD_result$cpt, y)
BS_CPD_LR
#> [1]  51 100 150
Hausdorff.dist(BS_CPD_LR, cpt_true)
#> [1] 1

Perform wild binary segmentation with user specified threshold parameter

tau_WBS = 4 # user specified threshold parameter for WBS
intervals = WBS.intervals(M = 300, lower = 1, upper = n)
WBS_object = WBS.univar(y, 1, n, intervals$Alpha, intervals$Beta, delta)
WBS_result = thresholdBS(WBS_object, tau_WBS)
WBS_result
#> $BS_tree_trimmed
#>                                       levelName
#> 1 Binary Segmentation Tree Trimmed with tau = 4
#> 2  °--N100                                     
#> 3      ¦--N51                                  
#> 4      °--N150                                 
#> 
#> $cpt_hat
#>   location     value level
#> 1      100 17.456657     1
#> 2       51 12.785814     2
#> 3      150  9.520455     2
cpt_WBS_hat = sort(WBS_result$cpt_hat[,1]) # estimated change points by WBS
cpt_WBS_hat
#> [1]  51 100 150
Hausdorff.dist(cpt_WBS_hat, cpt_true)
#> [1] 1
cpt_WBSlr_hat = local.refine.univar(cpt_WBS_hat, y) # perform local refinement
cpt_WBSlr_hat
#> [1]  51 100 150
Hausdorff.dist(cpt_WBSlr_hat, cpt_true)
#> [1] 1

Perform wild binary segmentation with automatically tuned threshold parameter

WBS_CPD_result = tuneBSunivar(WBS_object, y)
WBS_CPD_result$cpt
#> [1]  51 100 150
Hausdorff.dist(WBS_CPD_result$cpt, cpt_true)
#> [1] 1
WBS_CPD_LR = local.refine.univar(WBS_CPD_result$cpt, y)
WBS_CPD_LR
#> [1]  51 100 150
Hausdorff.dist(WBS_CPD_LR, cpt_true)
#> [1] 1