[Mboost-commits] r752 - in pkg/mboostDevel: R tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Oct 25 14:33:29 CEST 2013
Author: hofner
Date: 2013-10-25 14:33:29 +0200 (Fri, 25 Oct 2013)
New Revision: 752
Added:
pkg/mboostDevel/tests/regtest-inference.R
Modified:
pkg/mboostDevel/R/bmono.R
pkg/mboostDevel/R/inference.R
Log:
- changed default method in bmono
- added new experimental error.bound to stabsel
- added new tests for stabsel
Modified: pkg/mboostDevel/R/bmono.R
===================================================================
--- pkg/mboostDevel/R/bmono.R 2013-10-17 12:21:13 UTC (rev 751)
+++ pkg/mboostDevel/R/bmono.R 2013-10-25 12:33:29 UTC (rev 752)
@@ -2,7 +2,7 @@
bmono <- function(..., constraint = c("increasing", "decreasing",
"convex", "concave", "none",
"positive", "negative"),
- type = c("iterative", "quad.prog"),
+ type = c("quad.prog", "iterative"),
by = NULL, index = NULL, knots = 20, boundary.knots = NULL,
degree = 3, differences = 2, df = 4,
lambda = NULL, lambda2 = 1e6, niter = 10,
Modified: pkg/mboostDevel/R/inference.R
===================================================================
--- pkg/mboostDevel/R/inference.R 2013-10-17 12:21:13 UTC (rev 751)
+++ pkg/mboostDevel/R/inference.R 2013-10-25 12:33:29 UTC (rev 752)
@@ -1,11 +1,15 @@
stabsel <- function(object, cutoff, q, PFER,
- folds = cv(model.weights(object), type = "subsampling", B = 100),
- papply = mclapply, verbose = TRUE, FWER, ...) {
+ folds = cv(model.weights(object), type = "subsampling",
+ B = ifelse(error.bound = "MB", 100, 50)),
+ papply = mclapply, verbose = TRUE, FWER,
+ error.bound = c("MB", "SS"), ...) {
p <- length(variable.names(object))
ibase <- 1:p
+ error.bound <- match.args(error.bound)
+
## only two of the four arguments can be specified
if ((nmiss <- sum(missing(PFER), missing(cutoff),
missing(q), missing(FWER))) != 2) {
@@ -45,8 +49,18 @@
}
if (missing(cutoff)) {
- cutoff <- min(0.9, tmp <- (q^2 / (PFER * p) + 1) / 2)
- upperbound <- q^2 / p / (2 * cutoff - 1)
+ if (error.bound == "MB") {
+ cutoff <- min(0.9, tmp <- (q^2 / (PFER * p) + 1) / 2)
+ upperbound <- q^2 / p / (2 * cutoff - 1)
+ } else {
+ objective <- function(cutoff) {
+ PFER / p - minD(q, p, cutoff, B)
+ }
+ root <- uniroot(objective, lower = 0.5, upper = 0.9)$root
+ cutoff <- floor(root * 2 * B) / (2* B)
+ upperbound <- minD(q, p, cutoff, B) * p
+ }
+ upperbound <- signif(upperbound, 3)
if (verbose && tmp > 0.9 && upperbound - PFER > PFER/2) {
warning("Upper bound for PFER > ", PFER,
" for the given value of ", sQuote("q"),
@@ -55,8 +69,20 @@
}
if (missing(q)) {
- q <- ceiling(sqrt(PFER * (2 * cutoff - 1) * p))
- upperbound <- q^2 / p / (2 * cutoff - 1)
+ if (error.bound == "MB") {
+ q <- ceiling(sqrt(PFER * (2 * cutoff - 1) * p))
+ upperbound <- q^2 / p / (2 * cutoff - 1)
+ } else {
+ objective <- function(q) {
+ PFER / p - minD(q, p, cutoff, B)
+ }
+ root <- uniroot(objective, lower = 1,
+ upper = min(sqrt((B - 1) / (2 * B) * p^2),
+ (B - 1) / (2 * B) * p))$root
+ q <- ceiling(root)
+ upperbound <- minD(q, p, cutoff, B) * p
+ }
+ upperbound <- signif(upperbound, 3)
if (verbose && upperbound - PFER > PFER/2)
warning("Upper bound for PFER > ", PFER,
" for the given value of ", sQuote("cutoff"),
@@ -64,7 +90,12 @@
}
if (missing(PFER)) {
- upperbound <- PFER <- q^2 / p / (2 * cutoff - 1)
+ if (error.bound == "MB") {
+ upperbound <- PFER <- q^2 / p / (2 * cutoff - 1)
+ } else {
+ upperbound <- PFER <- minD(q, p, cutoff, B) * p
+ }
+ upperbound <- signif(upperbound, 3)
}
if (verbose && PFER >= p)
warning("Upper bound for PFER larger than the number of base-learners.")
@@ -75,9 +106,13 @@
xs[qq > q] <- xs[1]
xs
}
- ss <- cvrisk(object, fun = fun,
+ if (error.bound == "SS") {
+ ## use complementary pairs
+ folds <- cbind(folds, model.weights(object) - folds)
+ }
+ ss <- cvrisk(object, fun = fun,
folds = folds,
- papply = papply, ...)
+ papply = papply, ....)
if (verbose){
qq <- sapply(ss, function(x) length(unique(x)))
@@ -161,3 +196,71 @@
ret <- abs(ret) / length(ss)
ret
}
+
+
+### Modified version of the code accompanying the paper:
+### Shah, R. D. and Samworth, R. J. (2013), Variable selection with error
+### control: Another look at Stability Selection, J. Roy. Statist. Soc., Ser.
+### B, 75, 55-80. DOI: 10.1111/j.1467-9868.2011.01034.x
+###
+### Original code available from
+### http://www.statslab.cam.ac.uk/~rds37/papers/r_concave_tail.R
+### or
+### http://www.statslab.cam.ac.uk/~rjs57/r_concave_tail.R
+
+D <- function(theta, which, B, r) {
+ ## If pi = ceil{ B * 2 * eta} / B + 1/B,..., 1 return the tail probability.
+ ## If pi < ceil{ B * 2 * eta} / B return 1
+
+ if (which <= 0)
+ return(1)
+ ### pi muss ein vielfaches von 1/(2 * B) sein, oder?
+
+ s <- 1/r
+ thetaB <- theta * B
+ k_start <- (ceiling(2 * thetaB) + 1)
+ if(k_start > B)
+ stop("theta to large")
+
+ Find.a <- function(prev_a)
+ uniroot(Calc.a, lower = 0.0001, upper = prev_a,
+ tol = .Machine$double.eps^0.75)$root
+
+ Calc.a <- function(a) {
+ denom <- sum((a + 0:k)^s)
+ num <- sum((0:k) * (a + 0:k)^s)
+ num / denom - thetaB
+ }
+
+ OptimInt <- function(a) {
+ num <- (k + 1 - thetaB) * sum((a + 0:(t-1))^s)
+ denom <- sum((k + 1 - (0:k)) * (a + 0:k)^s)
+ 1 - num / denom
+ }
+
+ ## initialize a
+ a_vec <- rep(100000, B)
+
+ ## compute a values
+ for(k in k_start:B)
+ a_vec[k] <- Find.a(a_vec[k-1])
+
+ t <- which
+ cur_optim <- rep(0, B)
+ for (k in k_start:(B-1))
+ cur_optim[k] <- optimize(f=OptimInt, lower = a_vec[k+1],
+ upper = a_vec[k], maximum = TRUE)$objective
+ return(max(cur_optim))
+}
+
+minD <- function(q, p, pi, B, r = c(-1/2, -1/4)) {
+ which <- ceiling(signif(pi / (1/(2* B)), 10))
+ #if ((which) %% 1 != 0)
+ # stop(sQuote("pi"), " must be a multiple of 1/(2 * B)")
+
+ maxQ <- min(sqrt((B - 1) / (2 * B) * p^2),
+ (B - 1) / (2 * B) * p)
+ if (q > maxQ)
+ stop(sQuote("q"), " must be <= ", maxQ)
+ min(c(1, D(q^2 / p^2, which - B, B, r[1]), D(q / p, which , 2*B, r[2])))
+}
Added: pkg/mboostDevel/tests/regtest-inference.R
===================================================================
--- pkg/mboostDevel/tests/regtest-inference.R (rev 0)
+++ pkg/mboostDevel/tests/regtest-inference.R 2013-10-25 12:33:29 UTC (rev 752)
@@ -0,0 +1,131 @@
+require("mboostDevel")
+attach(asNamespace("mboostDevel"))
+
+### (Slightly) modified version of the code accompanying the paper:
+### Shah, R. D. and Samworth, R. J. (2013), Variable selection with error
+### control: Another look at Stability Selection, J. Roy. Statist. Soc., Ser.
+### B, 75, 55-80. DOI: 10.1111/j.1467-9868.2011.01034.x
+###
+### Original code available from
+### http://www.statslab.cam.ac.uk/~rds37/papers/r_concave_tail.R
+### or
+### http://www.statslab.cam.ac.uk/~rjs57/r_concave_tail.R
+r.TailProbs <- function(eta, B, r) {
+ ## If pi = ceil{ B * 2 * eta} / B + 1/B,..., 1 return the tail probability.
+ ## If pi < ceil{ B * 2 * eta} / B return 1
+
+ MAXa <- 100000
+ MINa <- 0.0001
+
+ s <- -1/r
+ etaB <- eta * B
+ k_start <- (ceiling(2 * etaB) + 1)
+ if(k_start > B)
+ stop("eta is too large")
+
+ a_vec <- rep(MAXa,B)
+
+ Find.a <- function(prev_a)
+ uniroot(Calc.a, lower = MINa, upper = prev_a,
+ tol = .Machine$double.eps^0.75)$root
+
+ Calc.a <- function(a) {
+ denom <- sum((a + 0:k)^(-s))
+ num <- sum((0:k) * (a + 0:k)^(-s))
+ num / denom - etaB
+ }
+
+ for(k in k_start:B)
+ a_vec[k] <- Find.a(a_vec[k-1])
+
+ OptimInt <- function(a) {
+ num <- (k + 1 - etaB) * sum((a + 0:(t-1))^(-s))
+ denom <- sum((k + 1 - (0:k)) * (a + 0:k)^(-s))
+ 1 - num / denom
+ }
+
+ output <- rep(1, B)
+
+ prev_k <- k_start
+ for(t in k_start:B) {
+ cur_optim <- rep(0, B)
+ for (k in prev_k:(B-1))
+ cur_optim[k] <- optimize(f=OptimInt, lower = a_vec[k+1],
+ upper = a_vec[k], maximum = TRUE)$objective
+ output[t] <- max(cur_optim)
+ prev_k <- which.max(cur_optim)
+ }
+ return(output)
+}
+
+pminD <- function(theta, B, r = c(-1/2, -1/4)) {
+ pmin(c(rep(1, B), r.TailProbs(theta^2, B, r[1])),
+ r.TailProbs(theta, 2*B, r[2]))
+}
+
+## test r-concave bound
+B <- 50
+x <- (1:(2 * B))/(2 * B)
+p <- 1000
+q <- 50
+theta <- q/p
+if (FALSE) {
+ plot(x, log(pminD(theta, B)), xlab = "pi")
+ abline(v = ceiling(2 * theta * 100) / 100)
+ Ds <- cbind(c(rep(1, B), r.TailProbs(theta^2, B, -1/2)),
+ r.TailProbs(theta, 2*B, -1/4))
+ round(log(Ds), 2)
+ lines(x, log(Ds[,1]), col = "red", lwd = 2)
+ lines(x, log(Ds[,2]), col = "blue", lwd = 2)
+}
+
+## r-concave bound of Shah & Samworth (2013)
+bound_ss <- (pminD(theta, B) * p)[40:100]
+plot(x[40:100], bound_ss, xlab = "pi", ylim = c(0, 50))
+## Bound of Meinshausen & Buehlmann (2010)
+points(x[40:100], q^2 / (2 * x[40:100] - 1) / p, col = "red")
+## now our implementation
+bound <- rep(NA, 61)
+for (i in 40:100) {
+ bound[i - 39] <- minD(q, p, i/100, B) * p
+}
+points((40:100)/100, bound, col = "green")
+stopifnot(bound == bound_ss)
+
+### computation of q from other values
+cutoff <- 0.6
+PFER <- 0.2
+B <- 50
+p <- 200
+objective <- function(q) {
+ PFER / p - minD(q, p, cutoff, B)
+}
+root <- uniroot(objective, lower = 1,
+ upper = min(sqrt((B - 1) / (2 * B) * p^2),
+ (B - 1) / (2 * B) * p))$root
+(q <- ceiling(root))
+# check:
+round(minD(q - 1, p, cutoff, B) * p, 3)
+round(minD(q, p, cutoff, B) * p, 3)
+
+
+### computation of cutoff from other values
+PFER <- 0.2
+B <- 50
+p <- 200
+q <- 7
+objective <- function(cutoff) {
+ PFER / p - minD(q, p, cutoff, B)
+}
+root <- uniroot(objective, lower = 0.5, upper = 0.9)$root
+(cutoff <- floor(root * 2 * B) / (2* B))
+# check:
+round(minD(q, p, cutoff, B) * p, 3)
+round(minD(q, p, cutoff + 1e-5, B) * p, 3)
+
+
+weights <- c(2, 3, 2, rep(1, 90), rep(0, 10))
+folds <- cv(weights, type = "subsampling", B = 10)
+head(folds)
+tail(folds)
+rowSums(cbind(folds, weights - folds))
More information about the Mboost-commits
mailing list