R Under development (unstable) (2023-07-19 r84711 ucrt) -- "Unsuffered Consequences" Copyright (C) 2023 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-7. 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.809084 > b <- bounds(y, r = 1e2, param = sd, lengths = c(1,4,9)) > b $bounds li ri lower upper 1 1 1 -5.17417003 -2.82582997 2 1 4 -3.08708502 -1.91291498 3 1 9 -0.39139001 0.39139001 4 2 2 -4.17417003 -1.82582997 5 2 5 -2.08708502 -0.91291498 6 3 3 -3.17417003 -0.82582997 7 3 6 -1.08708502 0.08708502 8 4 4 -2.17417003 0.17417003 9 4 7 -0.08708502 1.08708502 10 5 5 -1.17417003 1.17417003 11 5 8 0.91291498 2.08708502 12 6 6 -0.17417003 2.17417003 13 6 9 1.91291498 3.08708502 14 7 7 0.82582997 3.17417003 15 8 8 1.82582997 4.17417003 16 9 9 2.82582997 5.17417003 $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.17417003 -2.82582997 2 1 2 -3.08708502 -1.91291498 3 1 4 -0.39139001 0.39139001 4 2 2 -4.17417003 -1.82582997 6 2 2 -3.17417003 -0.82582997 8 2 2 -2.17417003 0.17417003 5 2 3 -2.08708502 -0.91291498 7 2 4 -1.08708502 0.08708502 9 2 4 -0.08708502 1.08708502 10 3 3 -1.17417003 1.17417003 11 3 4 0.91291498 2.08708502 12 4 4 -0.17417003 2.17417003 13 4 4 1.91291498 3.08708502 14 4 4 0.82582997 3.17417003 15 4 4 1.82582997 4.17417003 16 4 4 2.82582997 5.17417003 $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.17417003 -2.82582997 4 1 1 -4.17417003 -1.82582997 2 1 2 -3.08708502 -1.91291498 5 1 3 -2.08708502 -0.91291498 3 1 6 -0.39139001 0.39139001 6 2 2 -3.17417003 -0.82582997 8 2 2 -2.17417003 0.17417003 7 2 4 -1.08708502 0.08708502 9 2 5 -0.08708502 1.08708502 10 3 3 -1.17417003 1.17417003 11 3 6 0.91291498 2.08708502 12 4 4 -0.17417003 2.17417003 13 4 6 1.91291498 3.08708502 14 5 5 0.82582997 3.17417003 15 6 6 1.82582997 4.17417003 16 6 6 2.82582997 5.17417003 $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 11.000000 21.000000 27.000000 45.000000 100.00000 f -4.428331 3.425267 -3.514253 -7.664246 -2.15797 lambda 12.844288 28.170038 14.073688 9.293430 16.11798 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: [ -7.692309 , 3.897414 ] cost: 64.36299 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: [ 8.95 , 30 ] cost: -2914.387 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.2 , 1 ] cost: 628.5494 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: [ -9.405179 , 4.4646 ] 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: [ -2.97145 , -2.97145 ] cost: 955.9334 > stopifnot(nrow(print(stepcand(x[1], max.cand = 1))) == 1) Fitted step function of family gauss containing 1 blocks domain: ( NA , 1 ] range: [ -5.332991 , -5.332991 ] 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] x 0.2676046 0.2676046 0.2676046 0.2676046 0.2676046 0.26760461 -0.4416904 y 0.2676046 0.2676046 0.2676046 0.2676046 0.2676046 0.05481611 -0.2998314 [,8] [,9] [,10] [,11] [,12] [,13] [,14] x -0.44169038 0.8576839 0.8576839 0.8576839 0.8576839 0.8576839 0.8576839 y -0.05187811 0.5978090 0.8576839 0.8576839 0.8576839 0.8576839 0.8576839 [,15] [,16] [,17] [,18] [,19] [,20] [,21] x 0.8576839 0.8576839 0.8576839 0.8576839 0.7047178 0.7047178 0.7047178 y 0.8576839 0.8576839 0.8576839 0.8117940 0.7353110 0.7047178 0.6868215 [,22] [,23] [,24] [,25] [,26] [,27] [,28] x 0.6450634 0.6450634 -0.6313149 -0.6313149 -0.6313149 -0.6313149 -0.6313149 y 0.6569943 0.2621499 -0.3760393 -0.6313149 -0.6313149 -0.6313149 -0.6313149 [,29] [,30] x -0.6313149 -0.6313149 y -0.6313149 -0.6313149 > 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.05481611 1 6 1.392839 0.3610659 8 7 8 -0.05187811 7 8 1.041130 0.4536561 18 9 18 0.81179402 9 18 9.312203 7.3550140 21 19 21 0.68682145 19 21 11.439054 8.8640470 23 22 23 0.26214991 22 23 12.358198 9.3644111 30 24 30 -0.63131492 24 30 8.194269 11.8971678 cumSumWe acumSum acumSumSq acumSumWe bcumSum bcumSumSq bcumSumWe 6 6 3.354306 2.282275 11 0.2676046 0.07161223 1 8 8 5.069674 3.753518 13 0.8028138 0.21483668 3 18 18 12.358198 9.364411 23 5.0696741 3.75351809 13 21 21 10.719529 10.302934 26 7.6427256 5.96038286 16 23 23 9.456899 11.100051 28 9.3122035 7.35501398 18 30 30 NA NA NA 11.3508437 9.90437515 25 number improve lXy lcXy rcXy rXy 6 2 1.776472 0.3225106 1.638939 1.070418 0.2460096 8 NA NA 2.1780515 3.354306 1.392839 -0.2165843 18 1 2.941654 2.2233262 11.439054 7.642726 1.5730019 21 NA NA 0.6177525 11.982159 10.047514 1.3168916 23 NA NA -1.4848163 10.719529 11.439054 0.7652914 30 0 0.000000 0.0000000 NA 8.194269 0.0000000 > sre <- s[refit = y] > sre[,] leftEnd rightEnd value leftIndex rightIndex cumSum cumSumSq 6 1 6 0.2676046 1 6 1.392839 0.3610659 8 7 8 -0.4416904 7 8 1.041130 0.4536561 18 9 18 0.8576839 9 18 9.312203 7.3550140 21 19 21 0.7047178 19 21 11.439054 8.8640470 23 22 23 0.6450634 22 23 12.358198 9.3644111 30 24 30 -0.6313149 24 30 8.194269 11.8971678 cumSumWe acumSum acumSumSq acumSumWe bcumSum bcumSumSq bcumSumWe 6 6 3.354306 2.282275 11 0.2676046 0.07161223 1 8 8 5.069674 3.753518 13 0.8028138 0.21483668 3 18 18 12.358198 9.364411 23 5.0696741 3.75351809 13 21 21 10.719529 10.302934 26 7.6427256 5.96038286 16 23 23 9.456899 11.100051 28 9.3122035 7.35501398 18 30 30 NA NA NA 11.3508437 9.90437515 25 number improve lXy lcXy rcXy rXy 6 2 1.776472 0.3225106 1.638939 1.070418 0.2460096 8 NA NA 2.1780515 3.354306 1.392839 -0.2165843 18 1 2.941654 2.2233262 11.439054 7.642726 1.5730019 21 NA NA 0.6177525 11.982159 10.047514 1.3168916 23 NA NA -1.4848163 10.719529 11.439054 0.7652914 30 0 0.000000 0.0000000 NA 8.194269 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.1074199 length : 28 jump after: 2 step response start : 0 1 - step response end: 1.218145e-08 [1] 1 Digitised 2-pole Bessel filter cut-off frequency: 0.1122112 length : 27 jump after: 2 step response start : 0 1 - step response end: -1.272462e-09 [1] 1.361654 Digitised 3-pole Bessel filter cut-off frequency: 0.2468525 length : 13 jump after: 2 step response start : 0 1 - step response end: -4.451979e-09 [1] 1.755672 Digitised 4-pole Bessel filter cut-off frequency: 0.09416214 length : 32 jump after: 4 step response start : 0 1 - step response end: 1.004753e-08 [1] 2.113918 Digitised 5-pole Bessel filter cut-off frequency: 0.1429999 length : 21 jump after: 3 step response start : 0 1 - step response end: -4.718401e-08 [1] 2.427411 Digitised 6-pole Bessel filter cut-off frequency: 0.08685868 length : 35 jump after: 5 step response start : 0 1 - step response end: 3.562814e-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 > if(!is.null(warnings())) warnings() > stopifnot(is.null(warnings())) > > proc.time() user system elapsed 1.17 0.09 1.25