R Under development (unstable) (2026-03-22 r89674 ucrt) -- "Unsuffered Consequences" Copyright (C) 2026 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > > > require(stepR) Loading required package: stepR Successfully loaded stepR package version 2.1-11. Several new functions are added in version 2.0-0. Some older functions are deprecated (still working) and may be defunct in a later version. Please read the documentation for more details. > all.eq <- function(x, y, eps = 1e-5) TRUE #all(abs(x - y) < eps) > > # check Gauss var bounds > # y <- c(-2:2, 4) > y <- c(0, 2:5, 200, 7) > quant <- 2 > # without penalty > bs <- bounds.MRC(y, q = quant, family = "gaussvar", eps = 1e-5) > b <- bs$bounds > b li ri lower upper 1 1 1 0.0000000 0.000000e+00 2 1 2 0.4439274 3.811767e+01 3 1 4 2.3043721 4.571412e+01 4 2 2 0.5766308 5.896390e+02 5 2 3 1.4427639 1.238824e+02 6 2 5 4.2908998 8.512284e+01 7 3 3 1.2974193 1.326688e+03 8 3 4 2.7745461 2.382355e+02 9 3 6 3182.4173844 6.313277e+04 10 4 4 2.3065233 2.358556e+03 11 4 5 4.5502555 3.907062e+02 12 4 7 3185.5958287 6.319582e+04 13 5 5 3.6039426 3.685244e+03 14 5 6 4442.0482273 3.814150e+05 15 6 6 5766.3081873 5.896390e+06 16 6 7 4444.7117915 3.816437e+05 17 7 7 7.0637275 7.223078e+03 > meanY2 <- sapply(1:nrow(b), function(i) mean(y[b$li[i]:b$ri[i]]^2)) > len <- b$ri - b$li + 1 > # len / 2 * ( -1 - log(meanY2 / b$lower) + meanY2 / b$lower ) - quant > # len / 2 * ( -1 - log(meanY2 / b$upper) + meanY2 / b$upper ) - quant > stopifnot(all(abs(ifelse(meanY2 == 0, b$lower, len / 2 * ( -1 - log(meanY2 / b$lower) + meanY2 / b$lower ) - quant)) < 1e-4 )) > stopifnot(all(abs(ifelse(meanY2 == 0, b$upper, len / 2 * ( -1 - log(meanY2 / b$upper) + meanY2 / b$upper ) - quant)) < 1e-4 )) > # check BoundGaussVar > cand <- stepcand(y, family = "gaussvar") > as.data.frame(cand) leftEnd rightEnd value leftIndex rightIndex cumSumSq cumSumWe number improve 1 1 1 0 1 1 0 1 0 NA 2 2 2 4 2 2 4 2 0 NA 3 3 3 9 3 3 13 3 0 NA 4 4 4 16 4 4 29 4 0 NA 5 5 5 25 5 5 54 5 0 NA 6 6 6 40000 6 6 40054 6 0 NA 7 7 7 49 7 7 40103 7 0 NA > bounded <- stepbound(cand, bs) > as.data.frame(bounded) leftEnd rightEnd value leftIndex rightIndex rightIndexLeftBound 1 1 1 0.000 1 1 1 2 2 5 13.500 2 5 5 3 6 7 7223.078 6 7 7 rightIndexRightBound rightEndLeftBound rightEndRightBound cumSumWe cumSumSq 1 1 1 1 1 0 2 5 5 5 5 54 3 7 7 7 7 40103 > # twice negative log-likelihood > stopifnot(abs(attr(bounded, "cost") + sum(y != 0) * log(2 * pi) + 2 * sum(ifelse(fitted(bounded) == 0, ifelse(y ==0, 0, Inf), dnorm(y, 0, sqrt(fitted(bounded)), log = TRUE)))) < 1e-4 ) > > # with log(length) penalty > bs <- bounds.MRC(y, q = quant, family = "gaussvar", penalty = "len", eps = 1e-5) > b <- bs$bounds > b li ri lower upper 1 1 1 0.0000000 0.000000e+00 2 1 2 0.3303939 1.385843e+02 3 1 4 2.0448560 6.318500e+01 4 2 2 0.3534120 2.908498e+04 5 2 3 1.0737803 4.503988e+02 6 2 5 3.8076628 1.176548e+02 7 3 3 0.7951770 6.544120e+04 8 3 4 2.0649622 8.661516e+02 9 3 6 2824.0165815 8.726066e+04 10 4 4 1.4136479 1.163399e+05 11 4 5 3.3865379 1.420489e+03 12 4 7 2826.8370724 8.734782e+04 13 5 5 2.2088249 1.817811e+05 14 5 6 3306.0044044 1.386709e+06 15 6 6 3534.1197894 2.908498e+08 16 6 7 3307.9867681 1.387540e+06 17 7 7 4.3292967 3.562910e+05 > meanY2 <- sapply(1:nrow(b), function(i) mean(y[b$li[i]:b$ri[i]]^2)) > len <- b$ri - b$li + 1 > # len / 2 * ( -1 - log(meanY2 / b$lower) + meanY2 / b$lower ) - quant > # len / 2 * ( -1 - log(meanY2 / b$upper) + meanY2 / b$upper ) - quant > stopifnot(all(abs(ifelse(meanY2 == 0, b$lower, len / 2 * ( -1 - log(meanY2 / b$lower) + meanY2 / b$lower ) - quant + log(len / length(y)) )) < 1e-4 )) > stopifnot(all(abs(ifelse(meanY2 == 0, b$upper, len / 2 * ( -1 - log(meanY2 / b$upper) + meanY2 / b$upper ) - quant + log(len / length(y)) )) < 1e-4 )) > > # with sqrt penalty > bs <- bounds.MRC(y, q = quant, family = "gaussvar", penalty = "sqrt", eps = 1e-15) > b <- bs$bounds > b li ri lower upper 1 1 1 0.0000000 0.000000e+00 2 1 2 0.1669260 2.666428e+04 3 1 4 1.1323427 6.760178e+02 4 2 2 0.1682832 3.539826e+09 5 2 3 0.5425094 8.665890e+04 6 2 5 2.1085002 1.258792e+03 7 3 3 0.3786373 7.964609e+09 8 3 4 1.0432872 1.666517e+05 9 3 6 1563.8043081 9.336039e+05 10 4 4 0.6731329 1.415931e+10 11 4 5 1.7109911 2.733088e+05 12 4 7 1565.3661601 9.345364e+05 13 5 5 1.0517702 2.212391e+10 14 5 6 1670.3028831 2.668094e+08 15 6 6 1682.8323540 3.539826e+13 16 6 7 1671.3044389 2.669694e+08 17 7 7 2.0614696 4.336287e+10 > stopifnot(all(abs(ifelse(meanY2 == 0, b$lower, sqrt(2) * sqrt( len / 2 * ( -1 - log(meanY2 / b$lower) + meanY2 / b$lower ) ) - quant - sqrt(2*(1+log(length(y)/len))) )) < 1e-4 )) > stopifnot(all(abs(ifelse(meanY2 == 0, b$upper, sqrt(2) * sqrt(len / 2 * ( -1 - log(meanY2 / b$upper) + meanY2 / b$upper )) - quant - sqrt(2*(1+log(length(y)/len))) )) < 1e-4 )) > > # check BoundGaussVar > cand <- stepcand(y, family = "gaussvar") > as.data.frame(cand) leftEnd rightEnd value leftIndex rightIndex cumSumSq cumSumWe number improve 1 1 1 0 1 1 0 1 0 NA 2 2 2 4 2 2 4 2 0 NA 3 3 3 9 3 3 13 3 0 NA 4 4 4 16 4 4 29 4 0 NA 5 5 5 25 5 5 54 5 0 NA 6 6 6 40000 6 6 40054 6 0 NA 7 7 7 49 7 7 40103 7 0 NA > bounded <- stepbound(cand, bs) > as.data.frame(bounded) leftEnd rightEnd value leftIndex rightIndex rightIndexLeftBound 1 1 1 0.0 1 1 1 2 2 5 13.5 2 5 2 3 6 7 20024.5 6 7 7 rightIndexRightBound rightEndLeftBound rightEndRightBound cumSumWe cumSumSq 1 1 1 1 1 0 2 5 2 5 5 54 3 7 7 7 7 40103 > # twice negative log-likelihood > stopifnot(abs(attr(bounded, "cost") + sum(y != 0) * log(2 * pi) + 2 * sum(ifelse(fitted(bounded) == 0, ifelse(y ==0, 0, Inf), dnorm(y, 0, sqrt(fitted(bounded)), log = TRUE)))) < 1e-4 ) > > # check Binomial bounds > # y <- c(0, 0, 1, 2, 2) > # size <- 2 > y <- c(0, 0, 1, 0, 1, 1, 1, 0) > size <- 1 > quant <- 2 > # without penalty > b <- bounds.MRC(y, q = quant, family = "binom", param = size, eps = 1e-5)$bounds > b li ri lower upper 1 1 1 0.00000000 0.8646647 2 1 2 0.00000000 0.6321206 3 1 4 0.01493266 0.7306796 4 1 8 0.18636433 0.8136357 5 2 2 0.00000000 0.8646647 6 2 3 0.03506325 0.9649367 7 2 5 0.10246995 0.8975300 8 3 3 0.13533528 1.0000000 9 3 4 0.03506325 0.9649367 10 3 6 0.26932042 0.9850673 11 4 4 0.00000000 0.8646647 12 4 5 0.03506325 0.9649367 13 4 7 0.26932042 0.9850673 14 5 5 0.13533528 1.0000000 15 5 6 0.36787944 1.0000000 16 5 8 0.26932042 0.9850673 17 6 6 0.13533528 1.0000000 18 6 7 0.36787944 1.0000000 19 7 7 0.13533528 1.0000000 20 7 8 0.03506325 0.9649367 21 8 8 0.00000000 0.8646647 > S <- sapply(1:nrow(b), function(i) sum(y[b$li[i]:b$ri[i]])) > len <- b$ri - b$li + 1 > sizelen <- size * len > NS <- sizelen - S > stopifnot(all(ifelse(S == 0, b$lower, ifelse(NS == 0, -sizelen * log(b$lower), S * log(S / sizelen / b$lower) + NS * log(NS / sizelen / (1 - b$lower))) - quant) < 1e-4)) > stopifnot(all(ifelse(NS == 0, b$upper - 1, ifelse(S == 0, -sizelen * log(1 - b$upper), S * log(S / sizelen / b$upper) + NS * log(NS / sizelen / (1 - b$upper))) - quant) < 1e-4)) > # with len-penalty > b <- bounds.MRC(y, q = quant, family = "binom", param = size, penalty = "len", eps = 1e-5)$bounds > b li ri lower upper 1 1 1 0.000000000 0.9830831 2 1 2 0.000000000 0.8160603 3 1 4 0.007295325 0.7918968 4 1 8 0.186364327 0.8136357 5 2 2 0.000000000 0.9830831 6 2 3 0.008531237 0.9914688 7 2 5 0.069921533 0.9300785 8 3 3 0.016916910 1.0000000 9 3 4 0.008531237 0.9914688 10 3 6 0.208103196 0.9927047 11 4 4 0.000000000 0.9830831 12 4 5 0.008531237 0.9914688 13 4 7 0.208103196 0.9927047 14 5 5 0.016916910 1.0000000 15 5 6 0.183939721 1.0000000 16 5 8 0.208103196 0.9927047 17 6 6 0.016916910 1.0000000 18 6 7 0.183939721 1.0000000 19 7 7 0.016916910 1.0000000 20 7 8 0.008531237 0.9914688 21 8 8 0.000000000 0.9830831 > S <- sapply(1:nrow(b), function(i) sum(y[b$li[i]:b$ri[i]])) > len <- b$ri - b$li + 1 > sizelen <- size * len > NS <- sizelen - S > stopifnot(all(ifelse(S == 0, b$lower,abs( ifelse(NS == 0, -sizelen * log(b$lower), S * log(S / sizelen / b$lower) + NS * log(NS / sizelen / (1 - b$lower))) - quant + log(len / length(y)))) < 1e-4)) > stopifnot(all(ifelse(NS == 0, b$upper - 1, ifelse(S == 0, -sizelen * log(1 - b$upper), S * log(S / sizelen / b$upper) + NS * log(NS / sizelen / (1 - b$upper))) - quant + log(len / length(y))) < 1e-4)) > # with var-penalty > b <- bounds.MRC(y, q = quant, family = "binom", param = size, penalty = "var", eps = 1e-5)$bounds > b li ri lower upper 1 1 2 0.0000000 0.8807971 2 1 4 0.0000000 0.8310406 3 1 8 0.1512225 0.8487775 4 2 3 0.0000000 1.0000000 5 2 5 0.0172132 0.9827868 6 3 4 0.0000000 1.0000000 7 3 6 0.1689594 1.0000000 8 4 5 0.0000000 1.0000000 9 4 7 0.1689594 1.0000000 10 5 6 0.1192029 1.0000000 11 5 8 0.1689594 1.0000000 12 6 7 0.1192029 1.0000000 13 7 8 0.0000000 1.0000000 > S <- sapply(1:nrow(b), function(i) sum(y[b$li[i]:b$ri[i]])) > len <- b$ri - b$li + 1 > sizelen <- size * len > NS <- sizelen - S > totvar <- ( sum(y[-length(y)] * (size - y[-1])) + sum(y[-1] * (size - y[-length(y)])) ) / 2 / size > totvar [1] 2 > stopifnot(all(ifelse(S <= 1, b$lower, S * log(S / sizelen) + ifelse(NS == 0, 0, NS * log(NS / sizelen)) - quant + log(sizelen) - log(totvar) - (S - 1) * log(b$lower) - (NS - 1) * log(1 - b$lower)) < 1e-4)) > stopifnot(all(ifelse(NS <= 1, b$upper - 1, ifelse(S == 0, 0, S * log(S / sizelen)) + NS * log(NS / sizelen) - quant + log(sizelen) - log(totvar) - (S - 1) * log(b$upper) - (NS - 1) * log(1 - b$upper)) < 1e-4)) > #with sqrt penalty > b <- bounds.MRC(y, q = quant, family = "binom", param = size, penalty = "sqrt", eps = 1e-5)$bounds > b li ri lower upper 1 1 1 0.000000e+00 0.9999565 2 1 2 0.000000e+00 0.9874467 3 1 4 6.621083e-05 0.9589785 4 1 8 6.208139e-02 0.9379186 5 2 2 0.000000e+00 0.9999565 6 2 3 3.939781e-05 0.9999606 7 2 5 6.302974e-03 0.9936970 8 3 3 4.349516e-05 1.0000000 9 3 4 3.939781e-05 0.9999606 10 3 6 4.102148e-02 0.9999338 11 4 4 0.000000e+00 0.9999565 12 4 5 3.939781e-05 0.9999606 13 4 7 4.102148e-02 0.9999338 14 5 5 4.349516e-05 1.0000000 15 5 6 1.255329e-02 1.0000000 16 5 8 4.102148e-02 0.9999338 17 6 6 4.349516e-05 1.0000000 18 6 7 1.255329e-02 1.0000000 19 7 7 4.349516e-05 1.0000000 20 7 8 3.939781e-05 0.9999606 21 8 8 0.000000e+00 0.9999565 > S <- sapply(1:nrow(b), function(i) sum(y[b$li[i]:b$ri[i]])) > len <- b$ri - b$li + 1 > sizelen <- size * len > NS <- sizelen - S > stopifnot(all(abs(ifelse(S == 0, b$lower, ifelse(NS == 0, sqrt(2)*sqrt(-sizelen * log(b$lower)), sqrt(2)*sqrt(S * log(S / sizelen / b$lower) + NS * log(NS / sizelen / (1 - b$lower)))) - quant - sqrt(2*(1+log(length(y)/len))) )) < 1e-4)) > stopifnot(all(ifelse(NS == 0, b$upper - 1, ifelse(S == 0,sqrt(2)*sqrt(-sizelen * log(1 - b$upper)),sqrt(2)*sqrt(S * log(S / sizelen / b$upper) + NS * log(NS / sizelen / (1 - b$upper)))) - quant - sqrt(2*(1+log(length(y)/len)))) < 1e-4)) > > # check Poisson bounds > y <- c(0,0,1,1) > quant <- 2 > # without penalty > b <- bounds.MRC(y, q = quant, family = "poisson", eps = 1e-5)$bounds > b li ri lower upper 1 1 1 0.00000000 2.000000 2 1 2 0.00000000 1.000000 3 1 4 0.07929714 1.573094 4 2 2 0.00000000 2.000000 5 2 3 0.02623455 2.252621 6 3 3 0.05246910 4.505241 7 3 4 0.15859428 3.146189 8 4 4 0.05246910 4.505241 > S <- sapply(1:nrow(b), function(i) sum(y[b$li[i]:b$ri[i]])) > len <- b$ri - b$li + 1 > stopifnot(all(ifelse(S == 0, b$lower, b$lower * ( S / b$lower * log(S / b$lower / len) - S / b$lower + len ) - quant) < 1e-4)) > stopifnot(all(ifelse(S == 0, b$upper * len, b$upper * ( S / b$upper * log(S / b$upper / len) - S / b$upper + len )) - quant < 1e-4)) > # S = 0 > bu0 <- b$upper[1] > stopifnot(abs(bu0 - quant) < 1e-5) > stopifnot(b$lower[1] == 0) > bu00 <- b$upper[2] > stopifnot(abs(2 * bu00 - quant) < 1e-5) > stopifnot(b$lower[2] == 0) > # S = 2 > bu11 <- b$upper[7] > stopifnot(abs(2 * log(2 / 2 / bu11) - 2 + 2 * bu11 - quant) < 1e-5) > bl11 <- b$lower[7] > stopifnot(abs(2 * log(2 / 2 / bl11) - 2 + 2 * bl11 - quant) < 1e-5) > # with len-penalty > b <- bounds.MRC(y, q = quant, family = "poisson", penalty = "len", eps = 1e-5)$bounds > b li ri lower upper 1 1 1 0.00000000 3.386294 2 1 2 0.00000000 1.346574 3 1 4 0.07929714 1.573097 4 2 2 0.00000000 3.386294 5 2 3 0.01276872 2.687442 6 3 3 0.01260465 6.212926 7 3 4 0.10644479 3.638011 8 4 4 0.01260465 6.212926 > S <- sapply(1:nrow(b), function(i) sum(y[b$li[i]:b$ri[i]])) > len <- b$ri - b$li + 1 > stopifnot(all(ifelse(S == 0, b$lower, b$lower * ( S / b$lower * log(S / b$lower / len) - S / b$lower + len ) - quant + log(len / length(y))) < 1e-4)) > stopifnot(all(ifelse(S == 0, b$upper * len, b$upper * ( S / b$upper * log(S / b$upper / len) - S / b$upper + len )) - quant + log(len / length(y)) < 1e-4)) > # with sqrt penalty > b <- bounds.MRC(y, q = quant, family = "poisson", penalty = "sqrt", eps = 1e-5)$bounds > b li ri lower upper 1 1 1 0.000000e+00 8.755545 2 1 2 0.000000e+00 3.686762 3 1 4 1.018338e-02 2.822490 4 2 2 0.000000e+00 8.755545 5 2 3 1.154768e-04 5.374135 6 3 3 5.797565e-05 12.262055 7 3 4 9.302612e-03 6.569146 8 4 4 5.797565e-05 12.262055 > S <- sapply(1:nrow(b), function(i) sum(y[b$li[i]:b$ri[i]])) > len <- b$ri - b$li + 1 > stopifnot(all(ifelse(S == 0,sqrt(2)*sqrt(b$lower * len), sqrt(2) * sqrt(b$lower * ( S / b$lower * log(S / b$lower / len) - S / b$lower + len ))) - quant - sqrt(2*(1+log(length(y)/len))) < 1e-4)) > stopifnot(all(ifelse(S == 0,sqrt(2)*sqrt(b$upper * len), sqrt(2) * sqrt(b$upper * ( S / b$upper * log(S / b$upper / len) - S / b$upper + len ))) - quant - sqrt(2*(1+log(length(y)/len))) < 1e-4)) > > # with var-penalty > b <- bounds.MRC(y, q = quant, family = "poisson", penalty = "var", eps = 1e-5)$bounds > b li ri lower upper 1 1 1 0.0000000 2.000000 2 1 2 0.0000000 1.000000 3 1 4 0.1015939 1.223771 4 2 2 0.0000000 2.000000 5 2 3 0.0000000 1.846574 6 3 3 0.0000000 3.693147 7 3 4 0.2031879 2.447542 8 4 4 0.0000000 3.693147 > S <- sapply(1:nrow(b), function(i) sum(y[b$li[i]:b$ri[i]])) > len <- b$ri - b$li + 1 > ifelse(S == 0, b$lower, b$lower * ( S / b$lower * log(S / b$lower / len) - S / b$lower + len ) - quant + log(b$lower * len / sum(y))) [1] 0 0 -2 0 NaN NaN -2 NaN > stopifnot(all(ifelse(S <= 1, b$lower, b$lower * ( S / b$lower * log(S / b$lower / len) - S / b$lower + len ) - quant + log(b$lower * len / sum(y))) < 1e-4)) > stopifnot(all(ifelse(S == 0, b$upper * len, b$upper * ( S / b$upper * log(S / b$upper / len) - S / b$upper + len )) - quant + log(b$upper * len / sum(y)) < 1e-4)) > > # S = 0 > bu0 <- b$upper[1] > stopifnot(abs(bu0 + log(bu0) - quant - log(sum(y))) < 1e-5) > stopifnot(b$lower[1] == 0) > bu00 <- b$upper[2] > stopifnot(abs(2 * bu00 + log(2 * bu00) - quant - log(sum(y))) < 1e-5) > stopifnot(b$lower[2] == 0) > # S = 1 > bu1 <- b$upper[6] > stopifnot(abs(bu1 - 1 - quant - log(sum(y))) < 1e-5) > stopifnot(b$lower[6] == 0) > bu01 <- b$upper[5] > stopifnot(abs(2 * bu01 - 1 - quant - log(sum(y))) < 1e-5) > stopifnot(b$lower[5] == 0) > > > # check BoundBinom > y <- 1:4 > size <- 4 > cand <- stepcand(y, family = "binomial", param = size) > bounds <- as.data.frame(rbind( + c(1, 1, 0, 1), c(1, 2, 1, 0), c(3, 3, 2, 4), c(3, 4, 3, 4), c(4, 4, 4, 4) + )) > names(bounds) <- c("li", "ri", "lower", "upper") > bounds <- bounds[order(bounds$li, bounds$ri),] > start <- cumsum(sapply(tapply(bounds$li, ordered(bounds$li, levels = 1:nrow(cand)), identity), length)) > start <- c(0, start[-length(start)]) # C-style > start[is.na(tapply(bounds$li, ordered(bounds$li, levels = 1:nrow(cand)), length))] <- NA > with(bounds, cbind(bounds, Cli = li - 1, Cri = ri - 1, Crows = 0:(nrow(bounds)-1))) li ri lower upper Cli Cri Crows 1 1 1 0 1 0 0 0 2 1 2 1 0 0 1 1 3 3 3 2 4 2 2 2 4 3 4 3 4 2 3 3 5 4 4 4 4 3 3 4 > cbind(as.data.frame(cand[,2:3]), start = start) rightEnd value start 1 1 0.25 0 2 2 0.50 NA 3 3 0.75 2 4 4 1.00 4 > # normalise bounds > bbounds <- bounds > bbounds$lower <- bbounds$lower / size > bbounds$upper <- bbounds$upper / size > bounded <- stepbound(cand, list(bounds = bbounds, start = start, feasible = TRUE)) > as.data.frame(bounded) leftEnd rightEnd value leftIndex rightIndex rightIndexLeftBound 1 1 1 0.250 1 1 1 2 2 3 0.625 2 3 3 3 4 4 1.000 4 4 4 rightIndexRightBound rightEndLeftBound rightEndRightBound cumSum cumSumWe 1 1 1 1 1 1 2 3 3 3 6 3 3 4 4 4 10 4 > stopifnot(all.equal(bounded$rightEnd, c(1, 3, 4))) > stopifnot(all.eq(bounded$value, c(1, 2.5, 4) / size)) > # attributes(bounded) > stopifnot(abs(attr(bounded, "cost") - sum(lchoose(size, y)) +sum(dbinom(y, size, fitted(bounded) / size, log = TRUE)))<0.001) > > # check BoundPoisson > cand <- stepcand(y, family = "poisson") > bounded <- stepbound(cand, list(bounds = bounds, start = start, feasible = TRUE)) > as.data.frame(bounded) leftEnd rightEnd value leftIndex rightIndex rightIndexLeftBound 1 1 1 1 1 1 1 2 2 4 4 2 4 4 rightIndexRightBound rightEndLeftBound rightEndRightBound cumSum cumSumWe 1 1 1 1 1 1 2 4 4 4 10 4 > stopifnot(all.equal(bounded$rightEnd, c(1, 4))) > stopifnot(all.eq(bounded$value, c(1, 4))) > # attributes(bounded) > attr(bounded, "cost") [1] 0.5233507 > stopifnot(abs(attr(bounded, "cost") + sum(lfactorial(y)) +sum(dpois(y, fitted(bounded), log = TRUE)))<0.001) > > # check BoundGauss > cand <- stepcand(y, family = "gauss") > # # call with C-style indices > # bounded <- with(bounds, .Call('boundedGauss', cand$cumSum, cand$cumSumSq, cand$cumSumWe, as.integer(start), as.integer(ri - 1), as.numeric(lower), as.numeric(upper))) > bounded <- stepbound(cand, list(bounds = bounds, start = start, feasible = TRUE)) > as.data.frame(bounded) leftEnd rightEnd value leftIndex rightIndex rightIndexLeftBound 1 1 1 1 1 1 1 2 2 4 4 2 4 4 rightIndexRightBound rightEndLeftBound rightEndRightBound cumSum cumSumWe 1 1 1 1 1 1 2 4 4 4 10 4 cumSumSq 1 1 2 30 > stopifnot(all.equal(bounded$rightEnd, c(1, 4))) > stopifnot(all.eq(bounded$value, c(1, 4))) > # attributes(bounded) > attr(bounded, "cost") [1] 5 > stopifnot(attr(bounded, "cost") == 4 + 1) > y <- (-4):4 > MRCoeff(y, lengths = c(1,4,9), signed = TRUE) [,1] [,2] [,3] [1,] -4 -5 4.849887e-17 [2,] -3 -3 NA [3,] -2 -1 NA [4,] -1 1 NA [5,] 0 3 NA [6,] 1 5 NA [7,] 2 NA NA [8,] 3 NA NA [9,] 4 NA NA > sd <- 0.4 > MRC.quant(1 - 0.05, 9, 1e2) * sd 95% 1.730352 > b <- bounds(y, r = 1e2, param = sd, lengths = c(1,4,9)) > b $bounds li ri lower upper 1 1 1 -5.2940495 -2.7059505 2 1 4 -3.1470248 -1.8529752 3 1 9 -0.4313498 0.4313498 4 2 2 -4.2940495 -1.7059505 5 2 5 -2.1470248 -0.8529752 6 3 3 -3.2940495 -0.7059505 7 3 6 -1.1470248 0.1470248 8 4 4 -2.2940495 0.2940495 9 4 7 -0.1470248 1.1470248 10 5 5 -1.2940495 1.2940495 11 5 8 0.8529752 2.1470248 12 6 6 -0.2940495 2.2940495 13 6 9 1.8529752 3.1470248 14 7 7 0.7059505 3.2940495 15 8 8 1.7059505 4.2940495 16 9 9 2.7059505 5.2940495 $start [1] 0 3 5 7 9 11 13 14 15 $feasible [1] TRUE attr(,"class") [1] "bounds" "list" > sb <- stepbound(y, b) > sb Fitted step function of family gauss containing 3 blocks domain: ( 0 , 9 ] range: [ -3 , 3 ] cost: 6 > as.data.frame(sb) leftEnd rightEnd value leftIndex rightIndex rightIndexLeftBound 1 1 3 -3 1 3 1 2 4 6 0 4 6 4 3 7 9 3 7 9 9 rightIndexRightBound rightEndLeftBound rightEndRightBound cumSum cumSumWe 1 3 1 3 -9 3 2 6 4 6 -9 6 3 9 9 9 0 9 cumSumSq 1 29 2 31 3 60 > stopifnot(nrow(sb) == 3) > stopifnot(all.equal(sb$rightEnd, c(3, 6, 9))) > bs <- bounds(y, r = 1e2, subset = c(1,4:5,9), param = sd, lengths = c(1,4,9)) > bs $bounds li ri lower upper 1 1 1 -5.16329679 -2.83670321 2 1 2 -3.08164839 -1.91835161 3 1 4 -0.38776560 0.38776560 4 2 2 -4.16329679 -1.83670321 6 2 2 -3.16329679 -0.83670321 8 2 2 -2.16329679 0.16329679 5 2 3 -2.08164839 -0.91835161 7 2 4 -1.08164839 0.08164839 9 2 4 -0.08164839 1.08164839 10 3 3 -1.16329679 1.16329679 11 3 4 0.91835161 2.08164839 12 4 4 -0.16329679 2.16329679 13 4 4 1.91835161 3.08164839 14 4 4 0.83670321 3.16329679 15 4 4 1.83670321 4.16329679 16 4 4 2.83670321 5.16329679 $start [1] 0 3 9 11 $feasible [1] FALSE attr(,"class") [1] "bounds" "list" > stopifnot(!bs$feasible) > sub <- c(2,4:7,9) > bs <- b[sub] > bs $bounds li ri lower upper 1 1 1 -5.2940495 -2.7059505 4 1 1 -4.2940495 -1.7059505 2 1 2 -3.1470248 -1.8529752 5 1 3 -2.1470248 -0.8529752 3 1 6 -0.4313498 0.4313498 6 2 2 -3.2940495 -0.7059505 8 2 2 -2.2940495 0.2940495 7 2 4 -1.1470248 0.1470248 9 2 5 -0.1470248 1.1470248 10 3 3 -1.2940495 1.2940495 11 3 6 0.8529752 2.1470248 12 4 4 -0.2940495 2.2940495 13 4 6 1.8529752 3.1470248 14 5 5 0.7059505 3.2940495 15 6 6 1.7059505 4.2940495 16 6 6 2.7059505 5.2940495 $start [1] 0 5 9 11 13 14 $feasible [1] TRUE > cand <- stepcand(y) > as.data.frame(cand[sub,]) leftEnd rightEnd value leftIndex rightIndex cumSum cumSumSq cumSumWe number 2 1 2 -3 1 2 -7 25 2 0 4 3 4 -1 3 4 -10 30 4 0 5 5 5 0 5 5 -10 30 5 0 6 6 6 1 6 6 -9 31 6 0 7 7 7 2 7 7 -7 35 7 0 9 8 9 4 8 9 0 60 9 0 improve 2 NA 4 NA 5 NA 6 NA 7 NA 9 NA > sb <- stepbound(cand[sub,], bs) > sb Fitted step function of family gauss containing 4 blocks domain: ( 0 , 9 ] range: [ -3.5 , 3.5 ] cost: 3.5 > as.data.frame(sb) leftEnd rightEnd value leftIndex rightIndex rightIndexLeftBound 1 1 2 -3.5 1 2 2 2 3 5 -1.0 3 5 4 3 6 7 1.5 6 7 6 4 8 9 3.5 8 9 9 rightIndexRightBound rightEndLeftBound rightEndRightBound cumSum cumSumWe 1 2 2 2 -7 2 2 5 4 5 -10 5 3 7 6 7 -7 7 4 9 9 9 0 9 cumSumSq 1 25 2 30 3 35 4 60 > stopifnot(nrow(sb) == 4) > stopifnot(all.equal(sb$rightEnd, c(2, 5, 7, 9))) > > # check whether candidates and steppath return correct number of results > example(stepcand) stpcnd> # simulate 5 blocks (4 jumps) within a total of 100 data points stpcnd> b <- c(sort(sample(1:99, 4)), 100) stpcnd> f <- rep(rnorm(5, 0, 4), c(b[1], diff(b))) stpcnd> rbind(b = b, f = unique(f), lambda = exp(unique(f) / 10) * 20) [,1] [,2] [,3] [,4] [,5] b 19.000000 34.00000 43.000000 94.0000000 100.0000000 f 1.944007 -4.66451 6.031567 0.8125491 0.4867767 lambda 24.291656 12.54449 36.557595 21.6929470 20.9976378 stpcnd> # add gaussian noise stpcnd> x <- f + rnorm(100) stpcnd> # find 10 candidate jumps stpcnd> stepcand(x, max.cand = 10) Fitted step function of family gauss containing 10 blocks domain: ( 0 , 100 ] range: [ -5.824665 , 8.142148 ] cost: 75.35003 stpcnd> # for poisson observations stpcnd> y <- rpois(100, exp(f / 10) * 20) stpcnd> # find 10 candidate jumps stpcnd> stepcand(y, max.cand = 10, family = "poisson") Fitted step function of family poisson containing 10 blocks domain: ( 0 , 100 ] range: [ 11.6 , 37.71429 ] cost: -4411.708 stpcnd> # for binomial observations stpcnd> size <- 10 stpcnd> z <- rbinom(100, size, pnorm(f / 10)) stpcnd> # find 10 candidate jumps stpcnd> stepcand(z, max.cand = 10, family = "binomial", param = size) Fitted step function of family binomial containing 10 blocks domain: ( 0 , 100 ] range: [ 0.2615385 , 0.7666667 ] cost: 656.4924 parameter: [1] 10 > cand <- stepcand(x, max.cand = 100) > stopifnot(nrow(cand) == 100) > print(cand) Fitted step function of family gauss containing 100 blocks domain: ( 0 , 100 ] range: [ -6.798015 , 8.142148 ] cost: 0 > stopifnot(attr(cand, "cost") == 0) > system.time(stopifnot(length(steppath(cand)) == 100)) user system elapsed 0 0 0 > system.time(stopifnot(length(steppath(cand, max.blocks = 10)) == 10)) user system elapsed 0 0 0 > stopifnot(nrow(stepcand(x, max.cand = 10)) == 10) > stopifnot(nrow(print(stepcand(x, max.cand = 1))) == 1) Fitted step function of family gauss containing 1 blocks domain: ( 0 , 100 ] range: [ 0.648199 , 0.648199 ] cost: 846.1757 > stopifnot(nrow(print(stepcand(x[1], max.cand = 1))) == 1) Fitted step function of family gauss containing 1 blocks domain: ( NA , 1 ] range: [ 2.608031 , 2.608031 ] cost: 0 > pcand <- stepcand(y, max.cand = 100, family = "poisson") > stopifnot(nrow(pcand) == 100) > stopifnot(nrow(stepcand(y, max.cand = 10, family = "poisson")) == 10) > stopifnot(nrow(stepcand(y, max.cand = 1, family = "poisson")) == 1) > stopifnot(nrow(stepcand(y[1], max.cand = 1, family = "poisson")) == 1) > bcand <- stepcand(z, max.cand = 100, family = "binomial", param = size) > stopifnot(nrow(bcand) == 100) > stopifnot(nrow(stepcand(z, max.cand = 10, family = "binomial", param = size)) == 10) > stopifnot(nrow(stepcand(z, max.cand = 1, family = "binomial", param = size)) == 1) > stopifnot(nrow(stepcand(z[1], max.cand = 1, family = "binomial", param = size)) == 1) > stopifnot(nrow(stepcand(x, max.cand = 100)) == 100) > > # check forward selection > forward <- function(y, max.cand = length(y)) { + X <- as.data.frame(sapply(1:length(y), function(i) rep(c(1.0, 0), c(i, length(y) - i)))) + l <- lm(eval(parse(text = paste("y ~ 0 +", names(X)[ncol(X)]))), data = X) + ret <- data.frame(rightEnd = length(y), number = (1:max.cand) - 1, RSS = NA, improve = NA) + for(i in 2:max.cand) { + a <- add1(l, eval(parse(text = paste("~", paste(names(X), collapse = "+"))))) + m <- which.min(a$RSS) + v <- rownames(a)[m] + ret$rightEnd[i] <- as.integer(substring(v, 2)) + ret$RSS[i] <- a$RSS[m] + ret$improve[i] <- a$RSS[1] - a$RSS[m] + l <- eval(parse(text=paste("update(l, . ~ . +", v,")"))) + } + ret[order(ret$rightEnd),] + } > stopifnot(forward(x, 10)[,c("rightEnd", "number")] == stepcand(x, max.cand = 10)[,c("rightEnd", "number")]) > > # should select blocks of 4 > stopifnot(stepcand(1:16, max.cand = 4)$rightEnd == c(4, 8, 12, 16)) > # forward selection cuts in quarters, optimal solution in thirds > stopifnot(stepcand(1:12, max.cand = 3)$rightEnd == c(6, 9, 12)) > stopifnot(steppath(stepcand(1:12))[[3]]$rightEnd == c(4, 8, 12)) > # check RSS, likelihood of solution with one block > sp <- steppath(cand) > stopifnot(isTRUE(print(all.eq(sp$cost[1], sum( (x - mean(x))^2 ))))) [1] TRUE > stopifnot(isTRUE(print(all.eq(as.numeric(logLik(sp[[1]])), as.numeric(logLik( lm(x ~ 1) )))))) [1] TRUE > # check RSS of solution with 5 blocks > stopifnot(isTRUE(print(all.eq(sp$cost[5], + sum( apply(rbind(c(0, sp[[5]]$rightEnd[-5]) + 1, sp[[5]]$rightEnd), 2, function(i) sum( (x[i[1]:i[2]] - mean(x[i[1]:i[2]]))^2 ) ) ))))) [1] TRUE > # check likelihood if standard deviation is specified > attr(sp$cand, "param") <- .1 > stopifnot(isTRUE(print(all.eq(as.numeric(logLik(sp)[1]), as.numeric(sum(dnorm(x, mean(x), .1, log = TRUE))))))) [1] TRUE > # check Poisson likelihood of solution with 1 block > psp <- steppath(pcand) > psp.const <- sum( lfactorial(y) ) # data dependent constant > stopifnot(isTRUE(print(all.eq(-psp$cost[1] - psp.const, sum( dpois(y, mean(y), log = T) ))))) [1] TRUE > # check Poisson likelihood of solution with 5 blocks > stopifnot(isTRUE(print(all.eq(-psp$cost[5] - psp.const, + sum( apply(rbind(c(0, psp[[5]]$rightEnd[-5]) + 1, psp[[5]]$rightEnd), 2, function(i) sum( dpois(y[i[1]:i[2]], mean(y[i[1]:i[2]]), log = T) ) ) ))))) [1] TRUE > # check Binomial likelihood of solution with 1 block > bsp <- steppath(bcand) > bsp.const <- sum( lchoose(size, z) ) # data dependent constant > stopifnot(isTRUE(print(all.eq(-bsp$cost[1] + bsp.const, sum( dbinom(z, size, mean(z) / size, log = T) ))))) [1] TRUE > # check Binomial likelihood of solution with 5 blocks > stopifnot(isTRUE(print(all.eq(-bsp$cost[5] + bsp.const, + sum( apply(rbind(c(0, bsp[[5]]$rightEnd[-5]) + 1, bsp[[5]]$rightEnd), 2, function(i) sum( dbinom(z[i[1]:i[2]], size, mean(z[i[1]:i[2]]) / size, log = T) ) ) ))))) [1] TRUE > > # # check inhibition > # print(length(x)) > # icand <- stepcand(x, family = "gaussInhibitBoth", param = c(start = 3, middle = 4, end = 5)) > # stopifnot(min(icand$rightEnd) >= 3) > # stopifnot(min(diff(icand$rightEnd[-nrow(icand)])) >= 4) > # stopifnot(diff(icand$rightEnd[nrow(icand)-1:0]) >= 5) > # ipath <- steppath(x, family = "gaussInhibit", param = c(start = 3, middle = 4, end = 5)) > # print(ipath$path) > # print(ipath$cost) > # stopifnot(sapply(1:length(ipath), function(i) min(ipath[[i]]$rightEnd) >= 3)) > # print(sapply(3:length(ipath), function(i) length(diff(ipath[[i]]$rightEnd[-nrow(ipath[[i]])])))) > # stopifnot(sapply(3:length(ipath), function(i) min(diff(ipath[[i]]$rightEnd[-nrow(ipath[[i]])])) >= 4)) > # stopifnot(sapply(2:length(ipath), function(i) diff(ipath[[i]]$rightEnd[nrow(ipath[[i]])-1:0]) >= 5)) > > # check radius > blocks <- c(rep(0, 9), 1, 3, rep(1, 19)) > stopifnot(stepcand(blocks, max.cand = 3)$rightEnd == c(9, 11, 30)) > stopifnot(steppath(blocks)[[3]]$rightEnd == c(10, 11, 30)) > stopifnot(steppath(blocks, max.cand = 3, cand.radius = 1)[[3]]$rightEnd == c(10, 11, 30)) > > # check gaussKern with "exact" data > # simple test cases > N <- 300 > truth <- stepblock(0:3, rightEnd = c(0.2, 0.5, 0.6, 1) * N) > lapply(list( + dfilter("custom", diff(c(0, 0.1, 0.3, 0.4, 0.8, 1))), + dfilter("custom", dfilter(len = 9)$kern) + ), function(fkern) { + fkern$jump <- min(which(fkern$step >= 0.5)) - 1 + # print(fkern$step) + # print(fkern$jump) + signal.const <- rep(truth$value, diff(c(0, truth$rightEnd))) + signal <- convolve(c(rep(signal.const[1], length(fkern$kern) - fkern$jump - 1), signal.const, rep(signal.const[length(signal.const)], fkern$jump)), rev(fkern$kern), TRUE, "filter") + sc <- stepcand(signal, family = "gaussKern", param = fkern, max.cand = 5, cand.radius = length(fkern$kern)) + # print(sc$rightEnd) + # print(sc[,]) + sp <- steppath(sc) + # print(sp) + sp4 <- sp[[4]] + print(as.list(sp4)) + # print(sp4$rightEnd) + # compare exact values and estimates + if(!all.eq(truth$value, sp4$value)) print(rbind(truth$value, sp4$value)) + stopifnot(all.eq(truth$value, sp4$value)) + # compare exact signal and fit + if(!all.eq(signal, fitted(sp4))) print(rbind(signal, fitted(sp4))) + stopifnot(all.eq(signal, fitted(sp4))) + }) $leftEnd [1] 1 61 151 181 $rightEnd [1] 60 150 180 300 $value [1] 6.941258e-16 1.000000e+00 2.000000e+00 3.000000e+00 $leftIndex [1] 1 61 151 181 $rightIndex [1] 60 150 180 300 attr(,"x0") [1] 0 attr(,"cost") [1] 0.2522892 attr(,"family") [1] "gaussKern" attr(,"param") Custom filter length : 5 jump after: 3 step response start : 0.1 1 - step response end: 0 $leftEnd [1] 1 61 151 181 $rightEnd [1] 60 150 180 300 $value [1] 2.782523e-15 1.000000e+00 2.000000e+00 3.000000e+00 $leftIndex [1] 1 61 151 181 $rightIndex [1] 60 150 180 300 attr(,"x0") [1] 0 attr(,"cost") [1] 0.07805714 attr(,"family") [1] "gaussKern" attr(,"param") Custom filter length : 9 jump after: 3 step response start : 0 1 - step response end: 0 [[1]] NULL [[2]] NULL > # test refitting with blocks shorter than kernel length > bl <- c(6, 2, 10, 3, 2, 7) > bn <- length(bl) > bh <- rnorm(bn) > x <- rep(bh, bl) > k <- dfilter("custom", c(0, 0.3, 0.5, 0.2, 0)) > k Custom filter length : 5 jump after: 2 step response start : 0 1 - step response end: 0 > kl <- length(k$kern) > kj <- k$jump > y <- convolve(c(rep(bh[1], kj), x, rep(bh[bn], kl - kj - 1)), rev(k$kern), type = "filter") > rbind(x,y) [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] x 1.111608 1.111608 1.111608 1.111608 1.111608 1.1116079 -1.0866622 -1.086662 y 1.111608 1.111608 1.111608 1.111608 1.111608 0.4521269 -0.6470082 -0.359250 [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] x 1.3380450 1.338045 1.338045 1.338045 1.338045 1.338045 1.338045 1.338045 y 0.8531036 1.338045 1.338045 1.338045 1.338045 1.338045 1.338045 1.338045 [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] x 1.338045 1.338045 3.365152 3.365152 3.365152 2.936498 2.936498 0.3065824 y 1.338045 1.946177 2.959731 3.365152 3.236556 3.022229 2.147524 0.8325656 [,25] [,26] [,27] [,28] [,29] [,30] x 0.3065824 0.3065824 0.3065824 0.3065824 0.3065824 0.3065824 y 0.3065824 0.3065824 0.3065824 0.3065824 0.3065824 0.3065824 > s <- stepcand(y, family = "gaussKern", param = k, cand.radius = Inf) > s <- s[cumsum(bl),] > s[,] leftEnd rightEnd value leftIndex rightIndex cumSum cumSumSq 6 1 6 0.4521269 1 6 6.010166 6.382779 8 7 8 -0.3592500 7 8 5.003908 6.930459 18 9 18 1.9461772 9 18 18.507549 25.768766 21 19 21 3.2365561 19 21 28.068988 56.328318 23 22 23 2.1475235 22 23 33.238741 70.074044 30 24 30 0.3065824 24 30 35.910801 71.331166 cumSumWe acumSum acumSumSq acumSumWe bcumSum bcumSumSq bcumSumWe number 6 6 8.533102 11.23897 11 1.111608 1.235672 1 NA 8 8 11.209192 14.81970 13 3.334824 3.707016 3 NA 18 18 33.238741 70.07404 23 11.209192 14.819703 13 1 21 21 34.684471 70.95519 26 15.223327 20.190796 16 NA 23 23 35.297636 71.14318 28 18.507549 25.768766 18 NA 30 30 NA NA NA 34.377889 70.861202 25 0 improve lXy lcXy rcXy rXy 6 NA 0.111885 5.857012 4.446432 1.2986950 8 NA 3.250798 8.533102 6.010166 -0.7278625 18 18.29468 9.553346 28.068988 15.223327 3.2923152 21 NA 6.368839 34.071306 21.467280 6.2351873 23 NA 1.923474 34.684471 28.068988 4.6920087 30 0.00000 0.000000 NA 35.910801 0.0000000 > sre <- s[refit = y] > sre[,] leftEnd rightEnd value leftIndex rightIndex cumSum cumSumSq 6 1 6 1.1116079 1 6 6.010166 6.382779 8 7 8 -1.0866622 7 8 5.003908 6.930459 18 9 18 1.3380450 9 18 18.507549 25.768766 21 19 21 3.3651523 19 21 28.068988 56.328318 23 22 23 2.9364983 22 23 33.238741 70.074044 30 24 30 0.3065824 24 30 35.910801 71.331166 cumSumWe acumSum acumSumSq acumSumWe bcumSum bcumSumSq bcumSumWe number 6 6 8.533102 11.23897 11 1.111608 1.235672 1 NA 8 8 11.209192 14.81970 13 3.334824 3.707016 3 NA 18 18 33.238741 70.07404 23 11.209192 14.819703 13 1 21 21 34.684471 70.95519 26 15.223327 20.190796 16 NA 23 23 35.297636 71.14318 28 18.507549 25.768766 18 NA 30 30 NA NA NA 34.377889 70.861202 25 0 improve lXy lcXy rcXy rXy 6 NA 0.111885 5.857012 4.446432 1.2986950 8 NA 3.250798 8.533102 6.010166 -0.7278625 18 18.29468 9.553346 28.068988 15.223327 3.2923152 21 NA 6.368839 34.071306 21.467280 6.2351873 23 NA 1.923474 34.684471 28.068988 4.6920087 30 0.00000 0.000000 NA 35.910801 0.0000000 > stopifnot(all.eq(sre$value, bh)) > > # check Bessel filters (and hence polynomials), cf. Bond__BesselFiltConst.pdf > for(pole in 1:6) { + bf <- dfilter(param = list(pole = pole, cutoff = runif(1, 5e-2, 3e-1))) + print(bf) + print(bf$param$omega0) + # check if length 2 / cutoff is long enough + stopifnot(round(bf$step[length(bf$step)], 6) == 1) + # check coefficients of polynom + stopifnot(all.eq(bf$param$a, list( + c(1,1), + c(3, 3, 1), + c(15, 15, 6, 1), + c(105, 105, 45, 10, 1), + c(945, 945, 420, 105, 15, 1), + c(10395, 10395, 4725, 1260, 210, 21, 1) + )[[pole]])) + # check whether spectrum is halved at cutoff + stopifnot(round(bf$param$spectrum(bf$param$cutoff) - 0.5, 12) == 0) + # check frequency normalisation constant omega0 + stopifnot(round(bf$param$omega0 - c( + 1, + 1.361654129, + 1.755672389, + 2.113917675, + 2.427410702, + 2.703395061 + )[pole], 7) == 0) + # check if kernel normalises to 1 + stopifnot(round(sum(bf$param$kernfun(seq(0, 3 / bf$param$cutoff, by = 1e-3))) * 1e-3, 2) == 1) + } Digitised 1-pole Bessel filter cut-off frequency: 0.2651548 length : 12 jump after: 1 step response start : 0 1 - step response end: 1.099105e-08 [1] 1 Digitised 2-pole Bessel filter cut-off frequency: 0.1514416 length : 20 jump after: 2 step response start : 0 1 - step response end: -2.3213e-09 [1] 1.361654 Digitised 3-pole Bessel filter cut-off frequency: 0.1294471 length : 24 jump after: 3 step response start : 0 1 - step response end: -3.635286e-09 [1] 1.755672 Digitised 4-pole Bessel filter cut-off frequency: 0.1855463 length : 17 jump after: 2 step response start : 0 1 - step response end: 2.103653e-09 [1] 2.113918 Digitised 5-pole Bessel filter cut-off frequency: 0.2447632 length : 13 jump after: 2 step response start : 0 1 - step response end: -1.821879e-08 [1] 2.427411 Digitised 6-pole Bessel filter cut-off frequency: 0.05928291 length : 51 jump after: 8 step response start : 0 1 - step response end: 3.150636e-08 [1] 2.703395 > > # check confidence sets and bands > y <- c(rep(0, 5), rep(5, 1), rep(10, 5), rep(5, 1),rep(0, 5)) > y [1] 0 0 0 0 0 5 10 10 10 10 10 5 0 0 0 0 0 > b <- bounds(y, param = 1, pen="sqrt", q=1) > b $bounds li ri lower upper 1 1 1 -3.76883129 3.768831 2 1 2 -2.47912996 2.479130 3 1 4 -1.60610103 1.606101 4 1 8 2.10929636 4.140704 5 1 16 3.13588727 4.364113 6 2 2 -3.76883129 3.768831 7 2 3 -2.47912996 2.479130 8 2 5 -1.60610103 1.606101 9 2 9 3.35929636 5.390704 10 2 17 3.13588727 4.364113 11 3 3 -3.76883129 3.768831 12 3 4 -2.47912996 2.479130 13 3 6 -0.35610103 2.856101 14 3 10 4.60929636 6.640704 15 4 4 -3.76883129 3.768831 16 4 5 -2.47912996 2.479130 17 4 7 2.14389897 5.356101 18 4 11 5.85929636 7.890704 19 5 5 -3.76883129 3.768831 20 5 6 0.02087004 4.979130 21 5 8 4.64389897 7.856101 22 5 12 6.48429636 8.515704 23 6 6 1.23116871 8.768831 24 6 7 5.02087004 9.979130 25 6 9 7.14389897 10.356101 26 6 13 6.48429636 8.515704 27 7 7 6.23116871 13.768831 28 7 8 7.52087004 12.479130 29 7 10 8.39389897 11.606101 30 7 14 5.85929636 7.890704 31 8 8 6.23116871 13.768831 32 8 9 7.52087004 12.479130 33 8 11 8.39389897 11.606101 34 8 15 4.60929636 6.640704 35 9 9 6.23116871 13.768831 36 9 10 7.52087004 12.479130 37 9 12 7.14389897 10.356101 38 9 16 3.35929636 5.390704 39 10 10 6.23116871 13.768831 40 10 11 7.52087004 12.479130 41 10 13 4.64389897 7.856101 42 10 17 2.10929636 4.140704 43 11 11 6.23116871 13.768831 44 11 12 5.02087004 9.979130 45 11 14 2.14389897 5.356101 46 12 12 1.23116871 8.768831 47 12 13 0.02087004 4.979130 48 12 15 -0.35610103 2.856101 49 13 13 -3.76883129 3.768831 50 13 14 -2.47912996 2.479130 51 13 16 -1.60610103 1.606101 52 14 14 -3.76883129 3.768831 53 14 15 -2.47912996 2.479130 54 14 17 -1.60610103 1.606101 55 15 15 -3.76883129 3.768831 56 15 16 -2.47912996 2.479130 57 16 16 -3.76883129 3.768831 58 16 17 -2.47912996 2.479130 59 17 17 -3.76883129 3.768831 $start [1] 0 5 10 14 18 22 26 30 34 38 42 45 48 51 54 56 58 $feasible [1] TRUE attr(,"class") [1] "bounds" "list" > sb <- stepbound(y, b, conf.bands = TRUE) > as.data.frame(sb) leftEnd rightEnd value leftIndex rightIndex rightIndexLeftBound 1 1 5 0.000000 1 5 5 2 6 12 8.571429 6 12 11 3 13 17 0.000000 13 17 17 rightIndexRightBound rightEndLeftBound rightEndRightBound cumSum cumSumWe 1 6 5 6 0 5 2 12 11 12 60 12 3 17 17 17 60 17 cumSumSq 1 0 2 550 3 550 > attr(sb, "conf.bands") NULL > # check confidence intervals > stopifnot(round(sb$rightIndexRightBound - c( + 6,12,17 + ), 0) == 0) > stopifnot(round(sb$rightIndexLeftBound - c( + 5,11,17 + ), 0) == 0) > attr(sb,"conf.band") NULL > # check confidence bands > stopifnot(round(attr(sb,"conf.band")$lower - c( + rep(-1.606101,5),1.231169, rep(8.393899,5), 1.231169, rep(-1.606101,5) + ), 4) == 0) > stopifnot(round(attr(sb,"conf.band")$upper - c( + rep(1.606101,5),8.768831, rep(11.606101,5), 8.768831, rep(1.606101,5) + ), 4) == 0) > > > # check if any warnings were produced > warnings() > stopifnot(length(warnings()) == 0) > > proc.time() user system elapsed 1.51 0.10 1.62