R Under development (unstable) (2024-10-17 r87242 ucrt) -- "Unsuffered Consequences" Copyright (C) 2024 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-10. 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.73239 > b <- bounds(y, r = 1e2, param = sd, lengths = c(1,4,9)) > b $bounds li ri lower upper 1 1 1 -5.11855989 -2.88144011 2 1 4 -3.05927994 -1.94072006 3 1 9 -0.37285330 0.37285330 4 2 2 -4.11855989 -1.88144011 5 2 5 -2.05927994 -0.94072006 6 3 3 -3.11855989 -0.88144011 7 3 6 -1.05927994 0.05927994 8 4 4 -2.11855989 0.11855989 9 4 7 -0.05927994 1.05927994 10 5 5 -1.11855989 1.11855989 11 5 8 0.94072006 2.05927994 12 6 6 -0.11855989 2.11855989 13 6 9 1.94072006 3.05927994 14 7 7 0.88144011 3.11855989 15 8 8 1.88144011 4.11855989 16 9 9 2.88144011 5.11855989 $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.2115234 -2.7884766 2 1 2 -3.1057617 -1.8942383 3 1 4 -0.4038411 0.4038411 4 2 2 -4.2115234 -1.7884766 6 2 2 -3.2115234 -0.7884766 8 2 2 -2.2115234 0.2115234 5 2 3 -2.1057617 -0.8942383 7 2 4 -1.1057617 0.1057617 9 2 4 -0.1057617 1.1057617 10 3 3 -1.2115234 1.2115234 11 3 4 0.8942383 2.1057617 12 4 4 -0.2115234 2.2115234 13 4 4 1.8942383 3.1057617 14 4 4 0.7884766 3.2115234 15 4 4 1.7884766 4.2115234 16 4 4 2.7884766 5.2115234 $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.11855989 -2.88144011 4 1 1 -4.11855989 -1.88144011 2 1 2 -3.05927994 -1.94072006 5 1 3 -2.05927994 -0.94072006 3 1 6 -0.37285330 0.37285330 6 2 2 -3.11855989 -0.88144011 8 2 2 -2.11855989 0.11855989 7 2 4 -1.05927994 0.05927994 9 2 5 -0.05927994 1.05927994 10 3 3 -1.11855989 1.11855989 11 3 6 0.94072006 2.05927994 12 4 4 -0.11855989 2.11855989 13 4 6 1.94072006 3.05927994 14 5 5 0.88144011 3.11855989 15 6 6 1.88144011 4.11855989 16 6 6 2.88144011 5.11855989 $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 6.000000 8.000000 74.00000000 97.00000 100.000000 f -5.653255 1.683562 -0.09633049 -2.11576 2.515174 lambda 11.363503 23.667162 19.80826399 16.18615 25.719506 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.360078 , 2.214946 ] cost: 63.6128 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: [ 6 , 26.75 ] cost: -3693.576 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.15 , 0.9 ] cost: 668.8083 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.600806 , 2.34758 ] cost: 0 > stopifnot(attr(cand, "cost") == 0) > system.time(stopifnot(length(steppath(cand)) == 100)) user system elapsed 0.02 0.00 0.02 > 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.7537085 , -0.7537085 ] cost: 314.5523 > stopifnot(nrow(print(stepcand(x[1], max.cand = 1))) == 1) Fitted step function of family gauss containing 1 blocks domain: ( NA , 1 ] range: [ -4.458854 , -4.458854 ] 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.2577172 -0.2577172 -0.2577172 -0.2577172 -0.2577172 -0.2577172 -0.8086053 y -0.2577172 -0.2577172 -0.2577172 -0.2577172 -0.2577172 -0.4229836 -0.6984277 [,8] [,9] [,10] [,11] [,12] [,13] [,14] x -0.8086053 0.5161687 0.5161687 0.5161687 0.5161687 0.5161687 0.5161687 y -0.4111731 0.2512139 0.5161687 0.5161687 0.5161687 0.5161687 0.5161687 [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] x 0.5161687 0.5161687 0.5161687 0.5161687 1.266949 1.266949 1.266949 0.7140243 y 0.5161687 0.5161687 0.5161687 0.7414027 1.116793 1.266949 1.101071 0.8246092 [,23] [,24] [,25] [,26] [,27] [,28] x 0.7140243 -0.04404413 -0.04404413 -0.04404413 -0.04404413 -0.04404413 y 0.4866038 0.10756956 -0.04404413 -0.04404413 -0.04404413 -0.04404413 [,29] [,30] x -0.04404413 -0.04404413 y -0.04404413 -0.04404413 > 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.42298362 1 6 -1.711569 0.5110059 8 7 8 -0.41117313 7 8 -2.821170 1.1678704 18 9 18 0.74140267 9 18 2.300795 3.9120974 21 19 21 1.10107139 19 21 5.785608 7.9768406 23 22 23 0.48660379 22 23 7.096821 8.8936042 30 24 30 -0.04404413 24 30 6.940126 8.9168147 cumSumWe acumSum acumSumSq acumSumWe bcumSum bcumSumSq bcumSumWe 6 6 -1.5376191 1.763839 11 -0.2577172 0.06641814 1 8 8 -0.5052818 2.296699 13 -0.7731515 0.19925443 3 18 18 7.0968212 8.893604 23 -0.5052818 2.29669917 13 21 21 7.1163026 8.909055 26 1.0432241 3.09598942 16 23 23 7.0282143 8.912935 28 2.3007955 3.91209742 18 30 30 NA NA NA 7.1603467 8.90711527 25 number improve lXy lcXy rcXy rXy 6 NA NA -0.8455965 -2.569956 -1.030869 -0.6934912 8 NA NA 1.1099565 -1.537619 -1.711569 -0.9360061 18 NA NA 3.4838751 5.785608 1.043224 1.2585091 21 NA NA 1.5841821 7.204391 3.417588 2.2026205 23 NA NA 0.1439485 7.116303 5.785608 1.1867458 30 0 0 0.0000000 NA 6.940126 0.0000000 > sre <- s[refit = y] > sre[,] leftEnd rightEnd value leftIndex rightIndex cumSum cumSumSq 6 1 6 -0.25771718 1 6 -1.711569 0.5110059 8 7 8 -0.80860532 7 8 -2.821170 1.1678704 18 9 18 0.51616866 9 18 2.300795 3.9120974 21 19 21 1.26694871 19 21 5.785608 7.9768406 23 22 23 0.71402432 22 23 7.096821 8.8936042 30 24 30 -0.04404413 24 30 6.940126 8.9168147 cumSumWe acumSum acumSumSq acumSumWe bcumSum bcumSumSq bcumSumWe 6 6 -1.5376191 1.763839 11 -0.2577172 0.06641814 1 8 8 -0.5052818 2.296699 13 -0.7731515 0.19925443 3 18 18 7.0968212 8.893604 23 -0.5052818 2.29669917 13 21 21 7.1163026 8.909055 26 1.0432241 3.09598942 16 23 23 7.0282143 8.912935 28 2.3007955 3.91209742 18 30 30 NA NA NA 7.1603467 8.90711527 25 number improve lXy lcXy rcXy rXy 6 NA NA -0.8455965 -2.569956 -1.030869 -0.6934912 8 NA NA 1.1099565 -1.537619 -1.711569 -0.9360061 18 NA NA 3.4838751 5.785608 1.043224 1.2585091 21 NA NA 1.5841821 7.204391 3.417588 2.2026205 23 NA NA 0.1439485 7.116303 5.785608 1.1867458 30 0 0 0.0000000 NA 6.940126 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.1159441 length : 26 jump after: 1 step response start : 0 1 - step response end: 1.231482e-08 [1] 1 Digitised 2-pole Bessel filter cut-off frequency: 0.1964496 length : 16 jump after: 1 step response start : 0 1 - step response end: -7.321632e-10 [1] 1.361654 Digitised 3-pole Bessel filter cut-off frequency: 0.1462542 length : 21 jump after: 2 step response start : 0 1 - step response end: -6.836724e-09 [1] 1.755672 Digitised 4-pole Bessel filter cut-off frequency: 0.2345528 length : 13 jump after: 2 step response start : 0 1 - step response end: 3.680064e-08 [1] 2.113918 Digitised 5-pole Bessel filter cut-off frequency: 0.2235937 length : 14 jump after: 2 step response start : 0 1 - step response end: -2.911387e-08 [1] 2.427411 Digitised 6-pole Bessel filter cut-off frequency: 0.2978159 length : 11 jump after: 2 step response start : 0 1 - step response end: 2.622409e-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.65 0.17 1.81