[Mboost-commits] r710 - pkg/mboostDevel/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Apr 20 18:00:39 CEST 2013
Author: thothorn
Date: 2013-04-20 18:00:39 +0200 (Sat, 20 Apr 2013)
New Revision: 710
Modified:
pkg/mboostDevel/R/bkronecker.R
pkg/mboostDevel/R/bl.R
pkg/mboostDevel/R/helpers.R
Log:
add alternative monotone B-splines; very experimental
Modified: pkg/mboostDevel/R/bkronecker.R
===================================================================
--- pkg/mboostDevel/R/bkronecker.R 2013-04-17 14:54:04 UTC (rev 709)
+++ pkg/mboostDevel/R/bkronecker.R 2013-04-20 16:00:39 UTC (rev 710)
@@ -34,6 +34,9 @@
dpp <- function(weights) {
+ if (!is.null(attr(X$X1, "deriv")) || !is.null(attr(X$X2, "deriv")))
+ stop("fitting of derivatives of B-splines not implemented")
+
W <- matrix(weights, nrow = n1, ncol = n2)
### X = kronecker(X2, X1)
@@ -60,12 +63,22 @@
}
XtX <- XtX + lambda * K
+ ### nnls
+ constr <- (!is.null(attr(X$X1, "constraint"))) +
+ (!is.null(attr(X$X2, "constraint")))
+
+ if (constr == 2)
+ stop("only one dimension may be subject to constraints")
+ constr <- constr > 0
+
## matrizes of class dgeMatrix are dense generic matrices; they should
## be coerced to class matrix and handled in the standard way
if (is(XtX, "Matrix") && !extends(class(XtX), "dgeMatrix")) {
XtXC <- Cholesky(forceSymmetric(XtX))
mysolve <- function(y) {
Y <- matrix(y, nrow = n1) * W
+ if (constr)
+ return(nnls2D(X, as(XtXC, "matrix"), Y))
XWY <- as.vector(crossprod(X$X1, Y) %*% X$X2)
solve(XtXC, XWY) ## special solve routine from
## package Matrix
@@ -77,6 +90,8 @@
}
mysolve <- function(y) {
Y <- matrix(y, nrow = n1) * W
+ if (constr)
+ return(nnls2D(X, as(XtX, "matrix"), Y))
XWY <- crossprod(X$X1, Y) %*% X$X2
solve(XtX, matrix(as(XWY, "matrix"), ncol = 1),
LINPACK = FALSE)
Modified: pkg/mboostDevel/R/bl.R
===================================================================
--- pkg/mboostDevel/R/bl.R 2013-04-17 14:54:04 UTC (rev 709)
+++ pkg/mboostDevel/R/bl.R 2013-04-20 16:00:39 UTC (rev 710)
@@ -167,7 +167,8 @@
### hyper parameters for P-splines baselearner (including tensor product P-splines)
hyper_bbs <- function(mf, vary, knots = 20, boundary.knots = NULL, degree = 3,
differences = 2, df = 4, lambda = NULL, center = FALSE,
- cyclic = FALSE) {
+ cyclic = FALSE, constraint = c("none", "increasing", "decreasing"),
+ deriv = 0L) {
knotf <- function(x, knots, boundary.knots) {
if (is.null(boundary.knots))
@@ -197,8 +198,13 @@
ret[[n]] <- knotf(mf[[n]], if (is.list(knots)) knots[[n]] else knots,
if (is.list(boundary.knots)) boundary.knots[[n]]
else boundary.knots)
+ constraint <- match.arg(constraint)
+ if (cyclic & constraint != "none")
+ stop("constraints not implemented for cyclic B-splines")
+ stopifnot(is.numeric(deriv) & length(deriv) == 1)
list(knots = ret, degree = degree, differences = differences,
- df = df, lambda = lambda, center = center, cyclic = cyclic)
+ df = df, lambda = lambda, center = center, cyclic = cyclic,
+ constraint = constraint, deriv = deriv)
}
### model.matrix for P-splines baselearner (including tensor product P-splines)
@@ -209,12 +215,15 @@
X <- bsplines(mf[[i]],
knots = args$knots[[i]]$knots,
boundary.knots = args$knots[[i]]$boundary.knots,
- degree = args$degree)
+ degree = args$degree,
+ constraint = args$constraint,
+ deriv = args$deriv)
if (args$cyclic) {
X <- cbs(mf[[i]],
knots = args$knots[[i]]$knots,
boundary.knots = args$knots[[i]]$boundary.knots,
- degree = args$degree)
+ degree = args$degree,
+ deriv = args$deriv)
}
class(X) <- "matrix"
return(X)
@@ -291,6 +300,10 @@
} else {
K <- crossprod(K)
}
+ if (!is.null(attr(X, "constraint"))) {
+ D <- attr(X, "D")
+ K <- crossprod(D, K) %*% D
+ }
}
if (length(mm) == 2) {
suppressMessages(
@@ -467,7 +480,8 @@
### P-spline (and tensor-product spline) baselearner
bbs <- function(..., by = NULL, index = NULL, knots = 20, boundary.knots = NULL,
degree = 3, differences = 2, df = 4, lambda = NULL, center = FALSE,
- cyclic = FALSE) {
+ cyclic = FALSE, constraint = c("none", "increasing", "decreasing"),
+ deriv = 0) {
if (!is.null(lambda)) df <- NULL
@@ -541,13 +555,14 @@
ret$dpp <- bl_lin(ret, Xfun = X_bbs,
args = hyper_bbs(mf, vary, knots = knots, boundary.knots =
boundary.knots, degree = degree, differences = differences,
- df = df, lambda = lambda, center = center, cyclic = cyclic))
+ df = df, lambda = lambda, center = center, cyclic = cyclic,
+ constraint = constraint, deriv = deriv))
return(ret)
}
### cyclic B-splines
### adapted version of mgcv:cSplineDes from S.N. Wood
-cbs <- function (x, knots, boundary.knots, degree = 3) {
+cbs <- function (x, knots, boundary.knots, degree = 3, deriv) {
# require(splines)
nx <- names(x)
x <- as.vector(x)
@@ -564,10 +579,11 @@
(boundary.knots[2] - knots[(nKnots - ord + 1):(nKnots - 1)]),
knots)
ind <- x > xc
- X <- splineDesign(knots, x, ord, outer.ok = TRUE)
+ X <- splineDesign(knots, x, ord, derivs = rep(deriv, length(x)), outer.ok = TRUE)
x[ind] <- x[ind] - boundary.knots[2] + boundary.knots[1]
if (sum(ind)) {
- Xtmp <- splineDesign(knots, x[ind], ord, outer.ok = TRUE)
+ Xtmp <- splineDesign(knots, x[ind], ord, derivs = rep(deriv, length(ind)),
+ outer.ok = TRUE)
X[ind, ] <- X[ind, ] + Xtmp
}
## handling of NAs
@@ -580,11 +596,13 @@
attr(X, "degree") <- degree
attr(X,"knots") <- knots
attr(X,"boundary.knots") <- boundary.knots
+ if (deriv != 0)
+ attr(X, "deriv") <- deriv
dimnames(X) <- list(nx, 1L:ncol(X))
return(X)
}
-bsplines <- function(x, knots, boundary.knots, degree){
+bsplines <- function(x, knots, boundary.knots, degree, constraint, deriv){
nx <- names(x)
x <- as.vector(x)
## handling of NAs
@@ -600,17 +618,30 @@
## complete knot mesh
k <- c(bk_lower, knots, bk_upper)
## construct design matrix
- X <- splineDesign(k, x, degree + 1, outer.ok = TRUE)
+ X <- splineDesign(k, x, degree + 1, derivs = rep(deriv, length(x)),
+ outer.ok = TRUE)
## handling of NAs
if (nas) {
tmp <- matrix(NA, length(nax), ncol(X))
tmp[!nax, ] <- X
X <- tmp
}
+ ### constraints; experimental
+ D <- diag(ncol(X))
+ D[lower.tri(D)] <- 1
+ X <- switch(constraint, "none" = X,
+ "increasing" = X %*% D,
+ "decreasing" = -X %*% D)
## add attributes
attr(X, "degree") <- degree
- attr(X,"knots") <- knots
- attr(X,"boundary.knots") <- list(lower = bk_lower, upper = bk_upper)
+ attr(X, "knots") <- knots
+ attr(X, "boundary.knots") <- list(lower = bk_lower, upper = bk_upper)
+ if (constraint != "none")
+ attr(X, "constraint") <- constraint
+ if (constraint != "none")
+ attr(X, "D") <- D
+ if (deriv != 0)
+ attr(X, "deriv") <- deriv
dimnames(X) <- list(nx, 1L:ncol(X))
return(X)
}
@@ -639,6 +670,9 @@
dpp <- function(weights) {
+ if (!is.null(attr(X, "deriv")))
+ stop("fitting of derivatives of B-splines not implemented")
+
weights[!Complete.cases(mf)] <- 0
w <- weights
if (!is.null(index))
@@ -653,17 +687,25 @@
## be coerced to class matrix and handled in the standard way
if (is(X, "Matrix") && !extends(class(XtX), "dgeMatrix")) {
XtXC <- Cholesky(forceSymmetric(XtX))
- mysolve <- function(y)
- solve(XtXC, crossprod(X, y)) ## special solve routine from
- ## package Matrix
+ mysolve <- function(y) {
+ if (is.null(attr(X, "constraint")))
+ return(solve(XtXC, crossprod(X, y))) ## special solve routine from
+ ## package Matrix
+ ### non-negative LS only at the moment
+ return(nnls1D(as(XtX, "matrix"), as(X, "matrix"), y))
+ }
} else {
if (is(X, "Matrix")) {
## coerce Matrix to matrix
X <- as(X, "matrix")
XtX <- as(XtX, "matrix")
}
- mysolve <- function(y)
- solve(XtX, crossprod(X, y), LINPACK = FALSE)
+ mysolve <- function(y) {
+ if (is.null(attr(X, "constraint")))
+ return(solve(XtX, crossprod(X, y), LINPACK = FALSE))
+ ### non-negative LS only at the moment
+ return(nnls1D(XtX, X, y))
+ }
}
fit <- function(y) {
Modified: pkg/mboostDevel/R/helpers.R
===================================================================
--- pkg/mboostDevel/R/helpers.R 2013-04-17 14:54:04 UTC (rev 709)
+++ pkg/mboostDevel/R/helpers.R 2013-04-20 16:00:39 UTC (rev 710)
@@ -121,3 +121,42 @@
}
list(hatmatrix = B, trace = tr)
}
+
+nnls1D <- function(XtX, X, y) {
+
+ my <- min(y)
+ if (my < 0) {
+ ### first column is intercept
+ stopifnot(max(abs(X[,1] - 1)) < sqrt(.Machine$double.eps))
+ y <- y - my
+ }
+ stopifnot(require("nnls"))
+ cf <- nnls(XtX, crossprod(X, y))$x
+ if (my < 0) cf[1] <- cf[1] + my
+ cf
+}
+
+nnls2D <- function(X, XtX, Y) {
+
+ constr <- which(c(!is.null(attr(X$X1, "constraint")),
+ !is.null(attr(X$X2, "constraint"))))
+ if (length(constr) == 2)
+ stop("only one dimension may be subject to constraints")
+ my <- min(Y)
+ if (my < 0) {
+ ### first column is intercept
+ stopifnot(max(abs(X[[paste("X", constr, sep = "")]][,1] - 1)) <
+ sqrt(.Machine$double.eps))
+ Y <- Y - my
+ }
+ XWY <- as.vector(crossprod(X$X1, Y) %*% X$X2)
+
+ stopifnot(require("nnls"))
+ cf <- nnls(XtX, matrix(as(XWY, "matrix"), ncol = 1))$x
+ cf <- matrix(cf, nrow = ncol(X$X1))
+ if (my < 0) {
+ if (constr == 1) cf[1,] <- cf[1,] + my
+ if (constr == 2) cf[,1] <- cf[,1] + my
+ }
+ cf
+}
More information about the Mboost-commits
mailing list