MostRecentChangepoint/Comp_time_corrections.R at master · bardwell/MostRecentChangepoint · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
# Case 1 - Differing valus of K
rm(list=ls())
library(compiler)
enableJIT(3)

library(tbart)
library(ecp)

# Simulate data
source("sim_data.R")

###################
# load in functions - for segmentation
###################
# DOUBLE CUSUM
source("bs_DC.R")
# INDEP
source("PELT.R")
source("ind.R")
# MRC
source("mrc_mean.R")
source("multiple_mrc.R")
# MV
source("PELT_MV.R")

n = 500
N = 100

############################################
#set.seed(1)
############################################

# K = 1
true.chpts = vector("list",10)
ind.chpts = vector("list",10)
dcbs.chpts = vector("list",10)
mrc.chpts = vector("list",10)
agg.chpts = vector("list",10)
ecp.chpts = vector("list",10)

comp_time = matrix(nrow=10 , ncol = 5)

for (i in 1:10){

  set.seed(i)
  sim = simulate.data( K=5 , eps=1.0 )
  data = sim$data

  # true MRC's each series
  true.chpts[[i]] = sim$mrc[ sim$series.mrc ]

  # AGG
  start_time <- Sys.time()
  agg = apply( data , 2 , sum )
  pagg = PELT( agg , 200*log(dim(data)[2]) )
  agg.chpts[[i]] = rep( rev( pagg$cpts )[1] , N )
  end_time <- Sys.time()
  comp_time[i,1] = end_time - start_time

  # ECP
  start_time <- Sys.time()
  ecp.div = e.divisive( t(data) , alpha = 1)
  ecp.chpts[[i]] = rep( rev( ecp.div$estimates )[2] , N )
  end_time <- Sys.time()
  comp_time[i,2] = end_time - start_time

  # IND
  start_time <- Sys.time()
  ind.chpts[[i]] = ind( data )
  end_time <- Sys.time()
  comp_time[i,3] = end_time - start_time

  # DCBS
  start_time <- Sys.time()
  dcbs.chpts[[i]] = Bin_seg(data,10)
  end_time <- Sys.time()
  comp_time[i,4] = end_time - start_time

  # MRC
  start_time <- Sys.time()
  mrc = mrc.mean( data , beta = 1.5*log(n) )
  c = multiple.mrc( mrc , pmax=10 )
  p.hat = c$MDL
  mrc.chpts[[i]] = c$locs[[p.hat]][ c$affected[[p.hat]] ]
  end_time <- Sys.time()
  comp_time[i,5] = end_time - start_time

  print(i)
}

# agg , mv ecp
times = apply( comp_time , 2 , mean)

names(times) = c("AGG","ECP","IND","DCBS","MRC")
t( data.frame(times) )
library(xtable)
xtable(t( times ) )

## comp time -- MRC
comp_time = numeric(15)

for (i in 1:15){
print(i)
  set.seed(1)
  sim = simulate.data( K=5 , eps=1.0 )
  data = sim$data

  # MRC
  start_time <- Sys.time()
  mrc = mrc.mean( data , beta = 1.5*log(n) )
  c = multiple.mrc( mrc , pmax=i )
  p.hat = c$MDL
  end_time <- Sys.time()
  comp_time[i] = end_time - start_time

}

pdf("Kmax_comp_cosr.pdf")
plot( 2:15 ,  comp_time[2:15]  , ylab = "Computational cost (seconds)" , xlab = expression(K[max]) )
dev.off()