example_univariate_mean

library(changepoints)

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

There are 3 methods implemented for this purpose:

  1. DP.univar: performs dynamic programming for univariate mean change points detection.
  2. BS.univar: performs standard binary segmentation for univariate mean change points detection.
  3. WBS.univar: performs wild binary segmentation for univariate mean change points detection.

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