[Mboost-commits] r748 - in pkg/mboostDevel: . R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Oct 8 16:42:50 CEST 2013
Author: hofner
Date: 2013-10-08 16:42:50 +0200 (Tue, 08 Oct 2013)
New Revision: 748
Modified:
pkg/mboostDevel/DESCRIPTION
pkg/mboostDevel/NAMESPACE
pkg/mboostDevel/R/bmono.R
pkg/mboostDevel/R/helpers.R
Log:
- first experimental version of monotonic P-splines fitted with lsei
(least squares with equality and inequality constraints)
Modified: pkg/mboostDevel/DESCRIPTION
===================================================================
--- pkg/mboostDevel/DESCRIPTION 2013-10-08 14:23:03 UTC (rev 747)
+++ pkg/mboostDevel/DESCRIPTION 2013-10-08 14:42:50 UTC (rev 748)
@@ -16,7 +16,7 @@
trees as base-learners for fitting generalized linear, additive
and interaction models to potentially high-dimensional data.
Depends: R (>= 2.14.0), methods, stats, parallel
-Imports: Matrix, survival, splines, lattice, nnls
+Imports: Matrix, survival, splines, lattice, nnls, limSolve
Suggests: party (>= 1.0-3), TH.data, MASS, fields, BayesX, gbm, mlbench,
RColorBrewer, rpart (>= 4.0-3)
LazyData: yes
Modified: pkg/mboostDevel/NAMESPACE
===================================================================
--- pkg/mboostDevel/NAMESPACE 2013-10-08 14:23:03 UTC (rev 747)
+++ pkg/mboostDevel/NAMESPACE 2013-10-08 14:42:50 UTC (rev 748)
@@ -7,6 +7,7 @@
importFrom(splines, bs, splineDesign)
importFrom(lattice, levelplot)
importFrom(nnls, nnls)
+importFrom(limSolve, lsei)
export(glmboost,
gamboost,
Modified: pkg/mboostDevel/R/bmono.R
===================================================================
--- pkg/mboostDevel/R/bmono.R 2013-10-08 14:23:03 UTC (rev 747)
+++ pkg/mboostDevel/R/bmono.R 2013-10-08 14:42:50 UTC (rev 748)
@@ -2,6 +2,7 @@
bmono <- function(..., constraint = c("increasing", "decreasing",
"convex", "concave", "none",
"positive", "negative"),
+ type = c("iterative", "lsei"),
by = NULL, index = NULL, knots = 20, boundary.knots = NULL,
degree = 3, differences = 2, df = 4,
lambda = NULL, lambda2 = 1e6, niter = 10,
@@ -33,6 +34,8 @@
}
##
+ type = match.arg(type)
+
if (length(mf) == 1 && (is.matrix(mf[[1]]) || is.data.frame(mf[[1]]))) {
mf <- as.data.frame(mf[[1]])
} else {
@@ -105,10 +108,11 @@
degree = degree, differences = differences,
df = df, lambda = lambda, center = FALSE)
args$constraint <- constraint
+ args$type <- type
args$lambda2 <- lambda2
args$niter <- niter
args$boundary.constraints <- boundary.constraints
- if(boundary.constraints){
+ if(boundary.constraints) {
if (is.null(cons.arg$n)){
## use 10% of the knots on each side per default
cons.arg$n <- sapply(args$knots,
@@ -137,6 +141,7 @@
intercept = intercept,
contrasts.arg = contrasts.arg)
args$constraint <- constraint
+ args$type <- type
args$lambda2 <- lambda2
args$niter <- niter
## <FIXME> Was machen wir bei cat. Effekten? Da müsste das doch auch gehen!
@@ -169,7 +174,7 @@
diff_order <- which_diff(args$constraint)
- D <- V <- lambda2 <- vector(mode = "list", length =2)
+ D <- V <- lambda2 <- vector(mode = "list", length = 2)
## set up difference matrix
if (diff_order > 0) {
@@ -296,32 +301,37 @@
y <- y * weights
}
- for (i in 1:args$niter){
- coef <- mysolve(y, V)
- ## compare old and new V
- tmp1 <- do.call(args$constraint[[1]],
- args=list(D[[1]] %*% coef))
- if (lambda2[[2]] != 0)
- tmp2 <- do.call(args$constraint[[2]],
- args=list(D[[2]] %*% coef))
+ if (args$type == "iterative") {
+ for (i in 1:args$niter){
+ coef <- mysolve(y, V)
+ ## compare old and new V
+ tmp1 <- do.call(args$constraint[[1]],
+ args=list(D[[1]] %*% coef))
+ if (lambda2[[2]] != 0)
+ tmp2 <- do.call(args$constraint[[2]],
+ args=list(D[[2]] %*% coef))
- if ( all( V[[1]] == tmp1 ) &&
- ( lambda2[[2]] == 0 || all( V[[2]] == tmp2 ) ) )
- break # if both are equal: done!
- #if (args$boundary.constraints &&
- # all( V[[1]][-idxB, -idxB] == tmp1[-idxB, -idxB]) )
- # break # if both are equal (without V for boundary
- # # constraints): done!
- V[[1]] <- tmp1
- #if (args$boundary.constraints) {
- # V[[1]][idxFlat, idxFlat] <- diag(rep(1, length(idxFlat)))
- #}
- if (lambda2[[2]] != 0)
- V[[2]] <- tmp2
- if (i == args$niter)
- warning("no convergence of coef in bmono\n",
- "You could try increasing ", sQuote("niter"),
- " or ", sQuote("lambda2"))
+ if ( all( V[[1]] == tmp1 ) &&
+ ( lambda2[[2]] == 0 || all( V[[2]] == tmp2 ) ) )
+ break # if both are equal: done!
+ #if (args$boundary.constraints &&
+ # all( V[[1]][-idxB, -idxB] == tmp1[-idxB, -idxB]) )
+ # break # if both are equal (without V for boundary
+ # # constraints): done!
+ V[[1]] <- tmp1
+ #if (args$boundary.constraints) {
+ # V[[1]][idxFlat, idxFlat] <- diag(rep(1, length(idxFlat)))
+ #}
+ if (lambda2[[2]] != 0)
+ V[[2]] <- tmp2
+ if (i == args$niter)
+ warning("no convergence of coef in bmono\n",
+ "You could try increasing ", sQuote("niter"),
+ " or ", sQuote("lambda2"))
+ }
+ } else {
+ coef <- solveLSEI(XtX, crossprod(X, y),
+ constraint = args$constraint)
}
ret <- list(model = coef,
@@ -396,7 +406,7 @@
diag(c(as.numeric(diffs)) <= 0)
negative <- function(diffs)
- diag(c(as.numeric(diffs)) <= 0)
+ diag(c(as.numeric(diffs)) >= 0)
increasing <- function(diffs)
diag(c(as.numeric(diffs)) <= 0)
@@ -414,6 +424,7 @@
if (length(constraint) == 1) {
diff <- switch(constraint,
positive = 0,
+ negative = 0,
increasing = 1,
decreasing = 1,
convex = 2,
@@ -422,6 +433,7 @@
diff <- lapply(constraint, function(x)
switch(x,
positive = 0,
+ negative = 0,
increasing = 1,
decreasing = 1,
convex = 2,
Modified: pkg/mboostDevel/R/helpers.R
===================================================================
--- pkg/mboostDevel/R/helpers.R 2013-10-08 14:23:03 UTC (rev 747)
+++ pkg/mboostDevel/R/helpers.R 2013-10-08 14:42:50 UTC (rev 748)
@@ -159,11 +159,60 @@
XWY <- as.vector(crossprod(X$X1, Y) %*% X$X2)
cf <- nnls(XtX, matrix(as(XWY, "matrix"), ncol = 1))$x
cf <- matrix(cf, nrow = ncol(X$X1))
- if (constr == 1) cf[1,] <- cf[1,] + switch(attr(X[[Xc]], "Ts_constraint"),
+ if (constr == 1) cf[1,] <- cf[1,] + switch(attr(X[[Xc]], "Ts_constraint"),
"increasing" = my,
"decreasing" = -my)
- if (constr == 2) cf[,1] <- cf[,1] + switch(attr(X[[Xc]], "Ts_constraint"),
+ if (constr == 2) cf[,1] <- cf[,1] + switch(attr(X[[Xc]], "Ts_constraint"),
"increasing" = my,
"decreasing" = -my)
cf
}
+
+differences_LSEI <- function(constraint, diff_mat = FALSE, ncol = NULL) {
+ if (length(constraint) == 1) {
+ diff <- switch(constraint,
+ positive = 0,
+ negative = 0,
+ increasing = 1,
+ decreasing = -1,
+ convex = 2,
+ concave = -2)
+ } else {
+ diff <- lapply(constraint, function(x)
+ switch(x,
+ positive = 0,
+ negative = 0,
+ increasing = 1,
+ decreasing = -1,
+ convex = 2,
+ concave = -2))
+ }
+ if (!diff_mat) {
+ return(diff)
+ } else {
+ if (diff != 0) {
+ D <- sign(diff) * diff(diag(ncol), differences = abs(diff))
+ } else {
+ D <- ifelse(constraint == "positive", 1, -1) * diag(ncol)
+ }
+ }
+}
+
+
+## least squares with equalities and inequalities
+solveLSEI <- function(XtX, Xty, D = NULL, constraint = "none") {
+ if (constraint == "none" && is.null(D))
+ return(solve(XtX, Xty, LINPACK = FALSE))
+
+ if (is.null(D)) {
+ ## i.e., if constraint != "none"
+ D <- differences_LSEI(constraint, diff_mat = TRUE, ncol = ncol(XtX))
+ }
+
+ cf <- lsei(A= XtX, B = Xty, G = D, H = rep(0, nrow(D)), type = 1,
+ E = matrix(0, ncol(D), ncol(D)), F = rep(0, ncol(D)),
+ tol = .Machine$double.eps,
+ tolrank = c(.Machine$double.eps, .Machine$double.eps),
+ fulloutput = TRUE)$X
+ cf
+}
More information about the Mboost-commits
mailing list