The R Project SVN R

## partly moved from ../man/smooth.spline.Rd , quite system-specific.
if(!dev.interactive(TRUE)) pdf("smooth.spline-test.pdf")

##-- artificial example
y18 <- c(1:3, 5, 4, 7:3, 2*(2:5), rep(10, 4))
## "truly 64 bit platform" {have seen "x86-64" instead of "x86_64")
(b.64 <- grepl("^x86.64", Sys.info()[["machine"]]) &&
     .Machine$sizeof.pointer > 4)## "truly 64 bit platform"
(Lb.64 <- b.64 && Sys.info()[["sysname"]] == "Linux" && .Machine$sizeof.pointer == 8)
## i386-Linux: Df ~= (even! > ) 18 : interpolating -- much smaller PRESS
## It is the too low 'low = -3' which "kills" the algo; low= -2.6 still ok
## On other platforms, e.g., x64, ends quite differently (and fine)
## typically with Df = 8.636
(s2. <- smooth.spline(y18, cv = TRUE,
                      control = list(trace=TRUE, tol = 1e-6,
                                     low = if(b.64) -3 else -2)))
plot(y18)
xx <- seq(1,length(y18), len=201)
lines(predict(s2., xx), col = 4)
mtext(deparse(s2.$call,200), side= 1, line= -1, cex= 0.8, col= 4)

(sdf8 <- smooth.spline(y18, df = 8, control=list(trace=TRUE)))# 11 iter.
sdf8$df - 8 # -0.0009159978
(sdf8. <- smooth.spline(y18, df = 8, control=list(tol = 1e-8)))# 14 iter.

## This gave error: "... spar 'way too large'" -- now sees in dpbfa() that it can't factorize
## --> and gives *warning* about too large spar only
## e <- try(smooth.spline(y18, spar = 50)) #>> error
## stopifnot(inherits(e, "try-error"))
ss50 <- try(smooth.spline(y18, spar = 50)) #>> warning only (in R >= 3.4.0) -- FIXME ??
   e <- try(smooth.spline(y18, spar = -9)) #>> error : .. too small', not on 32-bit
## if(Lb.64) stopifnot(inherits(e, "try-error"))
if(Lb.64) inherits(e, "try-error") else "not Linux 64-bit"
## I see (in 32 bit Windows),
b.64 || inherits(ss50, "try-error")  # TRUE .. always?

## "extreme" range of spar, i.e., 'lambda' directly  (" spar = c(lambda = *) "):
##  ---------------------  --> problem/bug for too large lambda
e10 <- c(-20, -10, -7, -4:4, 7, 10)
(lams <- setNames(10^e10, paste0("lambda = 10^", e10)))
lamExp <- as.expression(lapply(e10, function(E)
                substitute(lambda == 10^e, list(e = E))))
sspl <- lapply(lams, function(LAM) try(smooth.spline(y18, lambda = LAM)))
sspl
ok <- vapply(sspl, class, "") == "smooth.spline"
stopifnot(ok[e10 <= 7])
ssok <- sspl[ok]
ssGet  <- function(ch) t(sapply(ssok, `[` , ch))
ssGet1 <- function(ch)   sapply(ssok, `[[`, ch)
stopifnot(all.equal(ssGet1("crit"), ssGet1("cv.crit"), tol = 1e-10))# seeing rel.diff = 6.57e-12
## Interesting:  for really large lambda, solution "diverges" from the straight line
ssGet(c("lambda", "df", "crit", "pen.crit"))

plot(y18); lines(predict(s2., xx), lwd = 5, col = adjustcolor(4, 1/4))
invisible(lapply(seq_along(ssok), function(i) lines(predict(ssok[[i]], xx), col=i)))
i18 <- 1:18
abline(lm(y18 ~ i18), col = adjustcolor('tomato',1/2), lwd = 5, lty = 3)
## --> lambda = 10^10 is clearly wrong: a *line* but not the L.S. one
legend("topleft", lamExp[ok], ncol = 2, bty = "n", col = seq_along(ssok), lty=1)

##--- Explore 'all.knots' and 'keep.stuff'

s2   <- smooth.spline(y18, cv = TRUE, keep.stuff=TRUE)

s2.7  <- smooth.spline(y18, cv = TRUE, keep.stuff=TRUE, nknots = 7)
s2.11 <- smooth.spline(y18, cv = TRUE, keep.stuff=TRUE, nknots = 11)
plot(y18)
lines(predict(s2, xx), lwd = 5, col = adjustcolor(4, 1/4))
lines(predict(s2.7,  xx), lwd = 3, col = adjustcolor("red", 1/4))
lines(predict(s2.11, xx), lwd = 2, col = adjustcolor("forestgreen", 1/4))
## s2.11 is very close to 's2'

if(!requireNamespace("Matrix", quietly = TRUE) && !interactive())
    q("no")

if(Lb.64 && interactive()) ## extra checks (from above), but _not_ part of R checks
    stopifnot(inherits(e, "try-error"))
## in any case:
rbind("s-9_err" = inherits(e, "try-error"),
      "s+50_err"= inherits(ss50, "try-error"))


aux2Mat <- function(auxM) {
    stopifnot(is.list(auxM),
              identical(vapply(auxM, class, ""),
                        setNames(rep("numeric", 4), c("XWy", "XWX", "Sigma", "R"))))
    ## requireNamespace("Matrix")# want sparse matrices
    nk <- length(XWy <- auxM[["XWy"]])
    list(XWy = XWy,
         XWX =  Matrix::bandSparse(nk, k= 0:3, diagonals= matrix(auxM[[ "XWX" ]], nk,4), symmetric=TRUE),
         Sigma= Matrix::bandSparse(nk, k= 0:3, diagonals= matrix(auxM[["Sigma"]], nk,4), symmetric=TRUE))
}

## "Prove" basic property :
##
##     \hat{\beta} =  (X'W X + \lambda \Sigma)^{-1} X'W y
##     ---------------------------------------------------
##
chkB <- function(smspl, tol = 1e-10) {
    stopifnot(inherits(smspl, "smooth.spline"))
    if(!is.list(smspl$auxM))
        stop("need result of  smooth.spline(., keep.stuff = TRUE)")
    lM <- aux2Mat(smspl$auxM)
    beta.hat <- solve(lM$XWX + smspl$lambda * lM$Sigma, lM$XWy)
    all.equal(as.vector(beta.hat),
              smspl$fit$coef, tolerance = tol)
}

stopifnot(chkB(s2))
stopifnot(chkB(s2.7))
stopifnot(chkB(s2.11))

lM <- aux2Mat(s2$auxM)
A <- lM$XWX + s2$lambda * lM$Sigma
R <- Matrix::chol(A)
c. <- s2$fit$coef
stopifnot(all.equal(c., as.vector( solve(A, lM$XWy))) )

## c' Sigma c =
pen <- as.vector(c. %*% lM$Sigma %*% c.)
c(unscaled.penalty = pen,
  scaled.penalty   = s2$lambda * pen)

Sigma.tit <- quote(list(Sigma == Omega, "where"~~ Omega[list(j,k)] ==
                                      integral({B[j]*second}(t)~{B[k]*second}(t)~dt)))
Matrix::image(lM$XWX, main = quote({X*minute}*W*X))
Matrix::image(lM$Sigma, main = Sigma.tit)
Matrix::image(A, main = quote({X*minute}*W*X + lambda*Sigma))
Matrix::image(R, main = quote(R == chol({X*minute}*W*X + lambda*Sigma)))


## Specifying  'all.knots' ourselves

## 1) compatibly :
s2.7.k  <- smooth.spline(y18, cv = TRUE, keep.stuff=TRUE,
                         all.knots = s2.7$fit$knot[3+ 1:7])
ii <- names(s2.7) != "call"
stopifnot( all.equal(s2.7  [ii],
                     s2.7.k[ii]))

## 2) "free" but approximately in [0,1]
s2.9f  <- smooth.spline(y18, cv = TRUE, keep.stuff=TRUE,
                        all.knots = seq(0, 1, length.out = 9))
lines(predict(s2.9f, xx), lwd = 2, lty=3, col = adjustcolor("tomato", 1/2))
## knots partly outside [0,1]  --- is that correct ? (see below)
s2.7f  <- smooth.spline(y18, cv = TRUE, keep.stuff=TRUE,
                        all.knots = c(-1,1,3,5,7,9,12)/10)

if(FALSE) { ## not allowed (currently)
    ## knots partly *inside* [0,1] i.e. data outside knots
    s2.5f  <- smooth.spline(y18, cv = TRUE, keep.stuff=TRUE, control=list(trace=TRUE),
                            all.knots = c(1,3,5,7,9)/10)
    ## ------ OOOPS!  Segmentation fault ... "in attrib.c" {when returning from .Fortran()}
    lines(predict(s2.5f, xx), lwd = 2, lty=3, col = adjustcolor("brown", 1/2))
}
##' back-transform knots to "data-scale":
dScaledKnots <- function(smsp, drop.ends=TRUE) {
    stopifnot(inherits(smsp, "smooth.spline"))
    sf <- smsp$fit
    nk <- length(kk <- sf$knot)
    stopifnot((nk <- length(kk <- sf$knot)) >= 7)
    if(drop.ends) kk <- kk[4:(nk-3)]
    sf$min + sf$range * kk
}
pLines <- function(ss) {
    abline(v = dScaledKnots(ss), lty=3, col=adjustcolor("black", 1/2))
    abline(h = 0, v = range(ss$x), lty=4, lwd = 1.5, col="skyblue4")
}

## The following shows the data boundaries are used even when the knots are outside:
xe <- seq(-5, 25, length=256)
##
plot(y18, xlim=range(xe), ylim = c(-4,10)+.5, xlab="x")
lines(predict(s2.7f, xe), col=2, lwd = 2)
pLines(s2.7f)

str(m2 <- predict(s2.7f, x=xe, deriv=2)) # \hat{m''}(x)
plot(m2, type="l", col=2, lwd = 2,
     main = "m''(x) -- for  m(.) := smooth.spl(*, all.knots=c(..))",
     sub = "(knots shown as vertical dotted lines)")
pLines(s2.7f)

## same phenomenon (data boundaries, ...):
m1 <- predict(s2.7f, x=xe, deriv = 1) # \hat{m'}(x)
plot(m1, type="l", col=2, lwd = 2,
     main = "m'(x) -- for m(.) := smooth.spl(*, all.knots=c(..))",
     sub = "(knots shown as vertical dotted lines)")
pLines(s2.7f)