Update fastcpd 0.14.0

*   Correct mBIC.

DESCRIPTION

@@ -1,7 +1,7 @@
 Type: Package
 Package: fastcpd
 Title: Fast Change Point Detection via Sequential Gradient Descent
-Version: 0.13.2
+Version: 0.14.0
 Authors@R: c(
     person("Xingchi", "Li", , "[email protected]",
            role = c("aut", "cre", "cph"),

NEWS.md

@@ -1,3 +1,7 @@
+# fastcpd 0.14.0
+
+*   Correct mBIC.
+
 # fastcpd 0.13.2
 
 *   Remove package check in examples.

R/check.R

@@ -269,9 +269,9 @@ get_pruning_coef <- function(
     pruning_coef <- -Inf
   }
   if (!pruning_coef_is_set && cost_adjustment == "MBIC") {
-    pruning_coef <- pruning_coef + p * log(4 / n) / 2
+    pruning_coef <- pruning_coef + p * log(2)
   } else if (!pruning_coef_is_set && cost_adjustment == "MDL") {
-    pruning_coef <- pruning_coef + p * log2(4 / n) / 2
+    pruning_coef <- pruning_coef + p * log2(2)
   }
   pruning_coef
 }

R/fastcpd_wrappers.R

@@ -263,7 +263,8 @@ fastcpd.mean <- fastcpd_mean  # nolint: Conventional R function style
 #' function is similar to [fastcpd()] except that the data is by
 #' default a matrix or data frame or a vector with each row / element as an
 #' observation and thus a formula is not required here.
-#' @example tests/testthat/examples/fastcpd_meanvariance.R
+#' @example tests/testthat/examples/fastcpd_meanvariance_1.R
+#' @example tests/testthat/examples/fastcpd_meanvariance_2.R
 #' @seealso [fastcpd()]
 #'
 #' @md

cran-comments.md

@@ -1,5 +1,9 @@
 ## Updates since last CRAN release (0.13.1)
 
+### fastcpd 0.14.0
+
+*   Correct mBIC.
+
 ### fastcpd 0.13.2
 
 *   Remove package check in examples.

src/fastcpd_class_private.cc

@@ -141,7 +141,7 @@ void Fastcpd::create_segment_statistics() {
 double Fastcpd::get_cost_adjustment_value(const unsigned nrows) {
   double adjusted = 0;
   if (cost_adjustment == "MBIC" || cost_adjustment == "MDL") {
-    adjusted = p * std::log(nrows) / 2.0;
+    adjusted = p * std::log((double) nrows / data_n_rows) / 2.0;
   }
   if (cost_adjustment == "MDL") {
     adjusted *= std::log2(M_E);
@@ -517,7 +517,7 @@ CostResult Fastcpd::get_nll_meanvariance(
   const unsigned int segment_start,
   const unsigned int segment_end
 ) {
-  const rowvec data_diff =
+  rowvec data_diff =
     zero_data.row(segment_end + 1) - zero_data.row(segment_start);
   const unsigned int segment_length = segment_end - segment_start + 1;
 
@@ -526,8 +526,26 @@ CostResult Fastcpd::get_nll_meanvariance(
       data_diff.subvec(0, d - 1)).t() * (data_diff.subvec(0, d - 1)
     ) / segment_length
   ) / segment_length);
-  if (det_value <= 0) {
-    det_value = 1e-10;
+  if (segment_length <= d) {
+    unsigned int approximate_segment_start;
+    unsigned int approximate_segment_end;
+    if (segment_start >= d) {
+      approximate_segment_start = segment_start - d;
+    } else {
+      approximate_segment_start = 0;
+    }
+    if (segment_end < data_n_rows - d) {
+      approximate_segment_end = segment_end + d;
+    } else {
+      approximate_segment_end = data_n_rows - 1;
+    }
+    data_diff = zero_data.row(approximate_segment_end + 1) -
+      zero_data.row(approximate_segment_start);
+    det_value = det((
+    reshape(data_diff.subvec(d, p - 1), d, d) - (
+      data_diff.subvec(0, d - 1)).t() * (data_diff.subvec(0, d - 1)
+    ) / (approximate_segment_end - approximate_segment_start + 1)
+  ) / (approximate_segment_end - approximate_segment_start + 1));
   }
 
   return {
@@ -572,8 +590,23 @@ CostResult Fastcpd::get_nll_variance(
   double det_value = det(arma::reshape(
     zero_data.row(segment_end + 1) - zero_data.row(segment_start), d, d
   ) / segment_length);
-  if (det_value <= 0) {
-    det_value = 1e-10;
+  if (segment_length < d) {
+    unsigned int approximate_segment_start;
+    unsigned int approximate_segment_end;
+    if (segment_start >= d) {
+      approximate_segment_start = segment_start - d;
+    } else {
+      approximate_segment_start = 0;
+    }
+    if (segment_end < data_n_rows - d) {
+      approximate_segment_end = segment_end + d;
+    } else {
+      approximate_segment_end = data_n_rows - 1;
+    }
+    det_value = det(arma::reshape(
+      zero_data.row(approximate_segment_end + 1) -
+        zero_data.row(approximate_segment_start), d, d
+    ) / (approximate_segment_end - approximate_segment_start + 1));
   }
 
   return {

tests/testthat/examples/fastcpd_meanvariance_1.R

@@ -0,0 +1,13 @@
+set.seed(1)
+p <- 1
+result <- fastcpd.mv(
+  rbind(
+    mvtnorm::rmvnorm(300, mean = rep(0, p), sigma = diag(1, p)),
+    mvtnorm::rmvnorm(400, mean = rep(10, p), sigma = diag(1, p)),
+    mvtnorm::rmvnorm(300, mean = rep(0, p), sigma = diag(100, p)),
+    mvtnorm::rmvnorm(300, mean = rep(0, p), sigma = diag(1, p)),
+    mvtnorm::rmvnorm(400, mean = rep(10, p), sigma = diag(1, p)),
+    mvtnorm::rmvnorm(300, mean = rep(10, p), sigma = diag(100, p))
+  )
+)
+summary(result)

tests/testthat/examples/fastcpd_meanvariance_2.R

@@ -0,0 +1,13 @@
+set.seed(1)
+p <- 4
+result <- fastcpd.mv(
+  rbind(
+    mvtnorm::rmvnorm(300, mean = rep(0, p), sigma = diag(1, p)),
+    mvtnorm::rmvnorm(400, mean = rep(10, p), sigma = diag(1, p)),
+    mvtnorm::rmvnorm(300, mean = rep(0, p), sigma = diag(100, p)),
+    mvtnorm::rmvnorm(300, mean = rep(0, p), sigma = diag(1, p)),
+    mvtnorm::rmvnorm(400, mean = rep(10, p), sigma = diag(1, p)),
+    mvtnorm::rmvnorm(300, mean = rep(10, p), sigma = diag(100, p))
+  )
+)
+summary(result)

tests/testthat/examples/fastcpd_ts.txt

@@ -18,7 +18,7 @@ result <- fastcpd.ts(
   lower = c(-2, -2, -2, -2, 1e-10),
   upper = c(2, 2, 2, 2, Inf),
   line_search = c(1, 0.1, 1e-2),
-  trim = 0.04
+  trim = 0.05
 )
 summary(result)
 plot(result)

tests/testthat/test-coverage.R

@@ -34,103 +34,6 @@ testthat::test_that(
   }
 )
 
-testthat::test_that(
-  "1d mv", {
-    set.seed(1)
-    data <- c(
-      rnorm(300, 0, 1), rnorm(400, 10, 1), rnorm(300, 0, 50),
-      rnorm(300, 0, 1), rnorm(400, 10, 1), rnorm(300, 10, 50)
-    )
-    result_meanvariance <- fastcpd.meanvariance(data)
-    testthat::expect_equal(
-      result_meanvariance@cp_set, c(300, 700, 1000, 1300, 1700)
-    )
-  }
-)
-
-testthat::test_that(
-  "1d custom multiple epochs", {
-    set.seed(1)
-    p <- 1
-    x <- matrix(rnorm(300 * p, 0, 1), ncol = p)
-    theta <- rbind(rnorm(p, 0, 1), rnorm(p, 4, 1))
-    y <- c(
-      rbinom(150, 1, 1 / (1 + exp(-x[1:150, ] * theta[1, ]))),
-      rbinom(150, 1, 1 / (1 + exp(-x[151:300, ] * theta[2, ])))
-    )
-    logistic_loss <- function(data, theta) {
-      x <- data[, -1]
-      y <- data[, 1]
-      u <- x * c(theta)
-      nll <- -y * u + log(1 + exp(u))
-      nll[u > 10] <- -y[u > 10] * u[u > 10] + u[u > 10]
-      sum(nll)
-    }
-    logistic_loss_gradient <- function(data, theta) {
-      x <- data[nrow(data), -1]
-      y <- data[nrow(data), 1]
-      c(-(y - 1 / (1 + exp(-x * theta)))) * x
-    }
-    logistic_loss_hessian <- function(data, theta) {
-      x <- data[nrow(data), -1]
-      prob <- 1 / (1 + exp(-x * theta))
-      (x %o% x) * c((1 - prob) * prob)
-    }
-    result <- fastcpd(
-      formula = y ~ . - 1,
-      data = data.frame(y = y, x = x),
-      beta = "BIC",
-      cost = logistic_loss,
-      cost_gradient = logistic_loss_gradient,
-      cost_hessian = logistic_loss_hessian,
-      multiple_epochs = function(segment_length) {
-        if (segment_length < 10) 1 else 0
-      },
-      r.progress = FALSE
-    )
-    testthat::expect_equal(result@cp_set, 137)
-  }
-)
-
-testthat::test_that(
-  "2d custom", {
-    set.seed(1)
-    p <- 2
-    x <- matrix(rnorm(300 * p, 0, 1), ncol = p)
-    theta <- rbind(rnorm(p, 0, 1), rnorm(p, 4, 1))
-    y <- c(
-      rbinom(150, 1, 1 / (1 + exp(-x[1:150, ] * theta[1, ]))),
-      rbinom(150, 1, 1 / (1 + exp(-x[151:300, ] * theta[2, ])))
-    )
-    logistic_loss <- function(data, theta) {
-      x <- data[, -1]
-      y <- data[, 1]
-      u <- x %*% theta
-      nll <- -y * u + log(1 + exp(u))
-      nll[u > 10] <- -y[u > 10] * u[u > 10] + u[u > 10]
-      sum(nll)
-    }
-    logistic_loss_gradient <- function(data, theta) {
-      x <- data[nrow(data), -1]
-      y <- data[nrow(data), 1]
-      c(-(y - 1 / (1 + exp(-x %*% theta)))) * x
-    }
-    logistic_loss_hessian <- function(data, theta) {
-      x <- data[nrow(data), -1]
-      prob <- 1 / (1 + exp(-x %*% theta))
-      (x %o% x) * c((1 - prob) * prob)
-    }
-    result <- fastcpd(
-      formula = y ~ . - 1,
-      data = data.frame(y = y, x = x),
-      cost = logistic_loss,
-      cost_gradient = logistic_loss_gradient,
-      cost_hessian = logistic_loss_hessian
-    )
-    testthat::expect_equal(result@cp_set, 153)
-  }
-)
-
 testthat::test_that(
   "warm start", {
     set.seed(1)

tests/testthat/test-examples-data-transcriptome.R

@@ -12,8 +12,8 @@ testthat::test_that(
     testthat::expect_equal(
       result@cp_set,
       c(
-        177, 264, 393, 534, 578, 656, 788, 811, 869, 934, 960, 1051, 1141, 1286,
-        1319, 1367, 1566, 1657, 1724, 1906, 1972, 1996, 2041, 2141, 2200
+        178, 264, 401, 534, 601, 656, 788, 811, 869, 934, 971, 1055,
+        1142, 1286, 1319, 1386, 1657, 1724, 1906, 1972, 1996, 2041
       )
     )
   }

tests/testthat/test-examples-data-well_log.R

@@ -6,9 +6,9 @@ testthat::test_that(
     testthat::expect_equal(
       result@cp_set,
       c(
-        7, 19, 356, 445, 716, 792, 1034, 1070, 1215, 1368,
-        1428, 1526, 1684, 1866, 2047, 2409, 2469, 2531, 2591,
-        2775, 3490, 3533, 3673, 3744, 3855, 3886, 3945, 3963, 4035
+        7, 19, 65, 356, 445, 717, 792, 1034, 1070, 1215, 1368, 1428,
+        1526, 1684, 1866, 2047, 2409, 2469, 2531, 2591, 2775, 3166,
+        3314, 3490, 3533, 3673, 3744, 3855, 3886, 3945, 3963, 4035
       )
     )
   }

tests/testthat/test-examples-fastcpd_2.R

@@ -4,6 +4,6 @@ testthat::test_that(
     testthat::skip_if_not_installed("mvtnorm")
 
     source("examples/fastcpd_2.R")
-    testthat::expect_equal(huber_regression_result@cp_set, c(401, 726))
+    testthat::expect_equal(huber_regression_result@cp_set, c(418, 726))
   }
 )

tests/testthat/test-examples-fastcpd_3.R

@@ -9,8 +9,8 @@ testthat::test_that(
       collapse = "\n"
     )))
 
-    testthat::expect_equal(result_builtin@cp_set, 200)
-    testthat::expect_equal(result_custom@cp_set, 199)
-    testthat::expect_equal(result_custom_two_epochs@cp_set, 200)
+    testthat::expect_equal(result_builtin@cp_set, 202)
+    testthat::expect_equal(result_custom@cp_set, 202)
+    testthat::expect_equal(result_custom_two_epochs@cp_set, 202)
   }
 )

tests/testthat/test-examples-fastcpd_arima.R

@@ -8,6 +8,6 @@ testthat::test_that(
       collapse = "\n"
     )))
 
-    testthat::expect_equal(result@cp_set, 178)
+    testthat::expect_equal(result@cp_set, 184)
   }
 )

tests/testthat/test-examples-fastcpd_binomial.R

@@ -3,6 +3,6 @@ testthat::test_that(
     testthat::skip_if_not_installed("ggplot2")
 
     source("examples/fastcpd_binomial.R")
-    testthat::expect_equal(result@cp_set, 300)
+    testthat::expect_equal(result@cp_set, 302)
   }
 )

tests/testthat/test-examples-fastcpd_lm.R

@@ -4,7 +4,7 @@ testthat::test_that(
     testthat::skip_if_not_installed("mvtnorm")
 
     source("examples/fastcpd_lm.R")
-    testthat::expect_equal(result_lm@cp_set, c(99, 201))
+    testthat::expect_equal(result_lm@cp_set, c(97, 201))
     testthat::expect_equal(result_mlm@cp_set, 350)
   }
 )

tests/testthat/test-examples-fastcpd_mean-time_1.R

@@ -1,7 +1,7 @@
 testthat::test_that(
   "examples/fastcpd_mean-time_1.R", {
     source("examples/fastcpd_mean-time_1.R")
-    testthat::expect_equal(result@cp_set, 10006)
+    testthat::expect_equal(result@cp_set, 10007)
     testthat::expect_lt(result_time[3], 30)
   }
 )