| require("splines") |
| |
| ## Bug report PR#16549 - 'bad value from splineDesign' |
| ## Date: Wed, 30 Sep 2015 12:12:47 +0000 |
| ## https://bugs.r-project.org/bugzilla/show_bug.cgi?id=16549 |
| |
| ## Reporter: roconnor@health.usf.edu (extended original example code) |
| |
| kn8.01 <- c(0,0,0,0,1,1,1,1) |
| m1 <- rbind(c(1, 0,0,0)) |
| m.1 <- rbind(c(0,0,0, 1)) |
| stopifnot( |
| ## the first gave (0 0 0 0) instead of m1 : |
| all.equal(splineDesign(kn8.01, c(0,2), outer.ok = TRUE), rbind(m1, 0)), |
| all.equal(splineDesign(kn8.01, c( 0,1,2), outer.ok = TRUE), rbind(m1, m.1, 0)), |
| all.equal(splineDesign(kn8.01, c(-1,0,1), outer.ok = TRUE), rbind(0, m1, m.1)), |
| all.equal(splineDesign(kn8.01, 0, outer.ok = TRUE), m1), |
| all.equal(splineDesign(kn8.01, 0), m1), |
| TRUE) |
| |
| ## The original fix proposal introduced a new bug, visible here: |
| S <- splineDesign(c(-3, -3, -2, 0, 2, 3, 3), x= -3:3, outer.ok=TRUE) |
| ## (had a NaN in the lower-right corner) |
| stopifnot(all.equal(S, |
| rbind(0, c(22,3,0)/45, c(193, 6*23, 9)/360, |
| c(1,3,1)/5, |
| c(9, 6*23, 193)/360, |
| c(0,3,22)/45, 0), |
| tolerance = 1e-14)) |
| |
| chkSum.Ok <- TRUE ## for check |
| chkSum <- function(knots, n = 1 + 2^9, ord = 4) { |
| stopifnot(is.numeric(knots), !is.unsorted(knots), (n.k <- length(knots)) >= ord) |
| dk <- diff(rk <- range(uk <- unique(knots))) |
| if(dk == 0) dk <- uk[1]/4 |
| d.x <- dk / (2*n.k) |
| ## x: the unique knots |
| x <- sort(c(uk, seq.int(min(knots)-d.x, max(knots)+d.x, length.out = n))) |
| bb <- splineDesign(knots, x = x, ord = ord, outer.ok = TRUE) |
| is.x.in <- knots[ord] <= x & x <= knots[n.k-(ord-1)] |
| ## ~~~~ ~~~~ same as in splineDesign(*, outer.ok=TRUE) |
| ## |
| ### no longer needed: |
| ## work around "infelicity" (or not; not sure if to call "bug") |
| ## n.RHk := #{duplicated RHS boundary knots} |
| ## n.RHk <- match(FALSE, rev(duplicated(knots))) |
| ## if(ord > 1 && n.RHk > ord) { ## the (ord == 1) case "works" via spike = 1 |
| ## ## <==> knots[n.k -j ] == knots[n.k] for j = 0,1,..,(n.RHk-1) |
| ## ## then spl(x= RHS-knot, knots) == 0, even though "should" be 1 |
| ## ## "FIX it up" for here. FIXME: do in splineDesign() or it's C code ?! |
| ## iR <- which(x == knots[n.k]) |
| ## ## ncol(bb) == n.k - ord |
| ## j <- n.k - n.RHk # == ncol(bb) - (n.RHk - ord) |
| ## if(any(bb[iR, j] == 0)) bb[iR, j] <- 1 |
| ## } |
| sumB <- rowSums(bb) |
| if(any(iBad <- !is.finite(sumB))) { |
| chkSum.Ok <<- FALSE ## for check |
| cat("** _FIXME_ NON-finite values in sumB: ord = ", ord, "; |knots| =", n.k,"\n") |
| cat("knots <- "); dput(knots) |
| cat("non-finite at x = "); dput(x[iBad]) |
| } else if(length(bb)) { # only when bb[] is not 0-dimensional |
| eps <- 3*.Machine$double.eps |
| stopifnot(abs(1 - sumB[is.x.in]) <= 2*eps, 0 <= sumB+eps, sumB-2*eps <= 1) |
| ## TODO: now also check derivatives |
| } |
| invisible(bb) |
| } |
| |
| plotSplD <- function(knots, n = 2^10, ord = 4, type = "l", |
| ylim = c(0,1), ylab = "B-splines", ...) { |
| stopifnot(is.numeric(knots), !is.unsorted(knots), (n.k <- length(knots)) >= ord) |
| dk <- diff(range(uk <- unique(knots))) |
| if(dk == 0) dk <- uk[1]/4 |
| d.x <- dk / (2*n.k) |
| ## x: will contain the unique knots uk |
| x <- sort(c(uk, seq.int(min(knots)-d.x, max(knots)+d.x, length.out = n))) |
| bb <- splineDesign(knots, x = x, ord = ord, outer.ok = TRUE) |
| matplot(x, bb, type=type, ylim=ylim, ylab=ylab, ...) |
| sumB <- rowSums(bb) |
| abline(v = knots, lty = 3, col = "light gray") |
| abline(v = knots[c(ord, n.k-(ord-1))], lty = 3, col = "gray10") |
| lines(x, sumB, col = adjustcolor("red", 0.4), lwd = 3) |
| abline(h=1, lty=2, col = adjustcolor(1, 0.4), lwd = 2) |
| ## lty, col,: from matplot() |
| legend(mean(x), 0.98, legend = paste0("B_", 1:ncol(bb)), lty=1:5, col=1:6, |
| bty = "n", xjust = 0.5) |
| invisible(list(x=x, splineDesign = bb)) |
| } |
| |
| if(!dev.interactive(orNone=TRUE)) pdf("spline-tst.pdf") |
| |
| plotSplD(kn8.01) |
| chkSum (kn8.01) |
| |
| ## from ../man/splineDesign.Rd : |
| knots <- c(1,1.8,3:5,6.5,7,8.1,9.2,10) # 10 => 10-4 = 6 Basis splines |
| str(plotSplD(knots)) # cubic splines, adding to 1 in [ 4, 7 ] |
| str(plotSplD(knots, ord=3))# quadratic splines, adding to 1 in [ 3, 8.1] |
| str(plotSplD(knots, ord=2))# linear splines, adding to 1 in [1.8,9.2] |
| str(plotSplD(knots, ord=1))# constant splines, adding to 1 in [1, 10] |
| |
| chkSum(knots) |
| chkSum(knots, ord=3) |
| chkSum(knots, ord=2) |
| chkSum(knots, ord=1) |
| |
| ## cases that failed too {Linux lynne 4.1.6 x86_64} |
| chkSum(c(1:5, 9,9), ord=2) # ok for ord in {1,3,4} |
| chkSum(c(1:4, 9,9,9), ord=3) |
| chkSum(c(1:3, 9,9,9,9), ord=4) |
| ## These failed in R <= 3.2.2 (but not after first fix): |
| chkSum(c(0,0,0,0, 1:3), ord=4) |
| chkSum(c(0,0,0, 1:4), ord=3) |
| chkSum(c(0,0, 1:5), ord=2) |
| |
| ## This should be symmetric, but was not; |
| ## the graphic is maybe most convincing: |
| k6.01 <- c(0,0,0, 1,1,1) |
| round(with(plotSplD(k6.01, ord=2, n=16), cbind(x, splineDesign)), 4) |
| x8 <- (0:8)/8 |
| sp8 <- splineDesign(k6.01, x=x8, ord=2) |
| print.table(8*sp8, zero.print=".") |
| stopifnot(all.equal(8*sp8, |
| cbind(0, 8:0, 0:8, 0), tol = 1e-14)) |
| ## or just |
| splineDesign(k6.01, x=0:1, ord=2) |
| ## 0 1 0 0 |
| ## 0 0 1 0 --- [finally !] |
| stopifnot(identical(cbind(0, diag(2), 0), |
| splineDesign(k6.01, x=0:1, ord=2))) |
| |
| ## Further: |
| for(k in 0:8) |
| chkSum(c(0:k, 9,9,9), ord=2) |
| for(k in 0:8) |
| chkSum(c(0:k, 9,9,9,9), ord=3) |
| for(k in 0:8) |
| chkSum(c(0:k, 9,9,9,9,9), ord=4) |
| |
| ## look at two examples |
| r <- plotSplD(c(0:4, 8.9,9,9,9), ord=3)# B_6 is narrow spike |
| r. <- plotSplD(c(0:4, 9, 9,9,9), ord=3)# B_6 diverged to all zero. |
| if(interactive()) |
| round(with(r., cbind(x, splineDesign)), 4) |
| stopifnot(r.$splineDesign[, 6] == 0) |
| ## one could argue that B_6 should also have finite area, and hence be "a delta". |
| ## ==> sum_j B_j(x = 9) would then be Inf or undefined .. |
| ## more sensible: B_6 should be omitted and should have B_5(9) == 1 |
| ## now implemented: |
| stopifnot(identical(c(0, 0, 0, 0, 1, 0), |
| with(r., splineDesign[x == 9,]))) |
| |
| ## Everything is fine, if left-right mirrored / reversed : |
| r03 <- plotSplD(c(0,0,0,0, 1:5), ord=3) |
| ## here, B_1 is "delta" and should rather be omitted |
| stopifnot(identical(c(0, 1, 0, 0, 0, 0), |
| with(r03, splineDesign[x == 0,]))) |
| ## 0 1 0 0 0 0 -- as it should be .. |
| round(with(r03, cbind(x, splineDesign)[50:60,]), 4) |
| |
| r04 <- plotSplD(c(0,0,0,0, 1:5), ord=4) |
| ## here, B_1 is "delta" and should rather be omitted |
| stopifnot(identical(c(1, 0, 0, 0, 0), |
| with(r04, splineDesign[x == 0,]))) |
| round(with(r04, cbind(x, splineDesign)[50:60,]), 4) |
| |
| r0.4 <- plotSplD(c(0,0,0, 1:5), ord=4) # basis is nice and correct |
| |
| |
| set.seed(17) |
| |
| for(n in 1:1000) { |
| if(n %% 50 == 0) cat(sprintf("n = %4d\n",n)) |
| kn <- sort.int(round(10* rnorm(4 + rpois(1, lambda=4)))) |
| for(oo in 1:4) |
| ## if(inherits(r <- tryCatch(chkSum(kn, ord=oo), error=identity), "error")) |
| ## cat(r$message, "\n chkSum(", deparse(kn), ", ord = ", oo, ")\n") |
| chkSum(kn, ord = oo) |
| } |
| |
| ## One of the cases with NaN {when used ( . <= x & x <= . )}: |
| |
| bb <- chkSum(c(-14, -4, 3, 5, 6, 15, 15)) |
| stopifnot(is.finite(rowSums(bb))) |
| |
| |
| stopifnot(chkSum.Ok) |
| |
| proc.time() |
| |
| ###---- "Bug report" (to R-core) from Trevor Hastie --- |
| ###----> bs(*, Boundary.knots = .) |
| ### needing boundary ajustment for correct extrapolation |
| |
| ## Trevor's Example, slightly modified and extended: |
| x <- seq(1.5, 8.5, by = 1/4) |
| set.seed(13) |
| y <- x + .01*(x - 5)^3 + rnorm(x) |
| |
| fit0 <- lm(y ~ bs(x, degree=3, knots=4)) |
| fit0.<- lm(y ~ bs(x, degree=3, knots=4, Boundary.knots=c(1,9)))# *NOT* outside |
| fit1 <- lm(y ~ bs(x, degree=3, knots=4, Boundary.knots=c(1,8)))# warning |
| fit2 <- lm(y ~ bs(x, degree=3, knots=4, Boundary.knots=c(2,8)))# warning "2 x" |
| |
| jx <- seq(from=-2,to=12, by=0.1) |
| p0 <- predict(fit0, list(x=jx)) |
| p0.<- predict(fit0, list(x=jx)) |
| p1 <- predict(fit1, list(x=jx)) |
| p2 <- predict(fit2, list(x=jx)) |
| stopifnot(all.equal(p0, p0.,tol=1e-14), |
| all.equal(p0, p1, tol=1e-14), |
| all.equal(p0, p2, tol=1e-14)) |
| ## ^^ p1 and p2 differed from p0 in R <= 3.2.2 |
| ## See numerical fuzz: |
| all.equal(p0, p0., tol=0) |
| all.equal(p0, p1, tol=0) |
| all.equal(p0, p2, tol=0) |
| all.equal(p1, p2, tol=0)# interestingly almost the same |
| |
| ## formula ==> print for default method |
| ispl <- with(women, interpSpline( height, weight )) |
| stopifnot(identical(format(formula(ispl)), |
| "weight ~ height")) ## was wrongly .Primitive(\"~\")(wei... |
| |
| ###------------- Problems for small n ------- want at least good error messages --- |
| |
| ## This gives an error, but not a "human readable" one for k <= 3 |
| ## -- now done: the English error message is "must have at least 'ord'=4 points" |
| ## FIXME: Would like to get degree=0 (constant), 1, 2 interpolation spline here |
| ## (and could use degree = 0 for both n=0 and n=1) ==> would have *no* error here |
| for(n in 0:4) { |
| cat("n = ", n,": ") |
| x <- cumsum(pmax(1, round(10*runif(n)))) # pmax(1,*): x[] must be distinct |
| y <- (x - 2)^2 + round(8*rnorm(n))/4 |
| if(!inherits(sp <- try(interpSpline(x,y)), "try-error")) print(sp) |
| cat("-------------------------------\n") |
| } |
| ## for n=4: |
| stopifnot(inherits(sp, "polySpline"), is.matrix(sp$coefficients), |
| identical(dim(sp$coefficients), c(4L, 4L))) |
| |
| ## This also gives a "hard to read" error _FIXME_ |
| try(ns(1[0])) # n = 0 |
| ## Error in splineDesign(Aknots, x, ord) : |
| ## length of 'derivs' is larger than length of 'x' -- barely ok + 2 warnings |
| try(bs(1[0])) # n = 0 : same error etc as ns() |
| |
| (b1 <- bs(pi)) # n = 1 : works fine |
| (n1 <- ns(pi)) # n = 1 : now ok; gave error: "qr.default(.. NA ...)" |
| |
| ##' keep {dim, dimnames} but nothing else : |
| noAttr <- function(x) `mostattributes<-`(x, list(dim=dim(x), dimnames=dimnames(x))) |
| stopifnot( |
| identical(noAttr(b1), rbind(c(`1`=0, `2`=0, `3`=0))), |
| all.equal(noAttr(n1), matrix(0.400891862868637, 1, dimnames=list(NULL,"1")), |
| tol = 5e-15)) |
| |
| d1 <- data.frame(u = 2, Y = 5) ## data set with 1 observation |
| summary(mbs <- lm(Y ~ bs(u), data=d1)) # fine though many coef etc are NA |
| summary(mbs1 <- lm(Y ~ bs(u, degree=1), data=d1)) # ok |
| stopifnot( |
| identical(coef(mbs), setNames(c(5, NA, NA, NA), |
| c("(Intercept)", paste0("bs(u)", 1:3)))) |
| , |
| identical(coef(mbs1), c("(Intercept)" = 5, "bs(u, degree = 1)" = NA))) |
| if(FALSE) |
| summary(mbs0 <- lm(Y ~ bs(u, degree=0), data=d1)) # error: degree >= 1 |
| |
| ## ns() has no 'degree' argument! |
| summary(mns <- lm(Y ~ ns(u), data=d1)) ## gave Error; now ok |
| stopifnot(identical(coef(mns), c("(Intercept)" = 5, "ns(u)" = NA)) |
| ## perfect prediction: all residuals == 0 : |
| , identical(residuals(mns), c(`1` = 0)) |
| , identical(residuals(mbs), c(`1` = 0)) |
| , identical(residuals(mbs1),c(`1` = 0)) |
| ) |