compareKernelSmoothing <- function(y, bandwidth) { n <- length(y) kernel <- function(x) ifelse(abs(x) <= 1 + 1e-15, 3 / 4 * (1 - x^2), 0) L <- max(0, min(as.integer(bandwidth * n + 1e-14), n - 1)) K <- kernel(-L:L / (bandwidth * n)) ret <- numeric(n) for (i in 1:n) { ret[i] <- sum(K[max(L + 2 - i, 1):min(n - i + L + 1, 2 * L + 1)] * y[max(i - L, 1):min(i + L, n)]) / sum(K[max(L + 2 - i, 1):min(n - i + L + 1, 2 * L + 1)]) } ret } compareImSy <- function(y, bandwidth) { n <- length(y) kernel <- function(x) ifelse(abs(x) <= 1 + 1e-15, 3 / 4 * (1 - x^2), 0) L <- max(0, min(as.integer(bandwidth * n + 1e-14), n - 1)) K <- kernel(-L:L / (bandwidth * n)) y - .kernelSmoothing(y, K) } compareImSX <- function(y, bandwidth) { n <- length(y) kernel <- function(x) ifelse(abs(x) <= 1 + 1e-15, 3 / 4 * (1 - x^2), 0) L <- max(0, min(as.integer(bandwidth * n + 1e-14), n - 1)) K <- kernel(-L:L / (bandwidth * n)) Matrix::bandSparse(n, n - 1, k = (-L + 1):L - 1, diagonals = .createImSX(n = n, K = rev(cumsum(K)))) } compareXty <- function(y, bandwidth, x = NULL) { if (is.null(x)) { x <- compareImSX(y, bandwidth) } Matrix::t(x) %*% y / length(y) } compareXtX <- function(y, bandwidth, x = NULL) { if (is.null(x)) { x <- compareImSX(y, bandwidth) } Matrix::t(x) %*% x / length(y) } test_that("basic functions are giving correct results I", { testy <- rnorm(100) testbandwidth <- 0.05 testn <- length(testy) testkernel <- function(x) ifelse(abs(x) <= 1 + 1e-15, 3 / 4 * (1 - x^2), 0) testL <- max(0, min(as.integer(testbandwidth * testn + 1e-14), testn - 1)) testK <- testkernel(-testL:testL / (testbandwidth * testn)) expect_equal(.kernelSmoothing(testy, testK), compareKernelSmoothing(testy, testbandwidth)) testImSy <- testy - .kernelSmoothingEpanechnikov(testy, testbandwidth) expect_equal(testImSy, compareImSy(testy, testbandwidth), tolerance = 1e-6) testnull_dev <- mean(testImSy^2) testlistKernel <- list(cusumKernel = numeric(2 * testL + 1), Xty = numeric(testn - 1), XtX = matrix(0, 2 * testL - 1, 4 * testL - 2), isComputedXtX = matrix(FALSE, 2 * testL - 1, 4 * testL - 2), XtXgap = numeric(2 * testL), ImSX = matrix(0, 3 * testL, 2 * testL), isComputedImSX = logical(2 * testL)) .initializeKernel(testImSy, testbandwidth, testlistKernel$cusumKernel, testlistKernel$Xty, testlistKernel$XtX, testlistKernel$isComputedXtX, testlistKernel$XtXgap, testlistKernel$ImSX, testlistKernel$isComputedImSX) object <- .getXty(testImSy, testbandwidth, testlistKernel$cusumKernel, testlistKernel$Xty, testlistKernel$XtX, testlistKernel$isComputedXtX, testlistKernel$XtXgap, testlistKernel$ImSX, testlistKernel$isComputedImSX) expect_equal(object, as.numeric(compareXty(testy, testbandwidth)), tolerance = 1e-1) expect_equal(.getXtX(testImSy, testbandwidth), as.matrix(compareXtX(testy, testbandwidth)), tolerance = 1e-6) testthresh <- 1e-7 testmaxit <- 1e5 testlambda <- 0.01 object <- as.numeric(.lassoImSX(testImSy, testbandwidth, testlistKernel$cusumKernel, testlistKernel$Xty, testlistKernel$XtX, testlistKernel$isComputedXtX, testlistKernel$XtXgap, testlistKernel$ImSX, testlistKernel$isComputedImSX, testlambda, testthresh * testnull_dev, testmaxit)) testImSX <- compareImSX(testy, testbandwidth) expected <- as.numeric(glmnet::glmnet(x = testImSX, y = testImSy, lambda = testlambda, standardize = FALSE, intercept = FALSE, thresh = testthresh, maxit = testmaxit)$beta) expect_equal(object, expected, tolerance = 1e-2) testcps <- c(5, 20, 56, 95) testImSX <- cbind(rep(0, testn), testImSX) object <- .postProcessing(testcps - 2, testImSy, testbandwidth, testlistKernel$cusumKernel, testlistKernel$Xty, testlistKernel$XtX, testlistKernel$isComputedXtX, testlistKernel$XtXgap, testlistKernel$ImSX, testlistKernel$isComputedImSX) expected <- as.numeric(Matrix::solve(Matrix::t(testImSX[, testcps, drop = FALSE]) %*% testImSX[, testcps, drop = FALSE], as.numeric(Matrix::t(testImSX[, testcps, drop = FALSE]) %*% testImSy))) expect_equal(object, expected, tolerance = 1e-2) }) test_that("basic functions are giving correct results II", { testy <- rnorm(100) testbandwidth <- 0.15 testn <- length(testy) testkernel <- function(x) ifelse(abs(x) <= 1 + 1e-15, 3 / 4 * (1 - x^2), 0) testL <- max(0, min(as.integer(testbandwidth * testn + 1e-14), testn - 1)) testK <- testkernel(-testL:testL / (testbandwidth * testn)) expect_equal(.kernelSmoothing(testy, testK), compareKernelSmoothing(testy, testbandwidth)) testImSy <- testy - .kernelSmoothingEpanechnikov(testy, testbandwidth) expect_equal(testImSy, compareImSy(testy, testbandwidth), tolerance = 1e-6) testnull_dev <- mean(testImSy^2) testlistKernel <- list(cusumKernel = numeric(2 * testL + 1), Xty = numeric(testn - 1), XtX = matrix(0, 2 * testL - 1, 4 * testL - 2), isComputedXtX = matrix(FALSE, 2 * testL - 1, 4 * testL - 2), XtXgap = numeric(2 * testL), ImSX = matrix(0, 3 * testL, 2 * testL), isComputedImSX = logical(2 * testL)) .initializeKernel(testImSy, testbandwidth, testlistKernel$cusumKernel, testlistKernel$Xty, testlistKernel$XtX, testlistKernel$isComputedXtX, testlistKernel$XtXgap, testlistKernel$ImSX, testlistKernel$isComputedImSX) object <- .getXty(testImSy, testbandwidth, testlistKernel$cusumKernel, testlistKernel$Xty, testlistKernel$XtX, testlistKernel$isComputedXtX, testlistKernel$XtXgap, testlistKernel$ImSX, testlistKernel$isComputedImSX) expect_equal(object, as.numeric(compareXty(testy, testbandwidth)), tolerance = 1e-1) expect_equal(.getXtX(testImSy, testbandwidth), as.matrix(compareXtX(testy, testbandwidth)), tolerance = 1e-6) testthresh <- 1e-7 testmaxit <- 1e5 testlambda <- 0.01 object <- as.numeric(.lassoImSX(testImSy, testbandwidth, testlistKernel$cusumKernel, testlistKernel$Xty, testlistKernel$XtX, testlistKernel$isComputedXtX, testlistKernel$XtXgap, testlistKernel$ImSX, testlistKernel$isComputedImSX, testlambda, testthresh * testnull_dev, testmaxit)) testImSX <- compareImSX(testy, testbandwidth) expected <- as.numeric(glmnet::glmnet(x = testImSX, y = testImSy, lambda = testlambda, standardize = FALSE, intercept = FALSE, thresh = testthresh, maxit = testmaxit)$beta) expect_equal(object, expected, tolerance = 1e-1) testcps <- c(5, 11) testImSX <- cbind(rep(0, testn), testImSX) object <- .postProcessing(testcps - 2, testImSy, testbandwidth, testlistKernel$cusumKernel, testlistKernel$Xty, testlistKernel$XtX, testlistKernel$isComputedXtX, testlistKernel$XtXgap, testlistKernel$ImSX, testlistKernel$isComputedImSX) expected <- as.numeric(Matrix::solve(Matrix::t(testImSX[, testcps, drop = FALSE]) %*% testImSX[, testcps, drop = FALSE], as.numeric(Matrix::t(testImSX[, testcps, drop = FALSE]) %*% testImSy))) expect_equal(object, expected, tolerance = 1e-1) }) test_that("basic functions are giving correct results III", { testy <- rnorm(100) testbandwidth <- 0.01 testn <- length(testy) testkernel <- function(x) ifelse(abs(x) <= 1 + 1e-15, 3 / 4 * (1 - x^2), 0) testL <- max(0, min(as.integer(testbandwidth * testn + 1e-14), testn - 1)) testK <- testkernel(-testL:testL / (testbandwidth * testn)) expect_equal(.kernelSmoothing(testy, testK), compareKernelSmoothing(testy, testbandwidth)) testImSy <- testy - .kernelSmoothingEpanechnikov(testy, testbandwidth) expect_equal(testImSy, compareImSy(testy, testbandwidth), tolerance = 1e-6) testnull_dev <- mean(testImSy^2) testlistKernel <- list(cusumKernel = numeric(2 * testL + 1), Xty = numeric(testn - 1), XtX = matrix(0, 2 * testL - 1, 4 * testL - 2), isComputedXtX = matrix(FALSE, 2 * testL - 1, 4 * testL - 2), XtXgap = numeric(2 * testL), ImSX = matrix(0, 3 * testL, 2 * testL), isComputedImSX = logical(2 * testL)) .initializeKernel(testImSy, testbandwidth, testlistKernel$cusumKernel, testlistKernel$Xty, testlistKernel$XtX, testlistKernel$isComputedXtX, testlistKernel$XtXgap, testlistKernel$ImSX, testlistKernel$isComputedImSX) object <- .getXty(testImSy, testbandwidth, testlistKernel$cusumKernel, testlistKernel$Xty, testlistKernel$XtX, testlistKernel$isComputedXtX, testlistKernel$XtXgap, testlistKernel$ImSX, testlistKernel$isComputedImSX) expect_equal(object, as.numeric(compareXty(testy, testbandwidth)), tolerance = 1e-1) expect_equal(.getXtX(testImSy, testbandwidth), as.matrix(compareXtX(testy, testbandwidth)), tolerance = 1e-6) }) test_that("basic functions are giving correct results IV", { testy <- rnorm(100) testbandwidth <- 0.24 testn <- length(testy) testkernel <- function(x) ifelse(abs(x) <= 1 + 1e-15, 3 / 4 * (1 - x^2), 0) testL <- max(0, min(as.integer(testbandwidth * testn + 1e-14), testn - 1)) testK <- testkernel(-testL:testL / (testbandwidth * testn)) expect_equal(.kernelSmoothing(testy, testK), compareKernelSmoothing(testy, testbandwidth)) testImSy <- testy - .kernelSmoothingEpanechnikov(testy, testbandwidth) expect_equal(testImSy, compareImSy(testy, testbandwidth), tolerance = 1e-1) testnull_dev <- mean(testImSy^2) testlistKernel <- list(cusumKernel = numeric(2 * testL + 1), Xty = numeric(testn - 1), XtX = matrix(0, 2 * testL - 1, 4 * testL - 2), isComputedXtX = matrix(FALSE, 2 * testL - 1, 4 * testL - 2), XtXgap = numeric(2 * testL), ImSX = matrix(0, 3 * testL, 2 * testL), isComputedImSX = logical(2 * testL)) .initializeKernel(testImSy, testbandwidth, testlistKernel$cusumKernel, testlistKernel$Xty, testlistKernel$XtX, testlistKernel$isComputedXtX, testlistKernel$XtXgap, testlistKernel$ImSX, testlistKernel$isComputedImSX) object <- .getXty(testImSy, testbandwidth, testlistKernel$cusumKernel, testlistKernel$Xty, testlistKernel$XtX, testlistKernel$isComputedXtX, testlistKernel$XtXgap, testlistKernel$ImSX, testlistKernel$isComputedImSX) expect_equal(object, as.numeric(compareXty(testy, testbandwidth)), tolerance = 1e-1) expect_equal(.getXtX(testImSy, testbandwidth), as.matrix(compareXtX(testy, testbandwidth)), tolerance = 1e-1) testthresh <- 1e-7 testmaxit <- 1e5 testlambda <- 0.01 object <- as.numeric(.lassoImSX(testImSy, testbandwidth, testlistKernel$cusumKernel, testlistKernel$Xty, testlistKernel$XtX, testlistKernel$isComputedXtX, testlistKernel$XtXgap, testlistKernel$ImSX, testlistKernel$isComputedImSX, testlambda, testthresh * testnull_dev, testmaxit)) testImSX <- compareImSX(testy, testbandwidth) expected <- as.numeric(glmnet::glmnet(x = testImSX, y = testImSy, lambda = testlambda, standardize = FALSE, intercept = FALSE, thresh = testthresh, maxit = testmaxit)$beta) expect_equal(object, expected, tolerance = 1e-1) testcps <- c(5, 95) testImSX <- cbind(rep(0, testn), testImSX) object <- .postProcessing(testcps - 2, testImSy, testbandwidth, testlistKernel$cusumKernel, testlistKernel$Xty, testlistKernel$XtX, testlistKernel$isComputedXtX, testlistKernel$XtXgap, testlistKernel$ImSX, testlistKernel$isComputedImSX) expected <- as.numeric(Matrix::solve(Matrix::t(testImSX[, testcps, drop = FALSE]) %*% testImSX[, testcps, drop = FALSE], as.numeric(Matrix::t(testImSX[, testcps, drop = FALSE]) %*% testImSy))) expect_equal(object, expected, tolerance = 1e-1) }) test_that("basic functions are giving correct results V", { testy <- 1:100 testbandwidth <- 0.075 testn <- length(testy) testkernel <- function(x) ifelse(abs(x) <= 1 + 1e-15, 3 / 4 * (1 - x^2), 0) testL <- max(0, min(as.integer(testbandwidth * testn + 1e-14), testn - 1)) testK <- testkernel(-testL:testL / (testbandwidth * testn)) expect_equal(.kernelSmoothing(testy, testK), compareKernelSmoothing(testy, testbandwidth)) testImSy <- testy - .kernelSmoothingEpanechnikov(testy, testbandwidth) expect_equal(testImSy, compareImSy(testy, testbandwidth), tolerance = 1e-6) testnull_dev <- mean(testImSy^2) testlistKernel <- list(cusumKernel = numeric(2 * testL + 1), Xty = numeric(testn - 1), XtX = matrix(0, 2 * testL - 1, 4 * testL - 2), isComputedXtX = matrix(FALSE, 2 * testL - 1, 4 * testL - 2), XtXgap = numeric(2 * testL), ImSX = matrix(0, 3 * testL, 2 * testL), isComputedImSX = logical(2 * testL)) .initializeKernel(testImSy, testbandwidth, testlistKernel$cusumKernel, testlistKernel$Xty, testlistKernel$XtX, testlistKernel$isComputedXtX, testlistKernel$XtXgap, testlistKernel$ImSX, testlistKernel$isComputedImSX) object <- .getXty(testImSy, testbandwidth, testlistKernel$cusumKernel, testlistKernel$Xty, testlistKernel$XtX, testlistKernel$isComputedXtX, testlistKernel$XtXgap, testlistKernel$ImSX, testlistKernel$isComputedImSX) expect_equal(object, as.numeric(compareXty(testy, testbandwidth)), tolerance = 1e-1) expect_equal(.getXtX(testImSy, testbandwidth), as.matrix(compareXtX(testy, testbandwidth)), tolerance = 1e-6) testthresh <- 1e-7 testmaxit <- 1e5 testlambda <- 0.01 object <- as.numeric(.lassoImSX(testImSy, testbandwidth, testlistKernel$cusumKernel, testlistKernel$Xty, testlistKernel$XtX, testlistKernel$isComputedXtX, testlistKernel$XtXgap, testlistKernel$ImSX, testlistKernel$isComputedImSX, testlambda, testthresh * testnull_dev, testmaxit)) testImSX <- compareImSX(testy, testbandwidth) expected <- as.numeric(glmnet::glmnet(x = testImSX, y = testImSy, lambda = testlambda, standardize = FALSE, intercept = FALSE, thresh = testthresh, maxit = testmaxit)$beta) expect_equal(object, expected, tolerance = 1e-1) testcps <- c(5) testImSX <- cbind(rep(0, testn), testImSX) object <- .postProcessing(testcps - 2, testImSy, testbandwidth, testlistKernel$cusumKernel, testlistKernel$Xty, testlistKernel$XtX, testlistKernel$isComputedXtX, testlistKernel$XtXgap, testlistKernel$ImSX, testlistKernel$isComputedImSX) expected <- as.numeric(Matrix::solve(Matrix::t(testImSX[, testcps, drop = FALSE]) %*% testImSX[, testcps, drop = FALSE], as.numeric(Matrix::t(testImSX[, testcps, drop = FALSE]) %*% testImSy))) expect_equal(object, expected, tolerance = 1e-1) }) test_that("basic functions are giving correct results VI", { testy <- 100:1 testbandwidth <- 0.0346 testn <- length(testy) testkernel <- function(x) ifelse(abs(x) <= 1 + 1e-15, 3 / 4 * (1 - x^2), 0) testL <- max(0, min(as.integer(testbandwidth * testn + 1e-14), testn - 1)) testK <- testkernel(-testL:testL / (testbandwidth * testn)) expect_equal(.kernelSmoothing(testy, testK), compareKernelSmoothing(testy, testbandwidth)) testImSy <- testy - .kernelSmoothingEpanechnikov(testy, testbandwidth) expect_equal(testImSy, compareImSy(testy, testbandwidth), tolerance = 1e-6) testnull_dev <- mean(testImSy^2) testlistKernel <- list(cusumKernel = numeric(2 * testL + 1), Xty = numeric(testn - 1), XtX = matrix(0, 2 * testL - 1, 4 * testL - 2), isComputedXtX = matrix(FALSE, 2 * testL - 1, 4 * testL - 2), XtXgap = numeric(2 * testL), ImSX = matrix(0, 3 * testL, 2 * testL), isComputedImSX = logical(2 * testL)) .initializeKernel(testImSy, testbandwidth, testlistKernel$cusumKernel, testlistKernel$Xty, testlistKernel$XtX, testlistKernel$isComputedXtX, testlistKernel$XtXgap, testlistKernel$ImSX, testlistKernel$isComputedImSX) object <- .getXty(testImSy, testbandwidth, testlistKernel$cusumKernel, testlistKernel$Xty, testlistKernel$XtX, testlistKernel$isComputedXtX, testlistKernel$XtXgap, testlistKernel$ImSX, testlistKernel$isComputedImSX) expect_equal(object, as.numeric(compareXty(testy, testbandwidth)), tolerance = 1e-1) expect_equal(.getXtX(testImSy, testbandwidth), as.matrix(compareXtX(testy, testbandwidth)), tolerance = 1e-6) testthresh <- 1e-7 testmaxit <- 1e5 testlambda <- 0.01 object <- as.numeric(.lassoImSX(testImSy, testbandwidth, testlistKernel$cusumKernel, testlistKernel$Xty, testlistKernel$XtX, testlistKernel$isComputedXtX, testlistKernel$XtXgap, testlistKernel$ImSX, testlistKernel$isComputedImSX, testlambda, testthresh * testnull_dev, testmaxit)) testImSX <- compareImSX(testy, testbandwidth) expected <- as.numeric(glmnet::glmnet(x = testImSX, y = testImSy, lambda = testlambda, standardize = FALSE, intercept = FALSE, thresh = testthresh, maxit = testmaxit)$beta) expect_equal(object, expected, tolerance = 1e-1) testcps <- c(5, 6, 20, 21, 56, 95, 96) testImSX <- cbind(rep(0, testn), testImSX) object <- .postProcessing(testcps - 2, testImSy, testbandwidth, testlistKernel$cusumKernel, testlistKernel$Xty, testlistKernel$XtX, testlistKernel$isComputedXtX, testlistKernel$XtXgap, testlistKernel$ImSX, testlistKernel$isComputedImSX) expected <- as.numeric(Matrix::solve(Matrix::t(testImSX[, testcps, drop = FALSE]) %*% testImSX[, testcps, drop = FALSE], as.numeric(Matrix::t(testImSX[, testcps, drop = FALSE]) %*% testImSy))) expect_equal(object, expected, tolerance = 1e-2) }) test_that(".pelt works if there are no cps", { testy <- rnorm(100) object <- .pelt(testy, sd = 1) testpelt <- changepoint::cpt.mean(data = testy[6:95] / 1, penalty = "Manual", pen.value = 2 * log(90), method = "PELT") expect_identical(object$cps, testpelt@cpts[-length(testpelt@cpts)] + 6L) expect_equal(unique(object$est), testpelt@param.est$mean, tolerance = 1e-14) for (i in object$cps) { expect_gt(abs(object$est[i] - object$est[i - 1]), 0) } }) test_that(".pelt works if there is one cps", { testy <- rnorm(100) testy[1:50] <- testy[1:50] + 5 object <- .pelt(testy, sd = 1) testpelt <- changepoint::cpt.mean(data = testy[6:95] / 1, penalty = "Manual", pen.value = 2 * log(90), method = "PELT") expect_identical(object$cps, testpelt@cpts[-length(testpelt@cpts)] + 6L) expect_equal(unique(object$est), testpelt@param.est$mean, tolerance = 1e-14) for (i in object$cps) { expect_gt(abs(object$est[i] - object$est[i - 1]), 0) } }) test_that("sd in .pelt works if there are no cps", { testy <- rnorm(100) object <- .pelt(testy, sd = 0.1) testpelt <- changepoint::cpt.mean(data = testy[6:95] / 0.1, penalty = "Manual", pen.value = 2 * log(90), method = "PELT") expect_identical(object$cps, testpelt@cpts[-length(testpelt@cpts)] + 6L) expect_equal(unique(object$est), testpelt@param.est$mean * 0.1, tolerance = 1e-14) for (i in object$cps) { expect_gt(abs(object$est[i] - object$est[i - 1]), 0) } }) test_that("sd in .pelt works if there is one cps", { testy <- rnorm(100) testy[1:50] <- testy[1:50] + 5 object <- .pelt(testy, sd = 100) testpelt <- changepoint::cpt.mean(data = testy[6:95] / 100, penalty = "Manual", pen.value = 2 * log(90), method = "PELT") expect_identical(object$cps, testpelt@cpts[-length(testpelt@cpts)] + 6L) expect_equal(unique(object$est), testpelt@param.est$mean * 100, tolerance = 1e-14) for (i in object$cps) { expect_gt(abs(object$est[i] - object$est[i - 1]), 0) } })