SimulationsCrossvalidationCP/Simulation_ExampleOverestimation1e4Mod.R at main · FlorianPein/SimulationsCrossvalidationCP · 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
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138

setwd("~/Desktop/Florian-Rajen/Cross-validation/Simulations")
source("methods.R")
maxK <- 30

M <- 1e4
sigma <- 0.0001
seed <- 4e6 + 2e4

name <- "over1e4Mod"

n <- 202
Delta <- 1

signal <- c(rep(0, n / 2 + 1), rep(Delta, n / 2 - 1))

Y <- matrix(0, M, n)
for (i in 1:M) {
  if (i %% 10 == 1) {
    print(i)
  }

  set.seed(seed + i)
  Y[i, ] <- rnorm(n, signal, sigma)
}
R.matlab::writeMat(paste0("matData/", name, "/Data.mat"), y = Y)

#####
# run run.m in Matlab before continuing
#####

namesMethod <- c("arlotV2", "arlotV5")
abbr <- c("V2", "V5")
est <- list()
mse <- list()

est$arlotV2 <- as.numeric(R.matlab::readMat(paste0("matData/", name, "/numberV2.mat"))$numberjumpsV2)
est$arlotV5 <- as.numeric(R.matlab::readMat(paste0("matData/", name, "/numberV5.mat"))$numberjumpsV5)


for (i in seq_along(abbr)) {
  values <- R.matlab::readMat(paste0("matData/", name, "/values", abbr[i], ".mat"))[[paste0("values", abbr[i])]]
  jumps <- R.matlab::readMat(paste0("matData/", name, "/jumps", abbr[i], ".mat"))[[paste0("jumps", abbr[i])]]

  mse[[namesMethod[i]]] <- numeric(M)

  for (j in 1:M) {
    cps <- c(0, jumps[j, jumps[j, ] != 0], n)
    fit <- numeric(n)
    for (l in seq_along(cps)[-1]) {
      fit[(cps[l - 1] + 1):cps[l]] <- values[j, l - 1]
    }
    mse[[namesMethod[i]]][j] <- mean((signal - fit)^2)
  }
}

namesMethod2 <- c("copps", "CV1", "CVmod", "CV2fold", "CV5fold", "CV10fold", "CV20fold", "nonAdaptiveCV5fold",
                  "pelt", "binseg", "sbs", "wbs", "fdrseg", "smuce", "hsmuce", "msfpop", "huber", "biweight")
namesMethod <- c(namesMethod, namesMethod2)

q <- FDRSeg::simulQuantile(0.9, n = n, type = "fdrseg")

for (i in 1:M) {
  if (i %% 10 == 1) {
    print(i)
  }

  set.seed(seed + i)
  Y <- rnorm(n, signal, sigma)

  ret <- estimating(Y = Y, q = q, maxK = maxK, fit = TRUE)
  for (nameM in namesMethod2) {
    est[[nameM]][i] <- ret[[nameM]]
    mse[[nameM]][i] <- mean((signal - ret$fit[[nameM]])^2)
  }
}

saveRDS(est, file = paste("results/Est_", name, ".RDS", sep = ""))
saveRDS(mse, file = paste("results/Mse_", name, ".RDS", sep = ""))


K <- 1
library(xtable)
tab <- matrix(0, length(namesMethod), 4)
for (i in seq_along(namesMethod)) {
  namesM <- namesMethod[i]
  tab[i, ] <- c(sum(est[[namesM]] < K), sum(est[[namesM]] == K), sum(est[[namesM]] > K), sum(mse[[namesM]]) / 100)
}
tab <- tab / M * 100


# methods <- c("COPSS", "COPPS L1", "COPPS new", "2-fold CV", "5-fold CV", "10-fold CV", "20-fold CV", "LOOCV", "5-fold CV L2",
#              "PELT", "BinSeg", "SBS", "WBS", "FDRseg", "Smuce")
tab <- cbind(namesMethod, as.data.frame(tab))
colnames(tab) <- c("Method", "$\\hat{K} < K$","$\\hat{K} = K$", "$\\hat{K} > K$", "MISE")
tab <- xtable(tab, digits = rep(4, 6), auto = TRUE, display = c("g","s","g","g","g", "g"),
              align = c("l", "l|", "c", "c", "c", "c"),
              caption = "test")
saveRDS(tab, file = paste("results/", name, ".RDS", sep = ""))
print(tab, include.rownames = FALSE, file = paste("results/", name, ".tex", sep = ""),
      append = FALSE, sanitize.text.function=function(x){x})


######

setwd("~/Schreibtisch/cam/projects/crossvalidationCP/R")

tab1 <- readRDS(file = paste("results/ExampleOverestimationSd1.RDS", sep = ""))
tab2 <- readRDS(file = paste("results/ExampleOverestimationSd0p1.RDS", sep = ""))
tab3 <- readRDS(file = paste("results/ExampleOverestimationSd0p01.RDS", sep = ""))

tab <- cbind(tab1[, -1], tab2[, -1], tab3[, -1])
# methods <- c("COPSS", "COPPS L1", "COPPS new", "2-fold CV", "5-fold CV", "10-fold CV", "20-fold CV", "LOOCV",
#              "PELT", "BinSeg", "SBS", "WBS", "FDRseg", "Smuce")
# tab <- cbind(methods, as.data.frame(tab))
colnames(tab) <- c("Method", rep(c("$\\hat{K} < K$","$\\hat{K} = K$", "$\\hat{K} > K$"), 3))
tab <- xtable(tab, digits = rep(4, 2 + 3 * 3), auto = TRUE, display = c("g","s", rep(c("g","g","g"), 3)),
              align = c("l", "l|", rep(c("c", "c", "c"), 3)),
              caption = "Overeestimation example I.")
print(tab, include.rownames = FALSE, file = paste("results/ExampleOverestimation1.tex", sep = ""),
      append = FALSE, sanitize.text.function=function(x){x})


tab1 <- readRDS(file = paste("results/ExampleOverestimationSd0p001.RDS", sep = ""))
tab2 <- readRDS(file = paste("results/ExampleOverestimationSd0p0001.RDS", sep = ""))
tab3 <- readRDS(file = paste("results/ExampleOverestimationSd0p0001Even.RDS", sep = ""))

tab <- cbind(tab1[, -1], tab2[, -1], tab3[, -1])
methods <- c("COPSS", "COPPS L1", "COPPS new", "2-fold CV", "5-fold CV", "10-fold CV", "20-fold CV", "LOOCV",
             "PELT", "BinSeg", "SBS", "WBS", "FDRseg", "Smuce")
tab <- cbind(methods, as.data.frame(tab))
colnames(tab) <- c("Method", rep(c("$\\hat{K} < K$","$\\hat{K} = K$", "$\\hat{K} > K$"), 3))
tab <- xtable(tab, digits = rep(4, 2 + 3 * 3), auto = TRUE, display = c("g","s", rep(c("g","g","g"), 3)),
              align = c("l", "l|", rep(c("c", "c", "c"), 3)),
              caption = "Overeestimation example II.")
print(tab, include.rownames = FALSE, file = paste("results/ExampleOverestimation2.tex", sep = ""),
      append = FALSE, sanitize.text.function=function(x){x})