Nothing
tests/testthat/testBasicFunctions.R
In PCpluS: Piecewise Constant Plus Smooth Regression
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)
}
})
Try the PCpluS package in your browser
Any scripts or data that you put into this service are public.
PCpluS documentation built on April 4, 2025, 2:14 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
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)
}
})
Try the PCpluS package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.