Nothing
R/gps.R
In gps: General P-Splines
Csparse2LTB <- function (A) {
if (!inherits(A, "dMatrix") || !inherits(A, "CsparseMatrix")) {
stop("'A' is not a \"dCsparseMatrix\"!")
}
n <- A@Dim[1L]
k <- A@Dim[2L]
i <- A@i
j <- rep.int(1:k, diff.default(A@p))
b1 <- sum(i == 0L)
b <- b1 - 1L
if (k < n || k > n + b) {
stop(sprintf("nrow(A) <= ncol(A) <= nrow(A) + %d expected!", b))
}
l <- j - i
if (any(l > b1 | l < 1L)) {
stop("Nonzero entries detected outside the indended band region!")
}
N1 <- sum(seq.int(n, by = -1L, length.out = b1))
N2 <- sum(seq.int(b, by = -1L, length.out = k - n))
nnz.max <- N1 + N2
x <- A@x
nnz <- length(x)
if (nnz == nnz.max) {
.Call("C_Csparse2LTB", b1, n, k, x, PACKAGE = "gps")
} else {
At <- matrix(0, b1, k)
At[cbind(l, i + 1L)] <- A@x
At
}
}
LTB2Csparse <- function (L, k = n, symmetric = FALSE) {
b1 <- nrow(L); n <- ncol(L); b <- b1 - 1L
k <- as.integer(k)
if (symmetric) {
k <- n
} else if (k < n || k > n + b) {
stop(sprintf("%d <= k <= %d expected!", n, n + b))
}
r <- n + b - k
p <- rep.int(seq.int(b1, by = -1L, length.out = r + 1L), c(n - r, rep.int(1, r)))
i <- sequence.default(nvec = p, from = seq.int(0L, n - 1L))
p <- c(0L, cumsum(p))
if (r) {
x <- numeric(p[length(p)])
x[seq_len(b1 * (k - b))] <- L[, seq_len(k - b)]
columns <- seq.int(to = n, length.out = r)
for (j in columns) {
ind <- seq.int(p[j] + 1L, p[j + 1L])
x[ind] <- L[seq_len(length(ind)), j]
}
} else {
x <- c(L)
}
if (symmetric) {
methods::new("dsCMatrix", i = i, p = p, x = x, Dim = c(k, n), uplo = "L")
} else {
methods::new("dgCMatrix", i = i, p = p, x = x, Dim = c(k, n))
}
}
LTB2Dense <- function (L, k) {
.Call("C_LTB2Dense", L, k, PACKAGE = "gps")
}
SolveLTB <- function (transA = FALSE, A, b, overwrite = FALSE) {
.Call("C_SolveLTB", transA, A, b, overwrite, PACKAGE = "gps")
}
LPBTRF <- function (A, overwrite = FALSE) {
.Call("C_LPBTRF", A, overwrite, PACKAGE = "gps")
}
DDt <- function (Dt) {
.Call("C_DDt", Dt, PACKAGE = "gps")
}
DtD <- function (Dt) {
.Call("C_DtD", Dt, PACKAGE = "gps")
}
Dt2ThinQR <- function (D, form.Q = TRUE) {
Dt <- Csparse2LTB(D)
A <- DDt(Dt)
Rt <- LPBTRF(A, overwrite = TRUE)
if (form.Q) {
denseD <- as_matrix(D)
Qt <- SolveLTB(transA = FALSE, Rt, denseD, overwrite = TRUE)
} else {
Qt <- NULL
}
list(Qt = Qt, Rt = Rt)
}
IsMonoInc <- function (x, n = length(x), xi = 1L) {
.Call("C_IsMonoInc", x, n, xi, PACKAGE = "gps") > 0L
}
MonoKnots <- function (xt, d) {
xt <- as.double(xt)
K <- length(xt)
if (K < 2 * d) stop("length(xt) >= 2 * d required!", call. = FALSE)
flag <- IsMonoInc(xt, n = K - 2 * (d - 1), xi = d)
if (!flag) stop("Domain knots are not strictly ascending!", call. = FALSE)
lbnd <- xt[seq.int(from = 1L, length.out = d)]
rbnd <- xt[seq.int(to = K, length.out = d)]
if (is.unsorted(lbnd) || is.unsorted(rbnd)) {
stop("Boundary knots are not non-decreasing!", call. = FALSE)
}
xt
}
PlaceKnots <- function (x, d, k, domain = NULL, uniform = FALSE, periodic = FALSE) {
xu <- sort.int(unique.default(x))
nx <- length(xu)
if (is.null(domain)) {
a <- xu[1L]
b <- xu[nx]
domain <- c(a, b)
} else {
a <- domain[1L]
b <- domain[2L]
if (xu[1L] < a || xu[nx] > b) {
stop("'domain' does not contain all x-values!")
}
}
degree <- d - 1
if (uniform) {
xd <- seq.int(a, b, length.out = k + 2)
h <- xd[2L] - xd[1L]
laux <- seq.int(to = a - h, by = h, length.out = degree)
raux <- seq.int(from = b + h, by = h, length.out = degree)
} else {
prob <- seq.int(0, 1, length.out = k + 2)
xd <- quantile(xu, prob, names = FALSE)
xd[c(1, k + 2)] <- domain
laux <- rep.int(a, degree)
raux <- rep.int(b, degree)
}
if (periodic) xd else c(laux, xd, raux)
}
MakeGrid <- function (xd, n, rm.dup = FALSE) {
xd <- as.double(xd)
if (!IsMonoInc(xd)) stop("'xd' is not strictly ascending!")
if (n == 1) {
lp <- xd[-length(xd)]
rp <- xd[-1]
mp <- 0.5 * (lp + rp)
return(mp)
}
if (rm.dup && (n == 2)) return(xd)
.Call("C_MakeGrid", xd, n, rm.dup, PACKAGE = "gps")
}
Zero2NA <- function (Bsparse) {
if (!inherits(Bsparse, "dgCMatrix")) {
stop("'Bsparse' is not a \"dgCMatrix\"!")
}
B <- matrix(NA_real_, Bsparse@Dim[1], Bsparse@Dim[2])
i <- Bsparse@i + 1L
j <- rep.int(seq_len(Bsparse@Dim[2]), diff.default(Bsparse@p))
B[cbind(i, j)] <- Bsparse@x
B
}
ExBS <- function (d = 4, uniform = TRUE, clamped = FALSE) {
if (d < 2 || d > 4) stop("d = 2, 3, 4 required!", call. = FALSE)
if (uniform) {
lbnd <- as.double(seq.int(to = -1, length.out = 4, by = 2))
rbnd <- as.double(seq.int(from = 1, length.out = 4, by = 2))
clamped <- FALSE
} else if (clamped) {
lbnd <- rep.int(-1, 4)
rbnd <- c(1, 2, 2.5, 4)
} else {
lbnd <- c(-3, -2.25, -1.75, -1)
rbnd <- c(1, 2, 2.5, 4)
}
xt <- c(lbnd, rbnd)
if (clamped) {
x <- MakeGrid(xt[4:8], n = 21)
} else {
x <- MakeGrid(xt, n = 21)
}
Bsparse <- splines::splineDesign(xt, x, d, outer.ok = TRUE, sparse = TRUE)
B <- Zero2NA(Bsparse)
list(xt = xt, clamped = clamped, d = d, x = x, B = B)
}
PlotExBS <- function (object) {
xt <- object$xt
d <- object$d
x <- object$x
B <- object$B
p <- ncol(B)
if (object$clamped) {
ip <- c(rep.int(NA_real_, 4 - d), apply(B[, (4 - d + 1):p], 2, which.max))
} else {
ip <- apply(B, 2, which.max)
}
xp <- x[ip]
Bp <- B[cbind(ip, 1:p)]
ymax <- max(Bp, na.rm = TRUE) * 1.2
ylim <- c(0, ymax)
nonzero <- seq.int(5 - d, length.out = d)
col.Bspl <- rep.int(8, p); col.Bspl[nonzero] <- seq_len(d)
op <- par(xpd = TRUE)
on.exit(par(op))
matplot(x, B, type = "l", lty = 1, lwd = 2, col = col.Bspl, ylim = ylim,
xaxs = "i", yaxs = "i", xaxt = "n", yaxt = "n", bty = "n", ann = FALSE)
segments(xt[1], 0, xt[8], 0)
polygon(c(-1, 1, 1, -1), c(0, 0, ymax, ymax), col = gray(0.3, 0.2), border = NA)
col.knots <- rep.int(8, 8); col.knots[seq.int(5 - d, 4 + d)] <- 1
if (object$clamped) {
stepsize <- 0.2
ystack <- seq.int(to = stepsize, length.out = 3, by = -stepsize)
points.default(xt[4:8], numeric(5), pch = 19, col = col.knots[4:8], cex = 1.5)
points.default(xt[1:3], ystack, pch = 19, col = col.knots[1:3], cex = 1.5)
} else {
points.default(xt, numeric(8), pch = 19, col = col.knots, cex = 1.5)
}
lab.Bspl <- vector("expression", p)
for (j in 1:p) {
lab.Bspl[j] <- eval(substitute(expression(B[list(j, d)]), list(j = j, d = d)))
}
text.default(xp, Bp, labels = lab.Bspl, col = col.Bspl, adj = c(0.2, 0), cex = 2)
lab.knots <- vector("expression", 8)
for (j in 1:8) {
lab.knots[j] <- eval(substitute(expression(t[j]), list(j = j)))
}
if (object$clamped) {
text.default(xt[1:3], ystack, labels = lab.knots[1:3], col = col.knots[1:3],
cex = 2, pos = 4)
text.default(xt[4:8], numeric(5), labels = lab.knots[4:8], col = col.knots[4:8],
cex = 2, pos = 1)
} else {
text.default(xt, numeric(8), labels = lab.knots, col = col.knots,
cex = 2, pos = 1)
}
usr <- par("usr")
text.default(usr[2], usr[4], labels = c(paste0("order: ", d)),
cex = 2, adj = c(1, 1))
}
DemoBS <- function (uniform = TRUE, clamped = FALSE) {
if (uniform && clamped) {
warning("uniform = TRUE implies clamped = FALSE!")
}
spl3 <- ExBS(d = 4, uniform, clamped)
spl2 <- ExBS(d = 3, uniform, clamped)
spl1 <- ExBS(d = 2, uniform, clamped)
if (uniform) {
caption <- "uniform B-splines"
} else if (clamped) {
caption <- "non-uniform & clamped B-splines"
} else {
caption <- "non-uniform B-splines"
}
op <- par(mfrow = c(3, 1), mar = c(2, 0.75, 0, 0.75), oma = c(0, 0, 2, 0))
on.exit(par(op))
PlotExBS(spl3)
PlotExBS(spl2)
PlotExBS(spl1)
title(caption, outer = TRUE, cex.main = 1.5)
}
DemoKnots <- function (aligned = TRUE) {
xt1 <- seq.int(0, 1, length.out = 12)
xt2 <- c(0, 0.08, 0.2, 0.27, 0.31, 0.42, 0.54, 0.59, 0.69, 0.82, 0.85, 1)
xt3 <- rep.int(xt2[4:9], c(4, 1, 1, 1, 1, 4))
xg1 <- MakeGrid(xt1, n = 21)
xg2 <- MakeGrid(xt2, n = 21)
xg3 <- MakeGrid(xt3[4:9], n = 21)
B1sparse <- splines::splineDesign(xt1, xg1, outer.ok = TRUE, sparse = TRUE)
B2sparse <- splines::splineDesign(xt2, xg2, outer.ok = TRUE, sparse = TRUE)
B3sparse <- splines::splineDesign(xt3, xg3, sparse = TRUE)
B1max <- max(B1sparse@x)
B2max <- max(B2sparse@x)
B3max <- max(B3sparse@x)
B1 <- Zero2NA(B1sparse)
B2 <- Zero2NA(B2sparse)
B3 <- Zero2NA(B3sparse)
Bspl.col <- c(1:6, 8, 7)
op <- par(mfrow = c(3, 1), mar = c(0.5, 0.5, 1.5, 0.5), xpd = TRUE)
on.exit(par(op))
matplot(xg1, B1, ann = FALSE, axes = FALSE, type = "l", col = Bspl.col, bty = "n",
lty = 1, lwd = 2, xlim = c(0, 1), xaxs = "i", yaxs = "i")
abline(h = 0, col = 8)
points.default(xt1, numeric(12), pch = 19, cex = 1.5)
polygon(xt1[c(4, 9, 9, 4)], c(0, 0, B1max, B1max), col = gray(0.3, 0.2), border = NA)
title("(a) uniform cubic B-splines", cex.main = 1.5)
matplot(xg2, B2, ann = FALSE, axes = FALSE, type = "l", col = Bspl.col, bty = "n",
lty = 1, lwd = 2, xlim = c(0, 1), xaxs = "i", yaxs = "i")
abline(h = 0, col = 8)
points.default(xt2, numeric(12), pch = 19, cex = 1.5)
polygon(xt2[c(4, 9, 9, 4)], c(0, 0, B2max, B2max), col = gray(0.3, 0.2), border = NA)
title("(b) non-uniform cubic B-splines", cex.main = 1.5)
xlim3 <- if (aligned) c(0, 1) else xt3[c(4, 9)]
matplot(xg3, B3, ann = FALSE, axes = FALSE, type = "l", col = Bspl.col, bty = "n",
lty = 1, lwd = 2, xlim = xlim3, xaxs = "i", yaxs = "i")
abline(h = 0, col = 8)
if (aligned) {
segments(0, 0, xt3[4], 0, col = Bspl.col[1], lwd = 2)
segments(xt3[9], 0, 1, 0, col = Bspl.col[8], lwd = 2)
segments(xt3[4], 0, xt3[4], 1, col = Bspl.col[1], lty = 2, lwd = 2)
segments(xt3[9], 0, xt3[9], 1, col = Bspl.col[8], lty = 2, lwd = 2)
}
ystack <- seq.int(0, by = 0.1, length.out = 4)
points.default(xt3, c(ystack, numeric(4), ystack), pch = 19, cex = 1.5)
polygon(xt3[c(4, 9, 9, 4)], c(0, 0, B3max, B3max), col = gray(0.3, 0.2), border = NA)
title("(c) non-uniform & clamped cubic B-splines", cex.main = 1.5)
}
PlotNull <- function (x, y, k, shape1 = 3, shape2 = 3, gps = FALSE) {
d <- 4
m <- 2
if (shape1 < 1 && shape2 < 1) {
xd <- seq.int(-1/8, 1/8, length.out = k + 2)
xd <- sign(xd) * abs(xd) ^ (1 / 3) + 0.5
distr <- "U-quadratic(0, 1)"
} else {
xd <- qbeta(seq.int(0, 1, length.out = k + 2), shape1, shape2)
distr <- sprintf("Beta(%d, %d)", shape1, shape2)
}
xt <- c(numeric(d - 1), xd, rep.int(1, d - 1))
K <- length(xt)
xg <- seq.int(0, 1, length.out = 101)
B <- splines::splineDesign(xt, x, d, sparse = TRUE)
Bg <- splines::splineDesign(xt, xg, d, sparse = TRUE)
if (gps) {
ld <- ComputeLD(xt, d)
} else {
ld <- matrix(1, k + d, d - 1)
}
H <- NullGD(ld, m, orthonormal = FALSE)
X <- as_matrix(B %*% H)
Xg <- as_matrix(Bg %*% H)
XtX <- base::crossprod(X)
Xty <- base::crossprod(X, y)
U <- chol.default(XtX)
f <- forwardsolve(U, Xty, upper.tri = TRUE, transpose = TRUE)
b <- backsolve(U, f)
yg <- c(Xg %*% b)
ylim <- range(y, yg)
plot.default(x, y, col = 8, ylim = ylim, ann = FALSE, xaxt = "n", yaxt = "n")
abline(v = xd, lty = 2, col = 8)
lines.default(xg, yg, lwd = 2)
title(distr)
}
DemoNull <- function (n, k, gps = FALSE) {
x <- seq.int(0, 1, length.out = n)
y <- rnorm(n, mean = x, sd = 0.2 * sd(x))
op <- par(mfrow = c(2, 3), mar = c(0.25, 0.25, 1.5, 0.25), oma = c(0, 0, 1.5, 0))
on.exit(par(op))
PlotNull(x, y, k, 3, 5, gps)
PlotNull(x, y, k, 5, 5, gps)
PlotNull(x, y, k, 5, 3, gps)
PlotNull(x, y, k, 1, 3, gps)
PlotNull(x, y, k, 0.5, 0.5, gps)
PlotNull(x, y, k, 3, 1, gps)
label <- sprintf("limiting %s P-spline fit at ", c("standard", "general")[gps + 1L])
expr <- substitute(expression(paste(label, lambda == +infinity)), list(label = label))
title(eval(expr), outer = TRUE)
}
DemoPBS <- function (uniform = TRUE) {
if (uniform) xd <- seq.int(0, 0.6, by = 0.1)
else xd <- c(0, 0.75, 1.25, 2, 3.5, 4.5, 5)
x <- MakeGrid(xd, 21)
PBsparse <- pbsDesign(x, xd, 4, 0, sparse = TRUE)
PB <- Zero2NA(PBsparse)
a <- xd[1L]
b <- xd[7L]
period <- b - a
raux <- xd[2:4] + period
xt <- c(xd, raux)
x.wrapped <- x[1:63] + period
B <- PB
B[1:63, 4:6] <- NA_real_
B.wrapped <- PB[1:63, 4:6]
ip <- apply(PB, 2, which.max)
xp.PBSpl <- x[ip]
yp <- PB[cbind(ip, 1:6)]
xp.PBSpl.modified <- xp.PBSpl[4:6]
sub <- xp.PBSpl.modified < xd[4]
xp.PBSpl.modified[sub] <- xp.PBSpl.modified[sub] + period
xp.BSpl <- c(xp.PBSpl[1:3], xp.PBSpl.modified)
xlim <- xt[c(1L, 10L)]
ymax <- max(yp) * 1.15
pch.knots <- rep.int(c(19, 15, 19, 15), c(1, 3, 3, 3))
col.knots <- rep.int(c(1, 8), c(7, 3))
lab1.knots <- vector("expression", 10)
lab1.knots[1] <- expression(a)
for (j in 1:5) {
lab1.knots[j + 1] <- eval(substitute(expression(s[j]), list(j = j)))
}
lab1.knots[7] <- expression(b)
for (j in 1:3) {
lab1.knots[j + 7] <- eval(substitute(expression(s[j]^','), list(j = j)))
}
lab2.knots <- vector("expression", 10)
for (j in 1:10) {
lab2.knots[j] <- eval(substitute(expression(xi[j]), list(j = j)))
}
lab.BSpl <- vector("expression", 6)
for (j in 1:6) {
lab.BSpl[j] <- eval(substitute(expression(B[j]), list(j = j)))
}
op <- par(mfrow = c(2, 1), xpd = TRUE, mar = c(3, 1, 1.25, 1))
on.exit(par(op))
matplot(x, PB, type = "n", ann = FALSE, xlim = xlim, ylim = c(0, ymax),
xaxs = "i", yaxs = "i", xaxt = "n", yaxt = "n", bty = "n")
segments(xt[1], 0, xt[10], 0, col = 8)
polygon(c(a, b, b, a), c(0, 0, ymax, ymax), col = gray(0.3, 0.2), border = NA)
matlines(x, B, lty = rep.int(c(1, 2), c(3, 3)), lwd = 2, col = 1:6)
matlines(x.wrapped, B.wrapped, lty = 2, lwd = 2, col = 4:6)
text.default(xp.BSpl, yp, labels = lab.BSpl, col = 1:6, adj = c(0.2, 0), cex = 1.5)
points.default(xt, numeric(10), pch = pch.knots, col = col.knots)
text.default(xt, numeric(10), labels = lab1.knots, col = col.knots, pos = 1, cex = 1.5)
mtext(lab2.knots, side = 1, line = 2, at = xt, col = col.knots, cex = 1.5)
title("ordinary cubic B-splines")
matplot(x, PB, type = "n", ann = FALSE, xlim = xlim, ylim = c(0, ymax),
xaxs = "i", yaxs = "i", xaxt = "n", yaxt = "n", bty = "n")
segments(a, 0, b, 0, col = 8)
polygon(c(a, b, b, a), c(0, 0, ymax, ymax), col = gray(0.3, 0.2), border = NA)
matlines(x, B, col = 1:6, lty = rep.int(c(1, 2), c(3, 3)), lwd = 2)
matlines(x.wrapped - period, B.wrapped, lty = 2, lwd = 2, col = 4:6)
text.default(xp.PBSpl, yp, labels = lab.BSpl, col = 1:6, adj = c(0.2, 0), cex = 1.5)
points.default(xd, numeric(7), pch = 19)
text.default(xd, numeric(7), labels = lab1.knots[1:7], col = col.knots[1:7],
pos = 1, cex = 1.5)
mtext(lab2.knots[1:7], side = 1, line = 2, at = xd, col = col.knots[1:7], cex = 1.5)
title("periodic cubic B-splines")
}
GetE_periodic <- function (xd, d) {
if (d < 2) stop("d >= 2 required!", call. = FALSE)
degree <- d - 1
k2 <- length(xd)
if (k2 <= d) stop("length(xd) >= d + 1 required!", call. = FALSE)
a <- rep.int(xd[1], degree)
b <- rep.int(xd[k2], degree)
local.knots.a <- c(a, xd[1:(d + 1)])
local.knots.b <- c(xd[(k2 - d):k2], b)
qs <- 0:(degree - 1)
Ca <- splines::splineDesign(local.knots.a, a, d, derivs = qs)
Cb <- splines::splineDesign(local.knots.b, b, d, derivs = qs)
Ca <- Ca[, 1:degree]
Cb <- Cb[, 2:d]
det.a <- prod(Ca[seq.int(1, length.out = degree, by = d)])
det.b <- abs(prod(Cb[seq.int(degree, length.out = degree, by = degree - 1)]))
if (det.a >= det.b) {
E <- forwardsolve(Ca, Cb)
attr(E, "eliminate") <- "head"
} else {
perm <- degree:1
E <- forwardsolve(Cb[, perm], Ca)[perm, ]
attr(E, "eliminate") <- "tail"
}
E
}
AbsorbE_periodic <- function (X, E) {
p <- ncol(X)
degree <- nrow(E)
X.head <- X[, 1:degree]
X.tail <- X[, (p - degree + 1):p]
X.mid <- X[, (degree + 1):(p - degree)]
if (attr(E, "eliminate") == "head") {
PX <- cbind(X.mid, X.tail + X.head %*% E)
} else {
PX <- cbind(X.head + X.tail %*% E, X.mid)
}
PX
}
DemoPBS2 <- function (uniform = TRUE, scaling = TRUE) {
if (uniform) xd <- seq.int(0, 0.6, by = 0.1)
else xd <- c(0, 0.75, 1.25, 2, 3.5, 4.5, 5)
a <- xd[1]
b <- xd[7]
Aknots <- c(rep.int(a, 3), xd, rep.int(b, 3))
x <- gps::MakeGrid(xd, n = 21)
B <- splines::splineDesign(Aknots, x)
E <- GetE_periodic(xd, 4)
PB <- AbsorbE_periodic(B, E)
lab.knots <- vector("expression", 7)
for (j in 1:7) {
lab.knots[j] <- eval(substitute(expression(xi[j]), list(j = j)))
}
lab.BSpl <- vector("expression", 9)
for (j in 1:9) {
lab.BSpl[j] <- eval(substitute(expression(B[j]), list(j = j)))
}
lab.PBSpl <- vector("expression", 9)
for (j in 1:6) {
lab.PBSpl[j] <- eval(substitute(expression(tilde(B)[j]), list(j = j)))
}
ip.BSpl <- apply(B, 2, which.max)
xp.BSpl <- x[ip.BSpl]
yp.BSpl <- B[cbind(ip.BSpl, 1:9)]
ip.PBSpl <- apply(abs(PB), 2, which.max)
xp.PBSpl <- x[ip.PBSpl]
yp.PBSpl <- PB[cbind(ip.PBSpl, 1:6)]
positive <- yp.PBSpl > 0
negative <- yp.PBSpl < 0
if (scaling) {
big <- abs(yp.PBSpl) > 1
big.yp.PBspl <- yp.PBSpl[big]
fctr <- sign(big.yp.PBspl) / max(abs(big.yp.PBspl))
PB[, big] <- PB[, big] * rep(fctr, each = length(x))
yp.PBSpl[big] <- yp.PBSpl[big] * fctr
}
B[B == 0] <- NA
PB[PB == 0] <- NA
op <- par(mfrow = c(2, 1), xpd = TRUE, mar = c(2, 0.75, 1, 0.75))
on.exit(par(op))
matplot(x, B, type = "l", lty = rep(c(2, 1, 2), each = 3), lwd = 2, ann = FALSE,
xaxs = "i", yaxs = "i", xaxt = "n", yaxt = "n", bty = "n", ylim = c(0, 1.1))
segments(a, 0, b, 0, col = 8)
text.default(xp.BSpl, yp.BSpl, labels = lab.BSpl, col = 1:6, adj = c(0.5, 0), cex = 1.5)
points.default(xd, numeric(7), pch = 19)
ystack <- seq.int(0.1, by = 0.1, length.out = 3)
points.default(rep.int(c(a, b), c(3, 3)), c(ystack, ystack), pch = 19, col = 8)
text.default(xd, numeric(7), labels = lab.knots, pos = 1, cex = 1.5)
title("ordinary cubic B-splines")
if (attr(E, "eliminate") == "head") {
lty.PBSpl <- rep(1:2, each = 3)
} else {
lty.PBSpl <- rep(2:1, each = 3)
}
ylim.PBSpl <- range(-0.03, yp.PBSpl) * c(1, 1.1)
matplot(x, PB, type = "l", lty = lty.PBSpl, lwd = 2, ann = FALSE, ylim = ylim.PBSpl,
xaxs = "i", yaxs = "i", xaxt = "n", yaxt = "n", bty = "n")
segments(a, 0, b, 0, col = 8)
text.default(xp.PBSpl[positive], yp.PBSpl[positive], labels = lab.PBSpl[positive],
col = (1:6)[positive], adj = c(0.5, 0), cex = 1.5)
text.default(xp.PBSpl[negative], yp.PBSpl[negative], labels = lab.PBSpl[negative],
col = (1:6)[negative], pos = 1, cex = 1.5)
points.default(xd, numeric(7), pch = 19)
text.default(xd, numeric(7), labels = lab.knots, pos = 1, cex = 1.5)
title("periodic cubic B-splines")
}
DemoSpl <- function (uniform = TRUE) {
if (uniform) {
xd <- 1:6
xt <- seq.int(-2, 9)
xg <- MakeGrid(xt, 21)
b <- c(0.44, 1.11, 1.66, 0.25, 1.6, 1.43, 1.49, 2.52)
} else {
xd <- c(1, 2.06, 2.65, 3.22, 3.91, 6)
xt <- xd[c(1, 1, 1, 1, 2:5, 6, 6, 6, 6)]
xg <- MakeGrid(xd, 21)
b <- c(1.9, 1.84, 1.35, 1.67, 0.48, 1.85, 1.59, 1.35)
}
Bg <- splines::splineDesign(xt, xg, outer.ok = uniform, sparse = TRUE)
Bgdense <- Zero2NA(Bg)
yg <- (Bg %*% b)@x
if (uniform) {
i1 <- 21 * 3 + 1
i2 <- 21 * 8
x <- xg[i1:i2]
y <- yg[i1:i2]
} else {
x <- xg
y <- yg
}
ymax <- max(y)
yd <- y[c(1, 1:5 * 21)]
op <- par(xpd = TRUE, mar = c(2, 2, 1.5, 0.5))
on.exit(par(op))
plot.default(x, y, type = "n", xlim = c(-2, 9), ylim = c(0, ymax),
xaxs = "i", yaxs = "i", ann = FALSE, bty = "n")
polygon(c(1, 6, 6, 1), c(0, 0, ymax, ymax), col = gray(0.3, 0.2), border = NA)
lines.default(x, y, lwd = 2)
points.default(xd, yd, pch = c(17, 19, 19, 19, 19, 17))
matlines(xg, Bgdense, type = "l", lty = 1, lwd = 2, col = c(1:6, 8, 7))
points.default(xd[2:5], numeric(4), pch = 19)
if (uniform) {
points.default(c(-2:0, 7:9), numeric(6), pch = 15)
} else {
segments(-2, 0, 1, 0, lwd = 2, col = 1)
segments(6, 0, 9, 0, lwd = 2, col = 7)
segments(1, 1, 1, 0, lwd = 2, col = 1, lty = 2)
segments(6, 1, 6, 0, lwd = 2, col = 7, lty = 2)
ystack <- seq.int(0.1, by = 0.1, length.out = 3)
points.default(c(1, 1, 1, 6, 6, 6), c(ystack, ystack), pch = 15)
}
points.default(c(1, 6), numeric(2), pch = 17)
type <- c("non-uniform", "uniform")[uniform + 1L]
title(sprintf("cubic spline with %s knots", type))
dB <- splines::splineDesign(xt, rep.int(xd[1:5], 3), derivs = rep(1:3, each = 5))
pc <- c(dB %*% b) * rep(1 / factorial(1:3), each = 5)
lab1 <- paste0("f", 1:5)
lab2 <- c("constant", "linear", "quadratic", "cubic")
pc <- matrix(c(yd[1:5], pc), ncol = 4, dimnames = list(lab1, lab2))
list(xd = xd, bspl.coef = b, piecepoly.coef = pc)
}
CheckSupp <- function (x, xt, d) {
degree <- d - 1
K <- length(xt)
xd <- xt[seq.int(d, K - degree)]
nI <- length(xd) - 1L
id <- findInterval(x, xd, rightmost.closed = TRUE)
if (min(id) < 1L || max(id) > nI) {
stop("Domain does not contain all x-values!", call. = FALSE)
}
ni <- tabulate(id, nbins = nI)
.Call("C_CheckSupp", ni, d, PACKAGE = "gps") > 0L
}
gps2bspl <- function (x, xt, d, m, gps = TRUE, periodic = FALSE,
unique.x.provided = FALSE) {
if (d < 2) stop("d >= 2 required!", call. = FALSE)
if (m < 1 || m > d - 1) stop("1 <= m <= d - 1 required!", call. = FALSE)
xt <- MonoKnots(xt, d)
p <- length(xt) - d
xu <- if (unique.x.provided) x else unique.default(x)
if (p > length(xu)) {
stop("There are more B-splines than unique x-values!", call. = FALSE)
}
if (CheckSupp(xu, xt, d)) {
stop("There are no x-values on some B-spline's support!", call. = FALSE)
}
B <- splines::splineDesign(xt, x, d, sparse = TRUE)
D <- SparseD(xt, d, m, gps)[[1L]]
Dt <- Csparse2LTB(D)
S <- DtD(Dt)
if (ncol(S) > p) S <- S[, 1:p]
if (periodic) {
C <- GetPBC(xt, d, compact = FALSE, sparse = FALSE)
} else {
C <- matrix(0, nrow = 0, ncol = p)
}
list(B = B, Dt = Dt, S = S, xt = xt, d = d, m = m, C = C)
}
gps2red <- function (bspl, y, w = NULL, scalePen = TRUE) {
y <- as.double(y)
X <- bspl$B
n <- X@Dim[1L]
if (length(y) != n) stop("Incorrect number of y-values!")
if (!is.null(w)) {
if (length(w) != n) stop("Incorrect number of weights!")
w <- (1 / sum(w)) * w
w <- sqrt(w)
X@x <- w[X@i + 1L] * X@x
y <- w * y
}
if (scalePen) {
scalePen.fctr <- VecDot(X@x) / VecDot(bspl$Dt)
Dt <- VecScal(sqrt(scalePen.fctr), bspl$Dt, overwrite = FALSE)
S <- VecScal(scalePen.fctr, bspl$S, overwrite = FALSE)
} else {
scalePen.fctr <- 1
Dt <- bspl$Dt
S <- bspl$S
}
Z <- Matrix::crossprod(X)
z <- Matrix::crossprod(X, y)@x
Z <- Csparse2LTB(Z)
L <- try(LPBTRF(Z, overwrite = FALSE), silent = TRUE)
if (inherits(L, "try-error")) {
stop("B-spline design matrix does not have full column rank!")
}
f <- SolveLTB(transA = FALSE, L, z, overwrite = FALSE)
minRSS <- VecDot(y) - VecDot(f)
E <- .Call("C_FormE", L, Dt, PACKAGE = "gps")
list(n = n, Z = Z, z = z, L = L, f = f, minRSS = minRSS,
scalePen.fctr = scalePen.fctr, Dt = Dt, S = S, E = E)
}
gps2DR <- function (red, tol = 1e-6) {
t0 <- Sys.time()
q <- ncol(red$E)
dbar <- .Call("C_MeanDR", red$E, PACKAGE = "gps")
dq <- .Call("C_MinDR", red$E, tol, PACKAGE = "gps")
d1 <- .Call("C_MaxDR", red$L, red$Dt, tol, PACKAGE = "gps")
dq.cutoff <- d1 * (.Machine$double.eps / 2)
if (dq < dq.cutoff) {
dq <- dq.cutoff
warning("The smallest eigenvalue is erratic as E'E is numerically singular!")
}
approx.DR <- .Call("C_ApproxDR", q, dq, d1, dbar, tol, PACKAGE = "gps")
t1 <- Sys.time()
secs <- difftime(t1, t0, units = "secs")[[1L]]
list(mean = dbar, min = dq, max = d1, approx.values = approx.DR, secs = secs)
}
ExactDR <- function (E) {
t0 <- Sys.time()
A <- .Call("C_LAUUM", E, PACKAGE = "gps")
exact.DR <- base::eigen(A, symmetric = TRUE, only.values = TRUE)$values
q <- length(exact.DR)
d1 <- exact.DR[1L]
dq.cutoff <- d1 * (.Machine$double.eps / 2)
ind <- which(exact.DR < dq.cutoff)
if (length(ind)) {
exact.DR[ind] <- dq.cutoff
warning("Small eigenvalues are erratic as E'E is numerically singular!")
}
t1 <- Sys.time()
secs <- difftime(t1, t0, units = "secs")[[1L]]
list(values = exact.DR, secs = secs)
}
Rho2REDF <- function (DR, rho) {
.Call("C_Rho2REDF", DR, rho, PACKAGE = "gps")
}
REDF2Rho <- function (DR, redf, rho0, MaxNewton, tol = 1e-6) {
.Call("C_REDF2Rho", DR, redf, rho0, MaxNewton, tol, PACKAGE = "gps")
}
gps2RhoLim <- function (DR, kappa = 0.01) {
lam.min <- kappa / (1 - kappa) / DR$mean
rho.min <- log(lam.min)
lam.max <- (1 - kappa) / kappa / DR$min
rho.max <- log(lam.max)
approx.DR <- DR$approx.values
if (is.na(approx.DR[1L])) {
rho.max.heuristic <- NA_real_
} else {
redf <- length(approx.DR) * kappa
rho.mid <- 0.5 * (rho.min + rho.max)
MaxNewton <- 0.25 * (rho.max - rho.min)
rho.max.heuristic <- REDF2Rho(approx.DR, redf, rho.mid, MaxNewton)
}
list(min = rho.min, max = rho.max, max.heuristic = rho.max.heuristic)
}
Q1fun <- function (z, a = 0, b = 1) {
if (a >= b) stop("a < b required!")
h <- z * z - z
theta <- exp(a + (b - a) * z)
bounds <- c(0, b - a)
list(h = h, theta = theta, bounds = bounds)
}
Q2fun <- function (z, a = 0, b = 1) {
if (a >= b) stop("a < b required!")
y <- 1 - z
B0 <- y * y * y
B1 <- 3 * z * y * y
B2 <- 3 * z * z * y
B3 <- z * z * z
h <- B1 - B2
theta <- exp((B0 + B2) * a + (B2 + B3) * b)
bounds <- c(a, (2 * a + b) / 3)
list(h = h, theta = theta, bounds = bounds)
}
map01 <- function (x) (1 / (x[1L] - x[length(x)])) * (x - x[length(x)])
ApproxDR <- function (DR, Qfun, gamma.grid = seq.int(0, 1, 0.05), verbose = TRUE) {
if (!is.function(Qfun)) stop("'Qfun' is not a function!")
Qz <- Qfun(z = 0.5, a = 0, b = 1)
if (anyNA(match(c("h", "theta", "bounds"), names(Qz)))) {
stop("'Qfun' needs to return elements 'h', 'theta' and 'bounds'!")
}
q <- length(DR)
a <- log(DR[q])
b <- log(DR[1L])
dbar <- sum(DR) / q
logDR <- map01(log(DR))
step <- 1 / (q + 1)
p <- seq.int(from = step, by = step, length.out = q)
aggregated.DR <- numeric(q)
success <- 0
op <- par(mfrow = c(1, 2), mar = c(2, 2, 1.5, 0.5))
on.exit(par(op))
for (gamma in gamma.grid) {
z <- map01(log(1 - p) - gamma * log(p))
zlab <- sprintf("z(%.2f)", gamma)
Qz <- Qfun(z, a, b)
plot.default(p, logDR, type = "l", ann = FALSE)
lines.default(p, z, lty = 2)
title(sprintf("log(d) and %s", zlab))
plot.default(z, logDR, type = "l", ann = FALSE)
title(sprintf("log(d) ~ %s", zlab))
theta <- Qz$theta
h <- Qz$h
eta <- theta * h
bounds <- Qz$bounds
tmpDR <- .Call("C_RootApproxDR", theta, h, eta, dbar, bounds, PACKAGE = "gps")
if (tmpDR[1L]) {
success <- success + 1
aggregated.DR <- aggregated.DR + tmpDR
lines.default(z, map01(log(tmpDR)), lty = 2)
}
if (verbose) readline("Hit <Enter> to continue: ")
if (success && tmpDR[1L] == 0) break
}
list(success = success, aggregated.DR = aggregated.DR, gamma.exit = gamma)
}
DemoApproxDR <- function (DR, Qfun = NULL, verbose = TRUE) {
gamma.grid <- seq.int(0, 1, 0.05)
if (is.null(Qfun)) {
out1 <- ApproxDR(DR, Q1fun, gamma.grid, verbose)
gamma.exit <- out1$gamma.exit
if (gamma.exit < 1) {
gamma.grid <- seq.int(gamma.exit, 1, 0.05)
} else {
gamma.grid <- numeric(0)
}
out2 <- ApproxDR(DR, Q2fun, gamma.grid, verbose)
success <- out1$success + out2$success
approx.DR <- (out1$aggregated.DR + out2$aggregated.DR) / success
} else {
out <- ApproxDR(DR, Qfun, gamma.grid, verbose)
success <- out$success
approx.DR <- out$aggregated.DR / success
}
if (success) {
op <- par(mfrow = c(1, 2), mar = c(2, 2, 1.5, 0.5))
on.exit(par(op))
q <- length(DR)
step <- 1 / (q + 1)
p <- seq.int(from = step, by = step, length.out = q)
logDR <- map01(log(DR))
result <- list(min = DR[q], mean = sum(DR) / q, approx.values = approx.DR)
rho.lim <- gps2RhoLim(result)
rho.grid <- seq.int(rho.lim$min, rho.lim$max, length.out = 100)
redf <- Rho2REDF(DR, rho.grid)
approx.redf <- Rho2REDF(approx.DR, rho.grid)
redf.heuristic <- Rho2REDF(DR, rho.lim$max.heuristic)
mapping <- sprintf("loose: %.2f -> %.2f\nheuristic: %.2f -> %.2f",
rho.lim$max, redf[100L], rho.lim$max.heuristic, redf.heuristic)
plot.default(p, logDR, type = "l", ann = FALSE)
lines.default(p, map01(log(approx.DR)), lty = 2)
title("log(d)")
plot.default(rho.grid, redf, type = "l", ann = FALSE)
lines.default(rho.grid, approx.redf, lty = 2)
points.default(rho.lim$max.heuristic, redf.heuristic, pch = 19)
text.default(rho.lim$max, redf[1L], labels = mapping, adj = c(1, 1))
title("redf ~ rho")
} else {
warning("Unable to approximate Demmler-Reinsch eigenvalues!")
}
}
GridPWLS <- function (bspl, red, rho) {
.Call("C_GridPWLS", red$Z, red$L, red$S, red$z, t.default(bspl$C), rho,
PACKAGE = "gps")
}
GridGCV <- function (red, pwls) {
.Call("C_GridGCV", red$n, red$L, red$f, red$minRSS, pwls$beta, pwls$edf,
PACKAGE = "gps")
}
gps2GS <- function (x, y, w = NULL, xt, d = 4, m = 2, gps = TRUE,
periodic = FALSE, ng = 20, scalePen = TRUE) {
bspl <- gps2bspl(x, xt, d, m, gps, periodic)
red <- gps2red(bspl, y, w, scalePen)
DR <- gps2DR(red)
rho.lim <- gps2RhoLim(DR)
if (ng > 0) {
if (is.na(rho.lim$max.heuristic)) {
rho <- seq.int(rho.lim$min, rho.lim$max, length.out = ng)
} else {
rho <- seq.int(rho.lim$min, rho.lim$max.heuristic, length.out = ng)
}
t0 <- Sys.time()
pwls <- GridPWLS(bspl, red, rho)
pwls$gcv <- GridGCV(red, pwls)
t1 <- Sys.time()
pwls$secs <- difftime(t1, t0, units = "secs")[[1L]]
} else {
pwls <- NULL
}
S <- LTB2Csparse(bspl$S, symmetric = TRUE)
BWB <- LTB2Csparse(red$Z, symmetric = TRUE)
eqn <- list(B = bspl$B, BWB = BWB, BWy = red$z,
omega = red$scalePen.fctr, S = S, C = bspl$C)
list(eqn = eqn, eigen = DR, rho.lim = rho.lim, E = red$E, pwls = pwls)
}
DemoRhoLim <- function (fit, plot = TRUE) {
DR <- ExactDR(fit$E)
approx.DR <- fit$eigen[c("approx.values", "secs")]
names(approx.DR) <- c("values", "secs")
q <- length(DR$values)
rho.min.star <- fit$rho.lim$min
rho.max.star <- fit$rho.lim$max
redf.min <- 0.01 * q
redf.max <- 0.99 * q
rho.mid <- (rho.min.star + rho.max.star) / 2
MaxNewton <- 0.25 * (rho.max.star - rho.min.star)
rho.min <- REDF2Rho(DR$values, redf.max, rho.mid, MaxNewton)
rho.max <- REDF2Rho(DR$values, redf.min, rho.mid, MaxNewton)
ng <- 100L
rho.grid <- seq.int(rho.min.star, rho.max.star, length.out = ng)
redf.exact <- Rho2REDF(DR$values, rho.grid)
redf.min.star <- redf.exact[ng]
redf.max.star <- redf.exact[1L]
rho.max.heuristic <- fit$rho.lim$max.heuristic
if (is.na(rho.max.heuristic)) {
redf.approx <- rep.int(NA_real_, ng)
redf.min.heuristic <- NA_real_
warning("No approximated eigenvalues and heuristic upper bound for 'rho'!")
} else {
redf.approx <- Rho2REDF(approx.DR$values, rho.grid)
redf.min.heuristic <- Rho2REDF(DR$values, rho.max.heuristic)
}
map.exact <- sprintf("[%.1f, %.1f] -> [%.1f, %.1f]",
rho.min, rho.max, redf.min, redf.max)
map.wider <- sprintf("[%.1f, %.1f] -> [%.1f, %.1f]",
rho.min.star, rho.max.star, redf.min.star, redf.max.star)
map.heuristic <- sprintf("%.1f -> %.1f", rho.max.heuristic, redf.min.heuristic)
if (plot) {
op <- par(mfrow = c(1, 2), mar = c(2, 2, 1.5, 0.5))
on.exit(par(op))
p <- (1:q) / (q + 1)
plot.default(p, log(DR$values), type = "l", ann = FALSE)
lines.default(p, log(approx.DR$values), lty = 2)
text.default(p[q], log(DR$values[1L]), adj = c(1, 1),
labels = sprintf("exact: solid\napprox: dashed"))
title("log eigenvalues")
plot.default(rho.grid, redf.exact, type = "l", ann = FALSE)
lines.default(rho.grid, redf.approx, lty = 2)
segments(rho.min, q, rho.min, 0.75 * q, col = 8, lty = 2)
segments(rho.max, 0, rho.max, 0.25 * q, col = 8, lty = 2)
points.default(rho.max.heuristic, redf.min.heuristic, pch = 19)
text.default(rho.max.star, redf.max.star, adj = c(1, 1),
labels = sprintf("exact: %s\nwider: %s\nheuristic: %s",
map.exact, map.wider, map.heuristic))
title("redf ~ rho")
}
eigen <- list(approx = approx.DR, exact = DR)
limit <- list(rho = list(exact = c(rho.min, rho.max),
wider = c(rho.min.star, rho.max.star),
heuristic = rho.max.heuristic),
redf = list(exact = c(redf.min, redf.max),
wider = c(redf.min.star, redf.max.star),
heuristic = redf.min.heuristic))
mapping <- list(exact = map.exact, wider = map.wider, heuristic = map.heuristic)
redf <- list(rho = rho.grid, exact = redf.exact, approx = redf.approx)
list(eigen = eigen, limit = limit, mapping = mapping, redf = redf)
}
Diff <- function (x, k = 1L, n = length(x), xi = 1L) {
.Call("C_Diff", x, k, n, xi, PACKAGE = "gps")
}
SparseDelta <- function (r) {
r <- as.integer(r)
x <- rep.int(c(-1, 1), r)
i <- rep(seq.int(0L, r - 1L), each = 2L)
p <- c(0L, seq.int(1L, length.out = r, by = 2L), 2L * r)
methods::new("dgCMatrix", i = i, p = p, Dim = c(r, r + 1L), x = x)
}
SparseWtDelta <- function (h) {
r <- length(h)
x <- rep.int(c(-1, 1), r) * rep(1 / h, each = 2)
i <- rep(seq.int(0L, r - 1L), each = 2L)
p <- c(0L, seq.int(1L, length.out = r, by = 2L), 2L * r)
methods::new("dgCMatrix", i = i, p = p, Dim = c(r, r + 1L), x = x)
}
SparseGD <- function (xt, d, m) {
K <- length(xt)
D <- vector("list", m)
h <- Diff(x = xt, k = d - 1L, n = K - 2L, xi = 2L)
D[[1L]] <- SparseWtDelta(h)
i <- 2L
while (i <= m) {
h <- Diff(x = xt, k = d - i, n = K - 2L * i, xi = i + 1L)
D[[i]] <- SparseWtDelta(h) %*% D[[i - 1L]]
i <- i + 1L
}
D
}
DiffCoef <- function (b, xt, d, m) {
if (d < 2) stop("d >= 2 required!", call. = FALSE)
if (m < 1 || m >= d) stop("1 <= m <= d - 1 required!", call. = FALSE)
xt <- MonoKnots(xt, d)
K <- length(xt)
if (length(b) != K - d) {
stop("length(b) == length(xt) - d required!", call. = FALSE)
}
b <- as.double(b)
for (i in 1:m) {
h <- Diff(x = xt, k = d - i, n = K - 2L * i, xi = i + 1L)
b <- Diff(b) / h
}
b
}
ComputeLD <- function (xt, d) {
.Call("C_ComputeLD", xt, d, PACKAGE = "gps")
}
NullGD <- function (ld, m = 1, orthonormal = TRUE) {
basis <- .Call("C_NullGD", ld, m, PACKAGE = "gps")
if (orthonormal && (m > 1)) {
Q <- qr.Q(qr.default(basis[, m:1]))[, m:1]
i <- sequence.default(1:(m - 1))
j <- rep.int(2:m, 1:(m - 1))
Q[cbind(i, j)] <- 0
basis <- Q
}
basis
}
QuadWts <- function (ord) {
if (ord == 1) {
return(2)
}
p <- seq.int(0, ord - 1)
P <- outer(seq.int(-1, 1, length.out = ord), p, "^")
Pinv <- solve.default(P)
pow <- outer(1:ord, p, "+")
H <- (1 - (-1) ^ pow) / pow
base::crossprod(Pinv, H %*% Pinv)
}
SbarBlocks <- function (xt, d, m) {
if (d < 2) stop("d >= 2 required!", call. = FALSE)
if (m < 0 || m >= d) stop("0 <= m <= d - 1 required!", call. = FALSE)
xt <- MonoKnots(xt, d)
K <- length(xt)
ord <- d - m
if (ord == 1) {
h <- Diff(x = xt, n = K - 2 * (d - 1), xi = d)
return(h)
}
W <- QuadWts(ord)
xd <- xt[seq.int(d, K - d + 1)]
xg <- MakeGrid(xd, ord)
xt.local <- xt[seq.int(1 + m, K - m)]
B <- splines::splineDesign(xt.local, xg, ord, sparse = TRUE)
.Call("C_SbarBlocks", xd, W, B@x, PACKAGE = "gps")
}
GramBS <- function (xt, d) {
blocks <- SbarBlocks(xt, d, m = 0)
LTB <- .Call("C_SbarLTB", blocks, PACKAGE = "gps")
LTB2Csparse(LTB, symmetric = TRUE)
}
btSb <- function (b, xt, d, m) {
db <- DiffCoef(b, xt, d, m)
blocks <- SbarBlocks(xt, d, m)
.Call("C_btSb", blocks, db, PACKAGE = "gps")
}
SparseD <- function (xt, d, m = NULL, gps = TRUE) {
d <- as.integer(d)
if (d < 2L) stop("d >= 2 required!", call. = FALSE)
if (is.null(m)) {
m <- seq_len(d - 1L)
} else {
m <- sort.int(as.integer(m), method = "radix")
}
m.min <- m[1L]
m.max <- m[length(m)]
if (m.min < 1L || m.max >= d) {
stop("1 <= m <= d - 1 required!", call. = FALSE)
}
NAMES <- sprintf("order.%d", m)
xt <- MonoKnots(xt, d)
D <- SparseGD(xt, d, m.max)
D <- D[m]
names(D) <- NAMES
if (gps) return(D)
nm <- length(m)
Sbar <- vector("list", nm); names(Sbar) <- NAMES
K <- vector("list", nm); names(K) <- NAMES
for (i in 1:nm) {
m_i <- m[i]
D_i <- D[[i]]
blocks <- SbarBlocks(xt, d, m_i)
if (is.array(blocks)) {
Sbar.LTB <- .Call("C_SbarLTB", blocks, PACKAGE = "gps")
Sbar.Matrix <- LTB2Csparse(Sbar.LTB, symmetric = TRUE)
L.LTB <- LPBTRF(Sbar.LTB, overwrite = TRUE)
L.Matrix <- LTB2Csparse(L.LTB)
K_i <- Matrix::crossprod(L.Matrix, D_i)
} else {
k1 <- length(blocks)
Sbar.Matrix <- methods::new("ddiMatrix", diag = "N", x = blocks, Dim = c(k1, k1))
L.diag <- sqrt(blocks)
newx <- L.diag[D_i@i + 1L] * D_i@x
K_i <- methods::new("dgCMatrix", i = D_i@i, p = D_i@p, Dim = D_i@Dim, x = newx)
}
Sbar[[i]] <- Sbar.Matrix
K[[i]] <- K_i
}
sandwich <- list(D = D, Sbar = Sbar)
structure(K, sandwich = sandwich)
}
PriorCoef1 <- function (n, D) {
q <- D@Dim[1L]
S <- as_matrix(Matrix::crossprod(D))
ei <- base::eigen(S, symmetric = TRUE)
qs <- seq_len(q)
d <- ei$values[qs]
U <- ei$vectors[, qs]
di <- 1 / d
e <- rnorm(n * q)
if (n > 1) e <- ChangeDim(e, c(q, n))
b <- U %*% (sqrt(di) * e)
if (n == 1) b <- c(b)
b
}
PriorCoef2 <- function (n, D) {
Dim <- D@Dim
q <- Dim[1L]
p <- Dim[2L]
m <- p - q
e <- rnorm(n * q)
if (n > 1) e <- ChangeDim(e, c(q, n))
Dt <- Csparse2LTB(D)
y <- SolveLTB(transA = TRUE, Dt, e, overwrite = TRUE)
R <- as_matrix(D[, seq.int(q + 1, p), drop = FALSE])
X <- SolveLTB(transA = TRUE, Dt, R, overwrite = TRUE)
XtX <- base::crossprod(X)
Xty <- base::crossprod(X, y)
U <- chol.default(XtX + diag(m))
f <- forwardsolve(U, Xty, upper.tri = TRUE, transpose = TRUE)
b <- backsolve(U, f)
r <- y - X %*% b
if (n > 1) rbind(r, b) else c(r, b)
}
PriorCoef3 <- function (n, D) {
Dim <- D@Dim
q <- Dim[1L]
p <- Dim[2L]
QR <- Dt2ThinQR(D, form.Q = FALSE)
Rt <- QR$Rt
e <- rnorm(n * q)
if (n > 1) e <- ChangeDim(e, c(q, n))
x <- SolveLTB(transA = FALSE, Rt, e, overwrite = TRUE)
x <- SolveLTB(transA = TRUE, Rt, x, overwrite = TRUE)
b <- Matrix::crossprod(D, x)
if (n == 1) b@x else as_matrix(b)
}
PriorCoef <- function (n, D) {
b <- PriorCoef3(n, D)
sig <- sqrt(MPinv0(D, only.diag = TRUE))
scale.fctr <- 1 / sig[1]
b <- VecScal(scale.fctr, b, overwrite = TRUE)
sig <- VecScal(scale.fctr, sig, overwrite = TRUE)
list(coef = b, sigma = sig)
}
MPinvUUt <- function (D, method = "qr") {
if (method == "qr") {
QR <- Dt2ThinQR(D, form.Q = TRUE)
Qt <- QR$Qt
Rt <- QR$Rt
SolveLTB(transA = TRUE, Rt, Qt, overwrite = TRUE)
} else if (method == "eigen") {
q <- D@Dim[1L]
S <- as_matrix(Matrix::crossprod(D))
ei <- base::eigen(S, symmetric = TRUE)
qs <- seq_len(q)
d <- ei$values[qs]
Q <- ei$vectors[, qs]
di <- 1 / d
sqrt(di) * t.default(Q)
} else {
stop("'method' must be either \"qr\" or \"eigen\"!", call. = FALSE)
}
}
MPinv0 <- function (D, only.diag = FALSE) {
Ut <- MPinvUUt(D, method = "qr")
if (only.diag) {
base::colSums(Ut * Ut)
} else {
base::crossprod(Ut)
}
}
MPinv <- function (D, only.diag = FALSE) {
V <- MPinv0(D, only.diag)
VecScal(1 / V[1L], V, overwrite = TRUE)
}
Prior2Cond <- function (k) {
rho <- seq.int(0.05, 1, 0.05)
N <- length(rho)
cond.X <- numeric(N)
for (i in 1:N) {
Dt <- matrix(c(rho[i], -(1 + rho[i]), 1), nrow = 3, ncol = k + 2)
R <- matrix(0, nrow = k + 2, ncol = 2L)
R[c(k + 1, k + 2, 2 * k + 4)] <- c(1, -(1 + rho[i]), 1)
X <- SolveLTB(transA = TRUE, Dt, R, overwrite = TRUE)
XtX <- base::crossprod(X)
ka <- try(kappa(X, exact = TRUE), silent = TRUE)
cond.X[i] <- if (inherits(ka, "try-error")) Inf else ka
}
cond.XtX <- cond.X ^ 2
op <- par(mfrow = c(1, 2), mgp = c(1.5, 0.5, 0), mar = c(2.5, 2.5, 1.5, 0.5))
on.exit(op)
plot.default(rho, log10(cond.X), type = "b", xlab = expression(rho),
ylab = "log10(condition number)", main = "X")
plot.default(rho, log10(cond.XtX), type = "b", xlab = expression(rho),
ylab = "log10(condition number)", main = "X'X")
list(rho = rho, X = cond.X, XtX = cond.XtX)
}
GetPBC <- function (xt, d, compact = TRUE, transpose = FALSE, sparse = TRUE) {
d <- as.integer(d)
if (d < 2L) stop("d >= 2 required!")
degree <- d - 1L
K <- length(xt)
if (length(xt) < degree + 2L * d) {
stop("length(xt) >= 3 * d - 1 required!")
}
p <- K - d
a <- xt[d]
b <- xt[K - degree]
nDeriv <- seq.int(0L, length.out = degree)
knots.a <- xt[seq.int(from = 1L, length.out = 2L * d)]
knots.b <- xt[seq.int(to = K, length.out = 2L * d)]
Ca <- splines::splineDesign(knots.a, rep.int(a, degree), d, nDeriv)
Cb <- splines::splineDesign(knots.b, rep.int(b, degree), d, nDeriv)
Ca <- Ca[, seq.int(1L, degree)]
Cb <- Cb[, seq.int(2L, d)]
C.compact <- cbind(Ca, -Cb, deparse.level = 0L)
Ct.compact <- t.default(C.compact)
if (compact) {
if (transpose) Ct.compact else C.compact
} else {
if (sparse) {
if (transpose) {
Ct.i <- rep.int(c(nDeriv, seq.int(to = p - 1L, length.out = degree)), degree)
Ct.p <- seq.int(0L, by = 2L * degree, length.out = d)
ChangeDim(Ct.compact, NULL)
Ct.x <- Ct.compact
methods::new("dgCMatrix", i = Ct.i, p = Ct.p, Dim = c(p, degree), x = Ct.x)
} else {
blocksize <- degree * degree
C.i <- rep.int(nDeriv, 2L * degree)
C.p <- c(seq.int(0, by = degree, length.out = d),
rep.int(blocksize, p - 2L * degree),
seq.int(blocksize + degree, by = degree, length.out = degree))
ChangeDim(C.compact, NULL)
C.x <- C.compact
methods::new("dgCMatrix", i = C.i, p = C.p, Dim = c(degree, p), x = C.x)
}
} else {
ind <- c(seq.int(1L, degree), seq.int(to = p, length.out = degree))
if (transpose) {
Ct <- matrix(0, p, degree)
Ct[ind, ] <- Ct.compact
Ct
} else {
C <- matrix(0, degree, p)
C[, ind] <- C.compact
C
}
}
}
}
NullPBC <- function (xt, d, compact = TRUE) {
d <- as.integer(d)
degree <- d - 1L
Ct <- GetPBC(xt, d, compact = TRUE, transpose = TRUE)
QR <- qr.default(Ct)
Q <- qr.Q(QR, complete = TRUE)
Q.compact <- Q[, seq.int(d, 2L * degree), drop = FALSE]
if (compact) {
Q.compact
} else {
p <- length(xt) - d
pp <- p - degree
N <- pp - degree
Q.compact.i <- c(seq.int(0L, length.out = degree),
seq.int(to = p - 1L, length.out = degree))
Q.i <- c(seq.int(degree, length.out = N), rep.int(Q.compact.i, degree))
Q.p <- c(seq.int(0, length.out = N),
seq.int(N, by = 2L * degree, length.out = d))
Q.x <- c(rep.int(1, N), Q.compact)
methods::new("dgCMatrix", i = Q.i, p = Q.p, Dim = c(p, pp), x = Q.x)
}
}
pbsDesign <- function (x, xd, d, nDeriv = 0, sparse = FALSE, wrap = TRUE) {
d <- as.integer(d)
if (d < 2L) stop("d >= 2 required!")
degree <- d - 1L
if (nDeriv < 0 || nDeriv > degree) stop("0 <= nDeriv <= d - 1 required!")
xd <- as.double(xd)
if (!IsMonoInc(xd)) stop("'xd' is not strictly ascending!")
k2 <- length(xd)
if (k2 <= d) stop("length(xd) >= d + 1 required!")
a <- xd[1L]
b <- xd[k2]
if (min(x) < a || max(x) > b) stop("domain does not contain all x-values!")
if (wrap) {
pbsDesign1(x, xd, d, nDeriv, sparse)
} else {
pbsDesign2(x, xd, d, nDeriv, sparse)
}
}
pbsDesign1 <- function (x, xd, d, nDeriv = 0, sparse = FALSE) {
degree <- d - 1L
k2 <- length(xd)
p <- k2 - 1L
a <- xd[1L]
b <- xd[k2]
period <- b - a
raux <- xd[2:d] + period
xt <- c(xd, raux)
b.aux <- raux[degree]
ind.x.wrapped <- which(x < xd[d])
if (length(ind.x.wrapped) == 0L) {
PB <- splines::splineDesign(xt, x, d, nDeriv, sparse = sparse)
} else {
nx <- length(x)
x.wrapped <- x[ind.x.wrapped] + period
x.wrapped[x.wrapped < b] <- b
x.wrapped[x.wrapped > b.aux] <- b.aux
x.padded <- c(x, x.wrapped)
xt0 <- c(rep.int(a, degree), xt, rep.int(b.aux, degree))
B <- splines::splineDesign(xt0, x.padded, d, nDeriv, sparse = TRUE)
B.p <- B@p[seq.int(d, length(B@p) - degree)]
ind <- seq.int(B.p[1L] + 1L, B.p[p - degree + 1L])
unchanged.i <- B@i[ind]
unchanged.x <- B@x[ind]
B.p.vital <- B.p[seq.int(to = p + 1L, length.out = d)]
ind <- seq.int(B.p.vital[1L] + 1L, B.p.vital[d])
unordered.i <- B@i[ind]
unordered.x <- B@x[ind]
ind <- (unordered.i >= nx)
unordered.i[ind] <- ind.x.wrapped[unordered.i[ind] + 1L - nx] - 1L
unordered.j <- rep.int(1:degree, diff.default(B.p.vital))
ind <- order(unordered.j, unordered.i, method = "radix")
reordered.i <- unordered.i[ind]
reordered.x <- unordered.x[ind]
PB.p <- B.p - B.p[1L]
PB.i <- c(unchanged.i, reordered.i)
PB.x <- c(unchanged.x, reordered.x)
if (sparse) {
PB <- methods::new("dgCMatrix", i = PB.i, p = PB.p, Dim = c(nx, p), x = PB.x)
} else {
PB.j <- rep.int(1:p, diff.default(PB.p))
PB <- matrix(0, nx, p)
PB[cbind(PB.i + 1L, PB.j)] <- PB.x
}
}
PB
}
pbsDesign2 <- function (x, xd, d, nDeriv = 0, sparse = FALSE) {
degree <- d - 1L
k2 <- length(xd)
xt <- c(rep.int(xd[1L], degree), xd, rep.int(xd[k2], degree))
B <- splines::splineDesign(xt, x, d, nDeriv, sparse = TRUE)
p <- ncol(B)
ind <- seq.int(B@p[d] + 1L, B@p[p + 1L - degree])
PB.i <- B@i[ind]
PB.x <- B@x[ind]
PB.p <- B@p[seq.int(d, p + 1L)]
PB.p <- PB.p - PB.p[1L]
IND <- seq.int(to = k2, length.out = degree)
PB.p[IND] <- PB.p[k2 - degree]
bnd_cols <- c(seq.int(1L, degree), seq.int(k2, p))
nnz_per_col <- B@p[bnd_cols + 1L] - B@p[bnd_cols]
active_col <- nnz_per_col > 0L
bnd_cols <- bnd_cols[active_col]
num_active_cols <- length(bnd_cols)
if (num_active_cols > 0L) {
Q <- NullPBC(xt, d)
Q <- Q[active_col, , drop = FALSE]
i_lst <- x_lst <- vector("list", num_active_cols)
for (i in 1:num_active_cols) {
j <- bnd_cols[i]
ind <- seq.int(B@p[j] + 1L, B@p[j + 1L])
i_lst[[i]] <- B@i[ind]
x_lst[[i]] <- B@x[ind]
}
active_row <- sort.int(unique.default(unlist(i_lst)), method = "radix")
num_active_rows <- length(active_row)
block.i <- unlist(lapply(i_lst, match, active_row))
block.j <- rep.int(1:num_active_cols, nnz_per_col[active_col])
block.x <- unlist(x_lst)
B.boundary <- matrix(0, num_active_rows, num_active_cols)
B.boundary[cbind(block.i, block.j)] <- block.x
B.transformed <- B.boundary %*% Q
PB.i <- c(PB.i, rep.int(active_row, degree))
PB.x <- c(PB.x, B.transformed)
PB.p[IND] <- PB.p[IND] + num_active_rows * (1:degree)
}
nr <- length(x)
nc <- k2 - 1L
if (sparse) {
PB <- methods::new("dgCMatrix", i = PB.i, p = PB.p, Dim = c(nr, nc), x = PB.x)
} else {
PB.j <- rep.int(1:nc, diff.default(PB.p))
PB <- matrix(0, nr, nc)
PB[cbind(PB.i + 1L, PB.j)] <- PB.x
}
PB
}
SparsePD <- function (xd, d, wrap = TRUE) {
d <- as.integer(d)
if (d < 2L) stop("d >= 2 required!")
degree <- d - 1L
xd <- as.double(xd)
if (!IsMonoInc(xd)) stop("'xd' is not strictly ascending!")
k2 <- length(xd)
if (k2 <= d) stop("length(xd) >= d + 1 required!")
if (wrap) {
SparsePD1(xd, d)
} else {
SparsePD2(xd, d)
}
}
SparsePD1 <- function (xd, d) {
degree <- d - 1L
k2 <- length(xd)
p <- k2 - 1L
a <- xd[1L]
b <- xd[k2]
period <- b - a
laux <- xd[(k2 - degree):p] - period
raux <- xd[2:d] + period
xt <- c(laux, xd, raux)
D <- SparseD(xt, d)
PD <- vector("list", degree)
for (m in 1:degree) {
PD[[m]] <- D[[m]][, 1:p]
PD[[m]][, 1:degree] <- PD[[m]][, 1:degree] + D[[m]][, p + 1:degree]
}
PD
}
SparsePD2 <- function (xd, d) {
degree <- d - 1L
k2 <- length(xd)
xt <- c(rep.int(xd[1L], degree), xd, rep.int(xd[k2], degree))
D <- SparseD(xt, d)
Q <- NullPBC(xt, d, compact = FALSE)
PD <- vector("list", degree)
for (m in 1:degree) PD[[m]] <- D[[m]] %*% Q
m <- 1L
while (m < degree) {
Qm <- NullPBC(xt[seq.int(1L + m, length(xt) - m)], d - m, compact = FALSE)
PD[[m]] <- Matrix::crossprod(Qm, PD[[m]])
m <- m + 1L
}
PD
}
DemoKernel <- function (bw = 1) {
sig <- 0.25 * bw / qnorm(0.75)
xg <- seq.int(-4 * sig, 4 * sig, length.out = 100)
yg <- dnorm(xg, sd = sig)
op <- par(mgp = c(1.5, 0.5, 0), mar = c(2.5, 2.5, 1.5, 0.5), xpd = TRUE)
on.exit(par(op))
plot.default(xg, yg, type = "l", xlab = "u", ylab = "K(u)",
lwd = 2, yaxs = "i", bty = "l")
title("kernel functions used in stats::ksmooth")
segments(-0.5 * bw, 1 / bw, 0.5 * bw, 1 / bw, lwd = 2, lty = 3)
segments(-0.5 * bw, 0, -0.5 * bw, 1 / bw, lty = 2, col = 8)
segments(0.5 * bw, 0, 0.5 * bw, 1 / bw, lty = 2, col = 8)
legend("topleft", legend = c("box", "gaussian"), lwd = 2, lty = c(3, 1))
}
ksVec <- function (x, r = 0.1) {
n <- length(x)
i <- seq_len(n)
bw <- n * r * qnorm(0.75) / 0.75
stats::ksmooth(i, x, "normal", bw, x.points = i)$y
}
buildB <- function (x, domain = NULL, k = 10) {
if (is.null(domain)) {
a <- min(x)
b <- max(x)
domain <- c(a, b)
} else {
a <- domain[1L]
b <- domain[2L]
if (min(x) < a || max(x) > b) {
stop("'domain' does not contain all x-values!", call. = FALSE)
}
}
xd <- seq.int(a, b, length.out = k + 2)
h <- xd[2L] - xd[1L]
laux <- seq.int(to = a - h, by = h, length.out = 3L)
raux <- seq.int(from = b + h, by = h, length.out = 3L)
xt <- c(laux, xd, raux)
B <- splines::splineDesign(xt, x, ord = 4L, sparse = TRUE)
list(B = B, xt = xt, domain = domain)
}
rsplRW1 <- function (x, domain = NULL, k = 26, r = 0.2, n = 1, plot = FALSE) {
spl <- buildB(x, domain, k)
B <- spl$B
D <- SparseD(spl$xt, d = 4, m = 1, gps = TRUE)[[1]]
Coef <- PriorCoef(n, D)
b0 <- Coef$coef
i <- seq_len(k + 4)
m <- mean.default(Coef$sigma)
if (n > 1) {
bm <- b0 + rep(runif(n, -m, m), each = k + 4)
b <- matrix(0, nrow = k + 4, ncol = n)
for (j in 1:n) b[, j] <- ksVec(bm[, j], r)
y <- as_matrix(B %*% b)
} else {
bm <- b0 + runif(1, -m, m)
b <- ksVec(bm, r)
y <- (B %*% b)@x
}
if (plot) {
ylim <- range(b0, bm)
op <- par(mfrow = c(2, 2), mgp = c(1.5, 0.5, 0),
mar = c(2.5, 2.5, 1.5, 0.5), oma = c(0, 0, 1.5, 0))
on.exit(par(op))
matplot(i, b0, type = "l", col = 1:7, ylim = ylim,
xlab = "index", ylab = "B-spline coefficients")
title("(1) prior B-spline coefficients")
matplot(i, bm, "l", col = 1:7, ylim = ylim,
xlab = "index", ylab = "B-spline coefficients")
title("(2) mean-adjusted coefficients")
matplot(i, b, "l", col = 1:7, ylim = ylim,
xlab = "index", ylab = "B-spline coefficients")
title("(3) smoothed coefficients")
matplot(x, y, "l", col = 1:7, ylim = ylim, ylab = "spline(s)")
title("(4) random cubic spline(s)")
title("rw(1) prior for B-spline coefficients", outer = TRUE)
}
list(y = y, b = b, xt = spl$xt, domain = spl$domain)
}
rsplRW2 <- function (x, domain = NULL, k = 26, r = 0.15, n = 1, plot = FALSE) {
spl <- buildB(x, domain, k)
B <- spl$B
D <- SparseD(spl$xt, d = 4, m = 2, gps = TRUE)[[1]]
Coef <- PriorCoef(n, D)
b0 <- Coef$coef
i <- seq_len(k + 4)
ic <- 0.5 * (k + 5)
s <- 2 / (k + 3)
m <- mean.default(Coef$sig)
if (n > 1) {
l0 <- base::tcrossprod(i - ic, runif(n, -s, s))
bl0 <- b0 + l0
bl <- bl0 + rep(runif(n, -m, m), each = k + 4)
b <- matrix(0, nrow = k + 4, ncol = n)
for (j in 1:n) b[, j] <- ksVec(bl[, j], r)
y <- as_matrix(B %*% b)
} else {
l0 <- runif(1, -s, s) * (i - ic)
bl0 <- b0 + l0
bl <- bl0 + runif(1, -m, m)
b <- ksVec(bl, r)
y <- (B %*% b)@x
}
if (plot) {
ylim <- range(b0, bl0, bl)
op <- par(mfrow = c(2, 3), mgp = c(1.5, 0.5, 0),
mar = c(2.5, 2.5, 1.5, 0.5), oma = c(0, 0, 1.5, 0))
on.exit(par(op))
matplot(i, b0, type = "l", col = 1:7, ylim = ylim,
xlab = "index", ylab = "B-spline coefficients")
title("(1) prior B-spline coefficients")
matplot(i, l0, type = "l", col = 1:7, ylim = ylim,
xlab = "index", ylab = "B-spline coefficients")
title("(2) zero-mean random trend(s)")
matplot(i, bl0, "l", col = 1:7, ylim = ylim,
xlab = "index", ylab = "B-spline coefficients")
title("(3) trend-adjusted coefficients")
matplot(i, bl, "l", col = 1:7, ylim = ylim,
xlab = "index", ylab = "B-spline coefficients")
title("(4) mean-adjusted coefficients")
matplot(i, b, "l", col = 1:7, ylim = ylim,
xlab = "index", ylab = "B-spline coefficients")
title("(5) smoothed coefficients")
matplot(x, y, "l", col = 1:7, ylim = ylim, ylab = "spline(s)")
title("(6) random cubic spline(s)")
title("rw(2) prior for B-spline coefficients", outer = TRUE)
}
list(y = y, b = b, xt = spl$xt, domain = spl$domain)
}
NullARI <- function (rho, k = 26) {
Dt <- matrix(0, nrow = k + 4, ncol = k + 2)
diagInd <- seq.int(1, by = k + 5, length.out = k + 2)
Dt[diagInd] <- rho
Dt[diagInd + 1L] <- -(1 + rho)
Dt[diagInd + 2L] <- 1
QR <- qr.default(Dt)
E <- matrix(0, nrow = k + 4, ncol = 2)
E[c(k + 3, 2 * k + 8)] <- 1
Q <- qr.qy(QR, E)
Z1 <- rep.int(1 / sqrt(k + 4), k + 4)
g <- solve.default(Q[1:2, ], Z1[1:2])
Z2 <- c(Q %*% c(-g[2L], g[1L]))
Z2 <- VecScal(sign(Z2[1L]), Z2, overwrite = TRUE)
cbind(Z1, Z2, deparse.level = 0L)
}
rsplARI <- function (x, domain = NULL, rho = 0.7, k = 26,
r = 0.175, n = 1, plot = FALSE) {
spl <- buildB(x, domain, k)
B <- spl$B
Dt <- matrix(c(rho, -(1 + rho), 1), nrow = 3, ncol = k + 2)
D <- Matrix::t(LTB2Csparse(Dt, k = k + 4))
Coef <- PriorCoef(n, D)
b0 <- Coef$coef
Z2 <- NullARI(rho, k)[, 2]
i <- seq_len(k + 4)
s <- 1 / Z2[1L]
m <- mean.default(Coef$sigma)
if (n > 1) {
l0 <- base::tcrossprod(Z2, runif(n, -s, s))
bl0 <- b0 + l0
bl <- bl0 + rep(runif(n, -m, m), each = k + 4)
b <- matrix(0, nrow = k + 4, ncol = n)
for (j in 1:n) b[, j] <- ksVec(bl[, j], r)
y <- as_matrix(B %*% b)
} else {
l0 <- runif(1, -s, s) * Z2
bl0 <- b0 + l0
bl <- bl0 + runif(1, -m, m)
b <- ksVec(bl, r)
y <- (B %*% b)@x
}
if (plot) {
ylim <- range(b0, bl0, bl)
op <- par(mfrow = c(2, 3), mgp = c(1.5, 0.5, 0),
mar = c(2.5, 2.5, 1.5, 0.5), oma = c(0, 0, 1.5, 0))
on.exit(par(op))
matplot(i, b0, type = "l", col = 1:7, ylim = ylim,
xlab = "index", ylab = "B-spline coefficients")
title("(1) prior B-spline coefficients")
matplot(i, l0, type = "l", col = 1:7, ylim = ylim,
xlab = "index", ylab = "B-spline coefficients")
title("(2) zero-mean random trend(s)")
matplot(i, bl0, "l", col = 1:7, ylim = ylim,
xlab = "index", ylab = "B-spline coefficients")
title("(3) trend-adjusted coefficients")
matplot(i, bl, "l", col = 1:7, ylim = ylim,
xlab = "index", ylab = "B-spline coefficients")
title("(4) mean-adjusted coefficients")
matplot(i, b, "l", col = 1:7, ylim = ylim,
xlab = "index", ylab = "B-spline coefficients")
title("(5) smoothed coefficients")
matplot(x, y, "l", col = 1:7, ylim = ylim, ylab = "spline(s)")
title("(6) random cubic spline(s)")
title(sprintf("ARI(1, 1) prior for B-spline coefficients, with rho = %.2f", rho),
outer = TRUE)
}
list(y = y, b = b, xt = spl$xt, domain = spl$domain)
}
rsplm <- function (x, domain = NULL, m = 1.7, n = 1, plot = FALSE) {
if (m < 1 || m > 2) stop("m should be in [1, 2]!", call. = FALSE)
k <- 26
if (m == 1) {
rsplRW1(x, domain, k, 0.20, n, plot)
} else if (m == 2) {
rsplRW2(x, domain, k, 0.15, n, plot)
} else {
rsplARI(x, domain, m - 1, k, 0.175, n, plot)
}
}
rsplNLE <- function (x, domain = NULL, m = 1.7, n = 1000, plot = TRUE) {
y <- rsplm(x, domain, m, n)$y
k <- numeric(n)
for (j in 1:n) k[j] <- nle(y[, j])
f <- tabulate(k)
K <- 1:length(f)
N <- sum(f)
if (N < n) {
f <- c(n - N, f)
K <- c(0, K)
}
names(f) <- as.character(K)
p <- (1 / n) * f
if (plot) {
op <- par(mgp = c(1.5, 0.5, 0), mar = c(2.5, 2.5, 1.5, 0.5), xpd = TRUE)
xp <- barplot.default(p, xlab = "number of local extrema", ylab = "percentage")
text.default(xp, p, sprintf("%.2f", p), pos = 3, font = 2)
title(sprintf("m = %.2f", m))
}
p
}
rsplCMP <- function (x, domain = NULL, m = c(1, 1.7, 2), n = 1000) {
if (!is.matrix(x)) {
mchar <- sprintf("%.2f", m)
N <- length(m)
p.list <- vector("list", N)
for (l in 1:N) {
p.list[[l]] <- rsplNLE(x, domain, m[l], n, plot = FALSE) * 100
}
kl <- lapply(p.list, function (x) as.integer(names(x)))
k <- sort.int(unique.default(unlist(kl, use.names = FALSE)))
j <- lapply(kl, match, table = k)
i <- rep.int(1:N, lengths(j, use.names = FALSE))
p <- matrix(0, nrow = N, ncol = length(k), dimnames = list(mchar, as.character(k)))
p[cbind(i, unlist(j))] <- unlist(p.list)
} else {
p <- x
N <- nrow(p)
mchar <- rownames(p)
}
if (N > 3) {
warning("Can't draw more than 3 bars per group!", call. = FALSE)
} else {
my.col <- seq.int(2, by = 1, length.out = N)
op <- par(mgp = c(1.5, 0.5, 0), mar = c(2.5, 2.5, 1.5, 0.5), xpd = TRUE)
on.exit(par(op))
p <- p[, colSums(p < 0.5) < N]
xp <- barplot.default(p, beside = TRUE, col = my.col,
xlab = "number of local extrema", ylab = "percentage")
text.default(xp, p, sprintf("%.2f", p), pos = 3, col = my.col, font = 2)
legend("topright", legend = mchar, fill = my.col, text.col = my.col)
title("rsplm() with different m values")
}
p
}
rspl <- function (x, domain = NULL, n = 1) {
rsplARI(x, domain, rho = 0.7, k = 26, r = 0.175, n, plot = FALSE)
}
as_matrix <- function (A) {
if (is.matrix(A)) return(A)
sparse <- inherits(A, "dsparseMatrix")
dense <- inherits(A, "ddenseMatrix")
if (!sparse && !dense) {
stop("'A' is not a \"dsparseMatrix\" or \"ddenseMatrix\"!", call. = FALSE)
}
nnz <- length(A@x)
nr <- A@Dim[1]
nc <- A@Dim[2]
if (nnz == nr * nc) {
denseA <- matrix(A@x, nr, nc)
} else if (inherits(A, "dMatrix") && inherits(A, "CsparseMatrix")) {
i <- A@i
j <- rep.int(seq.int(0L, nc - 1L), diff.default(A@p))
denseA <- matrix(0, nr, nc)
denseA[j * nr + (i + 1L)] <- A@x
if (inherits(A, "dsCMatrix")) denseA[i * nc + (j + 1L)] <- A@x
} else {
stop("as_matrix() does not support this matrix A yet!", call. = FALSE)
}
denseA
}
VecDot <- function (x, y = NULL) {
if (is.null(y)) y <- x
.Call("C_VecDot", x, y, PACKAGE = "gps")
}
VecScal <- function (a, x, overwrite = FALSE) {
.Call("C_VecScal", a, x, overwrite, PACKAGE = "gps")
}
ChangeDim <- function (x, Value) {
.Call("C_SetDim", x, Value, PACKAGE = "gps")
}
nle <- function (x) sum(Diff(sign(Diff(x))) != 0)
Try the gps package in your browser
Any scripts or data that you put into this service are public.
gps documentation built on Nov. 2, 2023, 6:08 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Csparse2LTB <- function (A) {
if (!inherits(A, "dMatrix") || !inherits(A, "CsparseMatrix")) {
stop("'A' is not a \"dCsparseMatrix\"!")
}
n <- A@Dim[1L]
k <- A@Dim[2L]
i <- A@i
j <- rep.int(1:k, diff.default(A@p))
b1 <- sum(i == 0L)
b <- b1 - 1L
if (k < n || k > n + b) {
stop(sprintf("nrow(A) <= ncol(A) <= nrow(A) + %d expected!", b))
}
l <- j - i
if (any(l > b1 | l < 1L)) {
stop("Nonzero entries detected outside the indended band region!")
}
N1 <- sum(seq.int(n, by = -1L, length.out = b1))
N2 <- sum(seq.int(b, by = -1L, length.out = k - n))
nnz.max <- N1 + N2
x <- A@x
nnz <- length(x)
if (nnz == nnz.max) {
.Call("C_Csparse2LTB", b1, n, k, x, PACKAGE = "gps")
} else {
At <- matrix(0, b1, k)
At[cbind(l, i + 1L)] <- A@x
At
}
}
LTB2Csparse <- function (L, k = n, symmetric = FALSE) {
b1 <- nrow(L); n <- ncol(L); b <- b1 - 1L
k <- as.integer(k)
if (symmetric) {
k <- n
} else if (k < n || k > n + b) {
stop(sprintf("%d <= k <= %d expected!", n, n + b))
}
r <- n + b - k
p <- rep.int(seq.int(b1, by = -1L, length.out = r + 1L), c(n - r, rep.int(1, r)))
i <- sequence.default(nvec = p, from = seq.int(0L, n - 1L))
p <- c(0L, cumsum(p))
if (r) {
x <- numeric(p[length(p)])
x[seq_len(b1 * (k - b))] <- L[, seq_len(k - b)]
columns <- seq.int(to = n, length.out = r)
for (j in columns) {
ind <- seq.int(p[j] + 1L, p[j + 1L])
x[ind] <- L[seq_len(length(ind)), j]
}
} else {
x <- c(L)
}
if (symmetric) {
methods::new("dsCMatrix", i = i, p = p, x = x, Dim = c(k, n), uplo = "L")
} else {
methods::new("dgCMatrix", i = i, p = p, x = x, Dim = c(k, n))
}
}
LTB2Dense <- function (L, k) {
.Call("C_LTB2Dense", L, k, PACKAGE = "gps")
}
SolveLTB <- function (transA = FALSE, A, b, overwrite = FALSE) {
.Call("C_SolveLTB", transA, A, b, overwrite, PACKAGE = "gps")
}
LPBTRF <- function (A, overwrite = FALSE) {
.Call("C_LPBTRF", A, overwrite, PACKAGE = "gps")
}
DDt <- function (Dt) {
.Call("C_DDt", Dt, PACKAGE = "gps")
}
DtD <- function (Dt) {
.Call("C_DtD", Dt, PACKAGE = "gps")
}
Dt2ThinQR <- function (D, form.Q = TRUE) {
Dt <- Csparse2LTB(D)
A <- DDt(Dt)
Rt <- LPBTRF(A, overwrite = TRUE)
if (form.Q) {
denseD <- as_matrix(D)
Qt <- SolveLTB(transA = FALSE, Rt, denseD, overwrite = TRUE)
} else {
Qt <- NULL
}
list(Qt = Qt, Rt = Rt)
}
IsMonoInc <- function (x, n = length(x), xi = 1L) {
.Call("C_IsMonoInc", x, n, xi, PACKAGE = "gps") > 0L
}
MonoKnots <- function (xt, d) {
xt <- as.double(xt)
K <- length(xt)
if (K < 2 * d) stop("length(xt) >= 2 * d required!", call. = FALSE)
flag <- IsMonoInc(xt, n = K - 2 * (d - 1), xi = d)
if (!flag) stop("Domain knots are not strictly ascending!", call. = FALSE)
lbnd <- xt[seq.int(from = 1L, length.out = d)]
rbnd <- xt[seq.int(to = K, length.out = d)]
if (is.unsorted(lbnd) || is.unsorted(rbnd)) {
stop("Boundary knots are not non-decreasing!", call. = FALSE)
}
xt
}
PlaceKnots <- function (x, d, k, domain = NULL, uniform = FALSE, periodic = FALSE) {
xu <- sort.int(unique.default(x))
nx <- length(xu)
if (is.null(domain)) {
a <- xu[1L]
b <- xu[nx]
domain <- c(a, b)
} else {
a <- domain[1L]
b <- domain[2L]
if (xu[1L] < a || xu[nx] > b) {
stop("'domain' does not contain all x-values!")
}
}
degree <- d - 1
if (uniform) {
xd <- seq.int(a, b, length.out = k + 2)
h <- xd[2L] - xd[1L]
laux <- seq.int(to = a - h, by = h, length.out = degree)
raux <- seq.int(from = b + h, by = h, length.out = degree)
} else {
prob <- seq.int(0, 1, length.out = k + 2)
xd <- quantile(xu, prob, names = FALSE)
xd[c(1, k + 2)] <- domain
laux <- rep.int(a, degree)
raux <- rep.int(b, degree)
}
if (periodic) xd else c(laux, xd, raux)
}
MakeGrid <- function (xd, n, rm.dup = FALSE) {
xd <- as.double(xd)
if (!IsMonoInc(xd)) stop("'xd' is not strictly ascending!")
if (n == 1) {
lp <- xd[-length(xd)]
rp <- xd[-1]
mp <- 0.5 * (lp + rp)
return(mp)
}
if (rm.dup && (n == 2)) return(xd)
.Call("C_MakeGrid", xd, n, rm.dup, PACKAGE = "gps")
}
Zero2NA <- function (Bsparse) {
if (!inherits(Bsparse, "dgCMatrix")) {
stop("'Bsparse' is not a \"dgCMatrix\"!")
}
B <- matrix(NA_real_, Bsparse@Dim[1], Bsparse@Dim[2])
i <- Bsparse@i + 1L
j <- rep.int(seq_len(Bsparse@Dim[2]), diff.default(Bsparse@p))
B[cbind(i, j)] <- Bsparse@x
B
}
ExBS <- function (d = 4, uniform = TRUE, clamped = FALSE) {
if (d < 2 || d > 4) stop("d = 2, 3, 4 required!", call. = FALSE)
if (uniform) {
lbnd <- as.double(seq.int(to = -1, length.out = 4, by = 2))
rbnd <- as.double(seq.int(from = 1, length.out = 4, by = 2))
clamped <- FALSE
} else if (clamped) {
lbnd <- rep.int(-1, 4)
rbnd <- c(1, 2, 2.5, 4)
} else {
lbnd <- c(-3, -2.25, -1.75, -1)
rbnd <- c(1, 2, 2.5, 4)
}
xt <- c(lbnd, rbnd)
if (clamped) {
x <- MakeGrid(xt[4:8], n = 21)
} else {
x <- MakeGrid(xt, n = 21)
}
Bsparse <- splines::splineDesign(xt, x, d, outer.ok = TRUE, sparse = TRUE)
B <- Zero2NA(Bsparse)
list(xt = xt, clamped = clamped, d = d, x = x, B = B)
}
PlotExBS <- function (object) {
xt <- object$xt
d <- object$d
x <- object$x
B <- object$B
p <- ncol(B)
if (object$clamped) {
ip <- c(rep.int(NA_real_, 4 - d), apply(B[, (4 - d + 1):p], 2, which.max))
} else {
ip <- apply(B, 2, which.max)
}
xp <- x[ip]
Bp <- B[cbind(ip, 1:p)]
ymax <- max(Bp, na.rm = TRUE) * 1.2
ylim <- c(0, ymax)
nonzero <- seq.int(5 - d, length.out = d)
col.Bspl <- rep.int(8, p); col.Bspl[nonzero] <- seq_len(d)
op <- par(xpd = TRUE)
on.exit(par(op))
matplot(x, B, type = "l", lty = 1, lwd = 2, col = col.Bspl, ylim = ylim,
xaxs = "i", yaxs = "i", xaxt = "n", yaxt = "n", bty = "n", ann = FALSE)
segments(xt[1], 0, xt[8], 0)
polygon(c(-1, 1, 1, -1), c(0, 0, ymax, ymax), col = gray(0.3, 0.2), border = NA)
col.knots <- rep.int(8, 8); col.knots[seq.int(5 - d, 4 + d)] <- 1
if (object$clamped) {
stepsize <- 0.2
ystack <- seq.int(to = stepsize, length.out = 3, by = -stepsize)
points.default(xt[4:8], numeric(5), pch = 19, col = col.knots[4:8], cex = 1.5)
points.default(xt[1:3], ystack, pch = 19, col = col.knots[1:3], cex = 1.5)
} else {
points.default(xt, numeric(8), pch = 19, col = col.knots, cex = 1.5)
}
lab.Bspl <- vector("expression", p)
for (j in 1:p) {
lab.Bspl[j] <- eval(substitute(expression(B[list(j, d)]), list(j = j, d = d)))
}
text.default(xp, Bp, labels = lab.Bspl, col = col.Bspl, adj = c(0.2, 0), cex = 2)
lab.knots <- vector("expression", 8)
for (j in 1:8) {
lab.knots[j] <- eval(substitute(expression(t[j]), list(j = j)))
}
if (object$clamped) {
text.default(xt[1:3], ystack, labels = lab.knots[1:3], col = col.knots[1:3],
cex = 2, pos = 4)
text.default(xt[4:8], numeric(5), labels = lab.knots[4:8], col = col.knots[4:8],
cex = 2, pos = 1)
} else {
text.default(xt, numeric(8), labels = lab.knots, col = col.knots,
cex = 2, pos = 1)
}
usr <- par("usr")
text.default(usr[2], usr[4], labels = c(paste0("order: ", d)),
cex = 2, adj = c(1, 1))
}
DemoBS <- function (uniform = TRUE, clamped = FALSE) {
if (uniform && clamped) {
warning("uniform = TRUE implies clamped = FALSE!")
}
spl3 <- ExBS(d = 4, uniform, clamped)
spl2 <- ExBS(d = 3, uniform, clamped)
spl1 <- ExBS(d = 2, uniform, clamped)
if (uniform) {
caption <- "uniform B-splines"
} else if (clamped) {
caption <- "non-uniform & clamped B-splines"
} else {
caption <- "non-uniform B-splines"
}
op <- par(mfrow = c(3, 1), mar = c(2, 0.75, 0, 0.75), oma = c(0, 0, 2, 0))
on.exit(par(op))
PlotExBS(spl3)
PlotExBS(spl2)
PlotExBS(spl1)
title(caption, outer = TRUE, cex.main = 1.5)
}
DemoKnots <- function (aligned = TRUE) {
xt1 <- seq.int(0, 1, length.out = 12)
xt2 <- c(0, 0.08, 0.2, 0.27, 0.31, 0.42, 0.54, 0.59, 0.69, 0.82, 0.85, 1)
xt3 <- rep.int(xt2[4:9], c(4, 1, 1, 1, 1, 4))
xg1 <- MakeGrid(xt1, n = 21)
xg2 <- MakeGrid(xt2, n = 21)
xg3 <- MakeGrid(xt3[4:9], n = 21)
B1sparse <- splines::splineDesign(xt1, xg1, outer.ok = TRUE, sparse = TRUE)
B2sparse <- splines::splineDesign(xt2, xg2, outer.ok = TRUE, sparse = TRUE)
B3sparse <- splines::splineDesign(xt3, xg3, sparse = TRUE)
B1max <- max(B1sparse@x)
B2max <- max(B2sparse@x)
B3max <- max(B3sparse@x)
B1 <- Zero2NA(B1sparse)
B2 <- Zero2NA(B2sparse)
B3 <- Zero2NA(B3sparse)
Bspl.col <- c(1:6, 8, 7)
op <- par(mfrow = c(3, 1), mar = c(0.5, 0.5, 1.5, 0.5), xpd = TRUE)
on.exit(par(op))
matplot(xg1, B1, ann = FALSE, axes = FALSE, type = "l", col = Bspl.col, bty = "n",
lty = 1, lwd = 2, xlim = c(0, 1), xaxs = "i", yaxs = "i")
abline(h = 0, col = 8)
points.default(xt1, numeric(12), pch = 19, cex = 1.5)
polygon(xt1[c(4, 9, 9, 4)], c(0, 0, B1max, B1max), col = gray(0.3, 0.2), border = NA)
title("(a) uniform cubic B-splines", cex.main = 1.5)
matplot(xg2, B2, ann = FALSE, axes = FALSE, type = "l", col = Bspl.col, bty = "n",
lty = 1, lwd = 2, xlim = c(0, 1), xaxs = "i", yaxs = "i")
abline(h = 0, col = 8)
points.default(xt2, numeric(12), pch = 19, cex = 1.5)
polygon(xt2[c(4, 9, 9, 4)], c(0, 0, B2max, B2max), col = gray(0.3, 0.2), border = NA)
title("(b) non-uniform cubic B-splines", cex.main = 1.5)
xlim3 <- if (aligned) c(0, 1) else xt3[c(4, 9)]
matplot(xg3, B3, ann = FALSE, axes = FALSE, type = "l", col = Bspl.col, bty = "n",
lty = 1, lwd = 2, xlim = xlim3, xaxs = "i", yaxs = "i")
abline(h = 0, col = 8)
if (aligned) {
segments(0, 0, xt3[4], 0, col = Bspl.col[1], lwd = 2)
segments(xt3[9], 0, 1, 0, col = Bspl.col[8], lwd = 2)
segments(xt3[4], 0, xt3[4], 1, col = Bspl.col[1], lty = 2, lwd = 2)
segments(xt3[9], 0, xt3[9], 1, col = Bspl.col[8], lty = 2, lwd = 2)
}
ystack <- seq.int(0, by = 0.1, length.out = 4)
points.default(xt3, c(ystack, numeric(4), ystack), pch = 19, cex = 1.5)
polygon(xt3[c(4, 9, 9, 4)], c(0, 0, B3max, B3max), col = gray(0.3, 0.2), border = NA)
title("(c) non-uniform & clamped cubic B-splines", cex.main = 1.5)
}
PlotNull <- function (x, y, k, shape1 = 3, shape2 = 3, gps = FALSE) {
d <- 4
m <- 2
if (shape1 < 1 && shape2 < 1) {
xd <- seq.int(-1/8, 1/8, length.out = k + 2)
xd <- sign(xd) * abs(xd) ^ (1 / 3) + 0.5
distr <- "U-quadratic(0, 1)"
} else {
xd <- qbeta(seq.int(0, 1, length.out = k + 2), shape1, shape2)
distr <- sprintf("Beta(%d, %d)", shape1, shape2)
}
xt <- c(numeric(d - 1), xd, rep.int(1, d - 1))
K <- length(xt)
xg <- seq.int(0, 1, length.out = 101)
B <- splines::splineDesign(xt, x, d, sparse = TRUE)
Bg <- splines::splineDesign(xt, xg, d, sparse = TRUE)
if (gps) {
ld <- ComputeLD(xt, d)
} else {
ld <- matrix(1, k + d, d - 1)
}
H <- NullGD(ld, m, orthonormal = FALSE)
X <- as_matrix(B %*% H)
Xg <- as_matrix(Bg %*% H)
XtX <- base::crossprod(X)
Xty <- base::crossprod(X, y)
U <- chol.default(XtX)
f <- forwardsolve(U, Xty, upper.tri = TRUE, transpose = TRUE)
b <- backsolve(U, f)
yg <- c(Xg %*% b)
ylim <- range(y, yg)
plot.default(x, y, col = 8, ylim = ylim, ann = FALSE, xaxt = "n", yaxt = "n")
abline(v = xd, lty = 2, col = 8)
lines.default(xg, yg, lwd = 2)
title(distr)
}
DemoNull <- function (n, k, gps = FALSE) {
x <- seq.int(0, 1, length.out = n)
y <- rnorm(n, mean = x, sd = 0.2 * sd(x))
op <- par(mfrow = c(2, 3), mar = c(0.25, 0.25, 1.5, 0.25), oma = c(0, 0, 1.5, 0))
on.exit(par(op))
PlotNull(x, y, k, 3, 5, gps)
PlotNull(x, y, k, 5, 5, gps)
PlotNull(x, y, k, 5, 3, gps)
PlotNull(x, y, k, 1, 3, gps)
PlotNull(x, y, k, 0.5, 0.5, gps)
PlotNull(x, y, k, 3, 1, gps)
label <- sprintf("limiting %s P-spline fit at ", c("standard", "general")[gps + 1L])
expr <- substitute(expression(paste(label, lambda == +infinity)), list(label = label))
title(eval(expr), outer = TRUE)
}
DemoPBS <- function (uniform = TRUE) {
if (uniform) xd <- seq.int(0, 0.6, by = 0.1)
else xd <- c(0, 0.75, 1.25, 2, 3.5, 4.5, 5)
x <- MakeGrid(xd, 21)
PBsparse <- pbsDesign(x, xd, 4, 0, sparse = TRUE)
PB <- Zero2NA(PBsparse)
a <- xd[1L]
b <- xd[7L]
period <- b - a
raux <- xd[2:4] + period
xt <- c(xd, raux)
x.wrapped <- x[1:63] + period
B <- PB
B[1:63, 4:6] <- NA_real_
B.wrapped <- PB[1:63, 4:6]
ip <- apply(PB, 2, which.max)
xp.PBSpl <- x[ip]
yp <- PB[cbind(ip, 1:6)]
xp.PBSpl.modified <- xp.PBSpl[4:6]
sub <- xp.PBSpl.modified < xd[4]
xp.PBSpl.modified[sub] <- xp.PBSpl.modified[sub] + period
xp.BSpl <- c(xp.PBSpl[1:3], xp.PBSpl.modified)
xlim <- xt[c(1L, 10L)]
ymax <- max(yp) * 1.15
pch.knots <- rep.int(c(19, 15, 19, 15), c(1, 3, 3, 3))
col.knots <- rep.int(c(1, 8), c(7, 3))
lab1.knots <- vector("expression", 10)
lab1.knots[1] <- expression(a)
for (j in 1:5) {
lab1.knots[j + 1] <- eval(substitute(expression(s[j]), list(j = j)))
}
lab1.knots[7] <- expression(b)
for (j in 1:3) {
lab1.knots[j + 7] <- eval(substitute(expression(s[j]^','), list(j = j)))
}
lab2.knots <- vector("expression", 10)
for (j in 1:10) {
lab2.knots[j] <- eval(substitute(expression(xi[j]), list(j = j)))
}
lab.BSpl <- vector("expression", 6)
for (j in 1:6) {
lab.BSpl[j] <- eval(substitute(expression(B[j]), list(j = j)))
}
op <- par(mfrow = c(2, 1), xpd = TRUE, mar = c(3, 1, 1.25, 1))
on.exit(par(op))
matplot(x, PB, type = "n", ann = FALSE, xlim = xlim, ylim = c(0, ymax),
xaxs = "i", yaxs = "i", xaxt = "n", yaxt = "n", bty = "n")
segments(xt[1], 0, xt[10], 0, col = 8)
polygon(c(a, b, b, a), c(0, 0, ymax, ymax), col = gray(0.3, 0.2), border = NA)
matlines(x, B, lty = rep.int(c(1, 2), c(3, 3)), lwd = 2, col = 1:6)
matlines(x.wrapped, B.wrapped, lty = 2, lwd = 2, col = 4:6)
text.default(xp.BSpl, yp, labels = lab.BSpl, col = 1:6, adj = c(0.2, 0), cex = 1.5)
points.default(xt, numeric(10), pch = pch.knots, col = col.knots)
text.default(xt, numeric(10), labels = lab1.knots, col = col.knots, pos = 1, cex = 1.5)
mtext(lab2.knots, side = 1, line = 2, at = xt, col = col.knots, cex = 1.5)
title("ordinary cubic B-splines")
matplot(x, PB, type = "n", ann = FALSE, xlim = xlim, ylim = c(0, ymax),
xaxs = "i", yaxs = "i", xaxt = "n", yaxt = "n", bty = "n")
segments(a, 0, b, 0, col = 8)
polygon(c(a, b, b, a), c(0, 0, ymax, ymax), col = gray(0.3, 0.2), border = NA)
matlines(x, B, col = 1:6, lty = rep.int(c(1, 2), c(3, 3)), lwd = 2)
matlines(x.wrapped - period, B.wrapped, lty = 2, lwd = 2, col = 4:6)
text.default(xp.PBSpl, yp, labels = lab.BSpl, col = 1:6, adj = c(0.2, 0), cex = 1.5)
points.default(xd, numeric(7), pch = 19)
text.default(xd, numeric(7), labels = lab1.knots[1:7], col = col.knots[1:7],
pos = 1, cex = 1.5)
mtext(lab2.knots[1:7], side = 1, line = 2, at = xd, col = col.knots[1:7], cex = 1.5)
title("periodic cubic B-splines")
}
GetE_periodic <- function (xd, d) {
if (d < 2) stop("d >= 2 required!", call. = FALSE)
degree <- d - 1
k2 <- length(xd)
if (k2 <= d) stop("length(xd) >= d + 1 required!", call. = FALSE)
a <- rep.int(xd[1], degree)
b <- rep.int(xd[k2], degree)
local.knots.a <- c(a, xd[1:(d + 1)])
local.knots.b <- c(xd[(k2 - d):k2], b)
qs <- 0:(degree - 1)
Ca <- splines::splineDesign(local.knots.a, a, d, derivs = qs)
Cb <- splines::splineDesign(local.knots.b, b, d, derivs = qs)
Ca <- Ca[, 1:degree]
Cb <- Cb[, 2:d]
det.a <- prod(Ca[seq.int(1, length.out = degree, by = d)])
det.b <- abs(prod(Cb[seq.int(degree, length.out = degree, by = degree - 1)]))
if (det.a >= det.b) {
E <- forwardsolve(Ca, Cb)
attr(E, "eliminate") <- "head"
} else {
perm <- degree:1
E <- forwardsolve(Cb[, perm], Ca)[perm, ]
attr(E, "eliminate") <- "tail"
}
E
}
AbsorbE_periodic <- function (X, E) {
p <- ncol(X)
degree <- nrow(E)
X.head <- X[, 1:degree]
X.tail <- X[, (p - degree + 1):p]
X.mid <- X[, (degree + 1):(p - degree)]
if (attr(E, "eliminate") == "head") {
PX <- cbind(X.mid, X.tail + X.head %*% E)
} else {
PX <- cbind(X.head + X.tail %*% E, X.mid)
}
PX
}
DemoPBS2 <- function (uniform = TRUE, scaling = TRUE) {
if (uniform) xd <- seq.int(0, 0.6, by = 0.1)
else xd <- c(0, 0.75, 1.25, 2, 3.5, 4.5, 5)
a <- xd[1]
b <- xd[7]
Aknots <- c(rep.int(a, 3), xd, rep.int(b, 3))
x <- gps::MakeGrid(xd, n = 21)
B <- splines::splineDesign(Aknots, x)
E <- GetE_periodic(xd, 4)
PB <- AbsorbE_periodic(B, E)
lab.knots <- vector("expression", 7)
for (j in 1:7) {
lab.knots[j] <- eval(substitute(expression(xi[j]), list(j = j)))
}
lab.BSpl <- vector("expression", 9)
for (j in 1:9) {
lab.BSpl[j] <- eval(substitute(expression(B[j]), list(j = j)))
}
lab.PBSpl <- vector("expression", 9)
for (j in 1:6) {
lab.PBSpl[j] <- eval(substitute(expression(tilde(B)[j]), list(j = j)))
}
ip.BSpl <- apply(B, 2, which.max)
xp.BSpl <- x[ip.BSpl]
yp.BSpl <- B[cbind(ip.BSpl, 1:9)]
ip.PBSpl <- apply(abs(PB), 2, which.max)
xp.PBSpl <- x[ip.PBSpl]
yp.PBSpl <- PB[cbind(ip.PBSpl, 1:6)]
positive <- yp.PBSpl > 0
negative <- yp.PBSpl < 0
if (scaling) {
big <- abs(yp.PBSpl) > 1
big.yp.PBspl <- yp.PBSpl[big]
fctr <- sign(big.yp.PBspl) / max(abs(big.yp.PBspl))
PB[, big] <- PB[, big] * rep(fctr, each = length(x))
yp.PBSpl[big] <- yp.PBSpl[big] * fctr
}
B[B == 0] <- NA
PB[PB == 0] <- NA
op <- par(mfrow = c(2, 1), xpd = TRUE, mar = c(2, 0.75, 1, 0.75))
on.exit(par(op))
matplot(x, B, type = "l", lty = rep(c(2, 1, 2), each = 3), lwd = 2, ann = FALSE,
xaxs = "i", yaxs = "i", xaxt = "n", yaxt = "n", bty = "n", ylim = c(0, 1.1))
segments(a, 0, b, 0, col = 8)
text.default(xp.BSpl, yp.BSpl, labels = lab.BSpl, col = 1:6, adj = c(0.5, 0), cex = 1.5)
points.default(xd, numeric(7), pch = 19)
ystack <- seq.int(0.1, by = 0.1, length.out = 3)
points.default(rep.int(c(a, b), c(3, 3)), c(ystack, ystack), pch = 19, col = 8)
text.default(xd, numeric(7), labels = lab.knots, pos = 1, cex = 1.5)
title("ordinary cubic B-splines")
if (attr(E, "eliminate") == "head") {
lty.PBSpl <- rep(1:2, each = 3)
} else {
lty.PBSpl <- rep(2:1, each = 3)
}
ylim.PBSpl <- range(-0.03, yp.PBSpl) * c(1, 1.1)
matplot(x, PB, type = "l", lty = lty.PBSpl, lwd = 2, ann = FALSE, ylim = ylim.PBSpl,
xaxs = "i", yaxs = "i", xaxt = "n", yaxt = "n", bty = "n")
segments(a, 0, b, 0, col = 8)
text.default(xp.PBSpl[positive], yp.PBSpl[positive], labels = lab.PBSpl[positive],
col = (1:6)[positive], adj = c(0.5, 0), cex = 1.5)
text.default(xp.PBSpl[negative], yp.PBSpl[negative], labels = lab.PBSpl[negative],
col = (1:6)[negative], pos = 1, cex = 1.5)
points.default(xd, numeric(7), pch = 19)
text.default(xd, numeric(7), labels = lab.knots, pos = 1, cex = 1.5)
title("periodic cubic B-splines")
}
DemoSpl <- function (uniform = TRUE) {
if (uniform) {
xd <- 1:6
xt <- seq.int(-2, 9)
xg <- MakeGrid(xt, 21)
b <- c(0.44, 1.11, 1.66, 0.25, 1.6, 1.43, 1.49, 2.52)
} else {
xd <- c(1, 2.06, 2.65, 3.22, 3.91, 6)
xt <- xd[c(1, 1, 1, 1, 2:5, 6, 6, 6, 6)]
xg <- MakeGrid(xd, 21)
b <- c(1.9, 1.84, 1.35, 1.67, 0.48, 1.85, 1.59, 1.35)
}
Bg <- splines::splineDesign(xt, xg, outer.ok = uniform, sparse = TRUE)
Bgdense <- Zero2NA(Bg)
yg <- (Bg %*% b)@x
if (uniform) {
i1 <- 21 * 3 + 1
i2 <- 21 * 8
x <- xg[i1:i2]
y <- yg[i1:i2]
} else {
x <- xg
y <- yg
}
ymax <- max(y)
yd <- y[c(1, 1:5 * 21)]
op <- par(xpd = TRUE, mar = c(2, 2, 1.5, 0.5))
on.exit(par(op))
plot.default(x, y, type = "n", xlim = c(-2, 9), ylim = c(0, ymax),
xaxs = "i", yaxs = "i", ann = FALSE, bty = "n")
polygon(c(1, 6, 6, 1), c(0, 0, ymax, ymax), col = gray(0.3, 0.2), border = NA)
lines.default(x, y, lwd = 2)
points.default(xd, yd, pch = c(17, 19, 19, 19, 19, 17))
matlines(xg, Bgdense, type = "l", lty = 1, lwd = 2, col = c(1:6, 8, 7))
points.default(xd[2:5], numeric(4), pch = 19)
if (uniform) {
points.default(c(-2:0, 7:9), numeric(6), pch = 15)
} else {
segments(-2, 0, 1, 0, lwd = 2, col = 1)
segments(6, 0, 9, 0, lwd = 2, col = 7)
segments(1, 1, 1, 0, lwd = 2, col = 1, lty = 2)
segments(6, 1, 6, 0, lwd = 2, col = 7, lty = 2)
ystack <- seq.int(0.1, by = 0.1, length.out = 3)
points.default(c(1, 1, 1, 6, 6, 6), c(ystack, ystack), pch = 15)
}
points.default(c(1, 6), numeric(2), pch = 17)
type <- c("non-uniform", "uniform")[uniform + 1L]
title(sprintf("cubic spline with %s knots", type))
dB <- splines::splineDesign(xt, rep.int(xd[1:5], 3), derivs = rep(1:3, each = 5))
pc <- c(dB %*% b) * rep(1 / factorial(1:3), each = 5)
lab1 <- paste0("f", 1:5)
lab2 <- c("constant", "linear", "quadratic", "cubic")
pc <- matrix(c(yd[1:5], pc), ncol = 4, dimnames = list(lab1, lab2))
list(xd = xd, bspl.coef = b, piecepoly.coef = pc)
}
CheckSupp <- function (x, xt, d) {
degree <- d - 1
K <- length(xt)
xd <- xt[seq.int(d, K - degree)]
nI <- length(xd) - 1L
id <- findInterval(x, xd, rightmost.closed = TRUE)
if (min(id) < 1L || max(id) > nI) {
stop("Domain does not contain all x-values!", call. = FALSE)
}
ni <- tabulate(id, nbins = nI)
.Call("C_CheckSupp", ni, d, PACKAGE = "gps") > 0L
}
gps2bspl <- function (x, xt, d, m, gps = TRUE, periodic = FALSE,
unique.x.provided = FALSE) {
if (d < 2) stop("d >= 2 required!", call. = FALSE)
if (m < 1 || m > d - 1) stop("1 <= m <= d - 1 required!", call. = FALSE)
xt <- MonoKnots(xt, d)
p <- length(xt) - d
xu <- if (unique.x.provided) x else unique.default(x)
if (p > length(xu)) {
stop("There are more B-splines than unique x-values!", call. = FALSE)
}
if (CheckSupp(xu, xt, d)) {
stop("There are no x-values on some B-spline's support!", call. = FALSE)
}
B <- splines::splineDesign(xt, x, d, sparse = TRUE)
D <- SparseD(xt, d, m, gps)[[1L]]
Dt <- Csparse2LTB(D)
S <- DtD(Dt)
if (ncol(S) > p) S <- S[, 1:p]
if (periodic) {
C <- GetPBC(xt, d, compact = FALSE, sparse = FALSE)
} else {
C <- matrix(0, nrow = 0, ncol = p)
}
list(B = B, Dt = Dt, S = S, xt = xt, d = d, m = m, C = C)
}
gps2red <- function (bspl, y, w = NULL, scalePen = TRUE) {
y <- as.double(y)
X <- bspl$B
n <- X@Dim[1L]
if (length(y) != n) stop("Incorrect number of y-values!")
if (!is.null(w)) {
if (length(w) != n) stop("Incorrect number of weights!")
w <- (1 / sum(w)) * w
w <- sqrt(w)
X@x <- w[X@i + 1L] * X@x
y <- w * y
}
if (scalePen) {
scalePen.fctr <- VecDot(X@x) / VecDot(bspl$Dt)
Dt <- VecScal(sqrt(scalePen.fctr), bspl$Dt, overwrite = FALSE)
S <- VecScal(scalePen.fctr, bspl$S, overwrite = FALSE)
} else {
scalePen.fctr <- 1
Dt <- bspl$Dt
S <- bspl$S
}
Z <- Matrix::crossprod(X)
z <- Matrix::crossprod(X, y)@x
Z <- Csparse2LTB(Z)
L <- try(LPBTRF(Z, overwrite = FALSE), silent = TRUE)
if (inherits(L, "try-error")) {
stop("B-spline design matrix does not have full column rank!")
}
f <- SolveLTB(transA = FALSE, L, z, overwrite = FALSE)
minRSS <- VecDot(y) - VecDot(f)
E <- .Call("C_FormE", L, Dt, PACKAGE = "gps")
list(n = n, Z = Z, z = z, L = L, f = f, minRSS = minRSS,
scalePen.fctr = scalePen.fctr, Dt = Dt, S = S, E = E)
}
gps2DR <- function (red, tol = 1e-6) {
t0 <- Sys.time()
q <- ncol(red$E)
dbar <- .Call("C_MeanDR", red$E, PACKAGE = "gps")
dq <- .Call("C_MinDR", red$E, tol, PACKAGE = "gps")
d1 <- .Call("C_MaxDR", red$L, red$Dt, tol, PACKAGE = "gps")
dq.cutoff <- d1 * (.Machine$double.eps / 2)
if (dq < dq.cutoff) {
dq <- dq.cutoff
warning("The smallest eigenvalue is erratic as E'E is numerically singular!")
}
approx.DR <- .Call("C_ApproxDR", q, dq, d1, dbar, tol, PACKAGE = "gps")
t1 <- Sys.time()
secs <- difftime(t1, t0, units = "secs")[[1L]]
list(mean = dbar, min = dq, max = d1, approx.values = approx.DR, secs = secs)
}
ExactDR <- function (E) {
t0 <- Sys.time()
A <- .Call("C_LAUUM", E, PACKAGE = "gps")
exact.DR <- base::eigen(A, symmetric = TRUE, only.values = TRUE)$values
q <- length(exact.DR)
d1 <- exact.DR[1L]
dq.cutoff <- d1 * (.Machine$double.eps / 2)
ind <- which(exact.DR < dq.cutoff)
if (length(ind)) {
exact.DR[ind] <- dq.cutoff
warning("Small eigenvalues are erratic as E'E is numerically singular!")
}
t1 <- Sys.time()
secs <- difftime(t1, t0, units = "secs")[[1L]]
list(values = exact.DR, secs = secs)
}
Rho2REDF <- function (DR, rho) {
.Call("C_Rho2REDF", DR, rho, PACKAGE = "gps")
}
REDF2Rho <- function (DR, redf, rho0, MaxNewton, tol = 1e-6) {
.Call("C_REDF2Rho", DR, redf, rho0, MaxNewton, tol, PACKAGE = "gps")
}
gps2RhoLim <- function (DR, kappa = 0.01) {
lam.min <- kappa / (1 - kappa) / DR$mean
rho.min <- log(lam.min)
lam.max <- (1 - kappa) / kappa / DR$min
rho.max <- log(lam.max)
approx.DR <- DR$approx.values
if (is.na(approx.DR[1L])) {
rho.max.heuristic <- NA_real_
} else {
redf <- length(approx.DR) * kappa
rho.mid <- 0.5 * (rho.min + rho.max)
MaxNewton <- 0.25 * (rho.max - rho.min)
rho.max.heuristic <- REDF2Rho(approx.DR, redf, rho.mid, MaxNewton)
}
list(min = rho.min, max = rho.max, max.heuristic = rho.max.heuristic)
}
Q1fun <- function (z, a = 0, b = 1) {
if (a >= b) stop("a < b required!")
h <- z * z - z
theta <- exp(a + (b - a) * z)
bounds <- c(0, b - a)
list(h = h, theta = theta, bounds = bounds)
}
Q2fun <- function (z, a = 0, b = 1) {
if (a >= b) stop("a < b required!")
y <- 1 - z
B0 <- y * y * y
B1 <- 3 * z * y * y
B2 <- 3 * z * z * y
B3 <- z * z * z
h <- B1 - B2
theta <- exp((B0 + B2) * a + (B2 + B3) * b)
bounds <- c(a, (2 * a + b) / 3)
list(h = h, theta = theta, bounds = bounds)
}
map01 <- function (x) (1 / (x[1L] - x[length(x)])) * (x - x[length(x)])
ApproxDR <- function (DR, Qfun, gamma.grid = seq.int(0, 1, 0.05), verbose = TRUE) {
if (!is.function(Qfun)) stop("'Qfun' is not a function!")
Qz <- Qfun(z = 0.5, a = 0, b = 1)
if (anyNA(match(c("h", "theta", "bounds"), names(Qz)))) {
stop("'Qfun' needs to return elements 'h', 'theta' and 'bounds'!")
}
q <- length(DR)
a <- log(DR[q])
b <- log(DR[1L])
dbar <- sum(DR) / q
logDR <- map01(log(DR))
step <- 1 / (q + 1)
p <- seq.int(from = step, by = step, length.out = q)
aggregated.DR <- numeric(q)
success <- 0
op <- par(mfrow = c(1, 2), mar = c(2, 2, 1.5, 0.5))
on.exit(par(op))
for (gamma in gamma.grid) {
z <- map01(log(1 - p) - gamma * log(p))
zlab <- sprintf("z(%.2f)", gamma)
Qz <- Qfun(z, a, b)
plot.default(p, logDR, type = "l", ann = FALSE)
lines.default(p, z, lty = 2)
title(sprintf("log(d) and %s", zlab))
plot.default(z, logDR, type = "l", ann = FALSE)
title(sprintf("log(d) ~ %s", zlab))
theta <- Qz$theta
h <- Qz$h
eta <- theta * h
bounds <- Qz$bounds
tmpDR <- .Call("C_RootApproxDR", theta, h, eta, dbar, bounds, PACKAGE = "gps")
if (tmpDR[1L]) {
success <- success + 1
aggregated.DR <- aggregated.DR + tmpDR
lines.default(z, map01(log(tmpDR)), lty = 2)
}
if (verbose) readline("Hit <Enter> to continue: ")
if (success && tmpDR[1L] == 0) break
}
list(success = success, aggregated.DR = aggregated.DR, gamma.exit = gamma)
}
DemoApproxDR <- function (DR, Qfun = NULL, verbose = TRUE) {
gamma.grid <- seq.int(0, 1, 0.05)
if (is.null(Qfun)) {
out1 <- ApproxDR(DR, Q1fun, gamma.grid, verbose)
gamma.exit <- out1$gamma.exit
if (gamma.exit < 1) {
gamma.grid <- seq.int(gamma.exit, 1, 0.05)
} else {
gamma.grid <- numeric(0)
}
out2 <- ApproxDR(DR, Q2fun, gamma.grid, verbose)
success <- out1$success + out2$success
approx.DR <- (out1$aggregated.DR + out2$aggregated.DR) / success
} else {
out <- ApproxDR(DR, Qfun, gamma.grid, verbose)
success <- out$success
approx.DR <- out$aggregated.DR / success
}
if (success) {
op <- par(mfrow = c(1, 2), mar = c(2, 2, 1.5, 0.5))
on.exit(par(op))
q <- length(DR)
step <- 1 / (q + 1)
p <- seq.int(from = step, by = step, length.out = q)
logDR <- map01(log(DR))
result <- list(min = DR[q], mean = sum(DR) / q, approx.values = approx.DR)
rho.lim <- gps2RhoLim(result)
rho.grid <- seq.int(rho.lim$min, rho.lim$max, length.out = 100)
redf <- Rho2REDF(DR, rho.grid)
approx.redf <- Rho2REDF(approx.DR, rho.grid)
redf.heuristic <- Rho2REDF(DR, rho.lim$max.heuristic)
mapping <- sprintf("loose: %.2f -> %.2f\nheuristic: %.2f -> %.2f",
rho.lim$max, redf[100L], rho.lim$max.heuristic, redf.heuristic)
plot.default(p, logDR, type = "l", ann = FALSE)
lines.default(p, map01(log(approx.DR)), lty = 2)
title("log(d)")
plot.default(rho.grid, redf, type = "l", ann = FALSE)
lines.default(rho.grid, approx.redf, lty = 2)
points.default(rho.lim$max.heuristic, redf.heuristic, pch = 19)
text.default(rho.lim$max, redf[1L], labels = mapping, adj = c(1, 1))
title("redf ~ rho")
} else {
warning("Unable to approximate Demmler-Reinsch eigenvalues!")
}
}
GridPWLS <- function (bspl, red, rho) {
.Call("C_GridPWLS", red$Z, red$L, red$S, red$z, t.default(bspl$C), rho,
PACKAGE = "gps")
}
GridGCV <- function (red, pwls) {
.Call("C_GridGCV", red$n, red$L, red$f, red$minRSS, pwls$beta, pwls$edf,
PACKAGE = "gps")
}
gps2GS <- function (x, y, w = NULL, xt, d = 4, m = 2, gps = TRUE,
periodic = FALSE, ng = 20, scalePen = TRUE) {
bspl <- gps2bspl(x, xt, d, m, gps, periodic)
red <- gps2red(bspl, y, w, scalePen)
DR <- gps2DR(red)
rho.lim <- gps2RhoLim(DR)
if (ng > 0) {
if (is.na(rho.lim$max.heuristic)) {
rho <- seq.int(rho.lim$min, rho.lim$max, length.out = ng)
} else {
rho <- seq.int(rho.lim$min, rho.lim$max.heuristic, length.out = ng)
}
t0 <- Sys.time()
pwls <- GridPWLS(bspl, red, rho)
pwls$gcv <- GridGCV(red, pwls)
t1 <- Sys.time()
pwls$secs <- difftime(t1, t0, units = "secs")[[1L]]
} else {
pwls <- NULL
}
S <- LTB2Csparse(bspl$S, symmetric = TRUE)
BWB <- LTB2Csparse(red$Z, symmetric = TRUE)
eqn <- list(B = bspl$B, BWB = BWB, BWy = red$z,
omega = red$scalePen.fctr, S = S, C = bspl$C)
list(eqn = eqn, eigen = DR, rho.lim = rho.lim, E = red$E, pwls = pwls)
}
DemoRhoLim <- function (fit, plot = TRUE) {
DR <- ExactDR(fit$E)
approx.DR <- fit$eigen[c("approx.values", "secs")]
names(approx.DR) <- c("values", "secs")
q <- length(DR$values)
rho.min.star <- fit$rho.lim$min
rho.max.star <- fit$rho.lim$max
redf.min <- 0.01 * q
redf.max <- 0.99 * q
rho.mid <- (rho.min.star + rho.max.star) / 2
MaxNewton <- 0.25 * (rho.max.star - rho.min.star)
rho.min <- REDF2Rho(DR$values, redf.max, rho.mid, MaxNewton)
rho.max <- REDF2Rho(DR$values, redf.min, rho.mid, MaxNewton)
ng <- 100L
rho.grid <- seq.int(rho.min.star, rho.max.star, length.out = ng)
redf.exact <- Rho2REDF(DR$values, rho.grid)
redf.min.star <- redf.exact[ng]
redf.max.star <- redf.exact[1L]
rho.max.heuristic <- fit$rho.lim$max.heuristic
if (is.na(rho.max.heuristic)) {
redf.approx <- rep.int(NA_real_, ng)
redf.min.heuristic <- NA_real_
warning("No approximated eigenvalues and heuristic upper bound for 'rho'!")
} else {
redf.approx <- Rho2REDF(approx.DR$values, rho.grid)
redf.min.heuristic <- Rho2REDF(DR$values, rho.max.heuristic)
}
map.exact <- sprintf("[%.1f, %.1f] -> [%.1f, %.1f]",
rho.min, rho.max, redf.min, redf.max)
map.wider <- sprintf("[%.1f, %.1f] -> [%.1f, %.1f]",
rho.min.star, rho.max.star, redf.min.star, redf.max.star)
map.heuristic <- sprintf("%.1f -> %.1f", rho.max.heuristic, redf.min.heuristic)
if (plot) {
op <- par(mfrow = c(1, 2), mar = c(2, 2, 1.5, 0.5))
on.exit(par(op))
p <- (1:q) / (q + 1)
plot.default(p, log(DR$values), type = "l", ann = FALSE)
lines.default(p, log(approx.DR$values), lty = 2)
text.default(p[q], log(DR$values[1L]), adj = c(1, 1),
labels = sprintf("exact: solid\napprox: dashed"))
title("log eigenvalues")
plot.default(rho.grid, redf.exact, type = "l", ann = FALSE)
lines.default(rho.grid, redf.approx, lty = 2)
segments(rho.min, q, rho.min, 0.75 * q, col = 8, lty = 2)
segments(rho.max, 0, rho.max, 0.25 * q, col = 8, lty = 2)
points.default(rho.max.heuristic, redf.min.heuristic, pch = 19)
text.default(rho.max.star, redf.max.star, adj = c(1, 1),
labels = sprintf("exact: %s\nwider: %s\nheuristic: %s",
map.exact, map.wider, map.heuristic))
title("redf ~ rho")
}
eigen <- list(approx = approx.DR, exact = DR)
limit <- list(rho = list(exact = c(rho.min, rho.max),
wider = c(rho.min.star, rho.max.star),
heuristic = rho.max.heuristic),
redf = list(exact = c(redf.min, redf.max),
wider = c(redf.min.star, redf.max.star),
heuristic = redf.min.heuristic))
mapping <- list(exact = map.exact, wider = map.wider, heuristic = map.heuristic)
redf <- list(rho = rho.grid, exact = redf.exact, approx = redf.approx)
list(eigen = eigen, limit = limit, mapping = mapping, redf = redf)
}
Diff <- function (x, k = 1L, n = length(x), xi = 1L) {
.Call("C_Diff", x, k, n, xi, PACKAGE = "gps")
}
SparseDelta <- function (r) {
r <- as.integer(r)
x <- rep.int(c(-1, 1), r)
i <- rep(seq.int(0L, r - 1L), each = 2L)
p <- c(0L, seq.int(1L, length.out = r, by = 2L), 2L * r)
methods::new("dgCMatrix", i = i, p = p, Dim = c(r, r + 1L), x = x)
}
SparseWtDelta <- function (h) {
r <- length(h)
x <- rep.int(c(-1, 1), r) * rep(1 / h, each = 2)
i <- rep(seq.int(0L, r - 1L), each = 2L)
p <- c(0L, seq.int(1L, length.out = r, by = 2L), 2L * r)
methods::new("dgCMatrix", i = i, p = p, Dim = c(r, r + 1L), x = x)
}
SparseGD <- function (xt, d, m) {
K <- length(xt)
D <- vector("list", m)
h <- Diff(x = xt, k = d - 1L, n = K - 2L, xi = 2L)
D[[1L]] <- SparseWtDelta(h)
i <- 2L
while (i <= m) {
h <- Diff(x = xt, k = d - i, n = K - 2L * i, xi = i + 1L)
D[[i]] <- SparseWtDelta(h) %*% D[[i - 1L]]
i <- i + 1L
}
D
}
DiffCoef <- function (b, xt, d, m) {
if (d < 2) stop("d >= 2 required!", call. = FALSE)
if (m < 1 || m >= d) stop("1 <= m <= d - 1 required!", call. = FALSE)
xt <- MonoKnots(xt, d)
K <- length(xt)
if (length(b) != K - d) {
stop("length(b) == length(xt) - d required!", call. = FALSE)
}
b <- as.double(b)
for (i in 1:m) {
h <- Diff(x = xt, k = d - i, n = K - 2L * i, xi = i + 1L)
b <- Diff(b) / h
}
b
}
ComputeLD <- function (xt, d) {
.Call("C_ComputeLD", xt, d, PACKAGE = "gps")
}
NullGD <- function (ld, m = 1, orthonormal = TRUE) {
basis <- .Call("C_NullGD", ld, m, PACKAGE = "gps")
if (orthonormal && (m > 1)) {
Q <- qr.Q(qr.default(basis[, m:1]))[, m:1]
i <- sequence.default(1:(m - 1))
j <- rep.int(2:m, 1:(m - 1))
Q[cbind(i, j)] <- 0
basis <- Q
}
basis
}
QuadWts <- function (ord) {
if (ord == 1) {
return(2)
}
p <- seq.int(0, ord - 1)
P <- outer(seq.int(-1, 1, length.out = ord), p, "^")
Pinv <- solve.default(P)
pow <- outer(1:ord, p, "+")
H <- (1 - (-1) ^ pow) / pow
base::crossprod(Pinv, H %*% Pinv)
}
SbarBlocks <- function (xt, d, m) {
if (d < 2) stop("d >= 2 required!", call. = FALSE)
if (m < 0 || m >= d) stop("0 <= m <= d - 1 required!", call. = FALSE)
xt <- MonoKnots(xt, d)
K <- length(xt)
ord <- d - m
if (ord == 1) {
h <- Diff(x = xt, n = K - 2 * (d - 1), xi = d)
return(h)
}
W <- QuadWts(ord)
xd <- xt[seq.int(d, K - d + 1)]
xg <- MakeGrid(xd, ord)
xt.local <- xt[seq.int(1 + m, K - m)]
B <- splines::splineDesign(xt.local, xg, ord, sparse = TRUE)
.Call("C_SbarBlocks", xd, W, B@x, PACKAGE = "gps")
}
GramBS <- function (xt, d) {
blocks <- SbarBlocks(xt, d, m = 0)
LTB <- .Call("C_SbarLTB", blocks, PACKAGE = "gps")
LTB2Csparse(LTB, symmetric = TRUE)
}
btSb <- function (b, xt, d, m) {
db <- DiffCoef(b, xt, d, m)
blocks <- SbarBlocks(xt, d, m)
.Call("C_btSb", blocks, db, PACKAGE = "gps")
}
SparseD <- function (xt, d, m = NULL, gps = TRUE) {
d <- as.integer(d)
if (d < 2L) stop("d >= 2 required!", call. = FALSE)
if (is.null(m)) {
m <- seq_len(d - 1L)
} else {
m <- sort.int(as.integer(m), method = "radix")
}
m.min <- m[1L]
m.max <- m[length(m)]
if (m.min < 1L || m.max >= d) {
stop("1 <= m <= d - 1 required!", call. = FALSE)
}
NAMES <- sprintf("order.%d", m)
xt <- MonoKnots(xt, d)
D <- SparseGD(xt, d, m.max)
D <- D[m]
names(D) <- NAMES
if (gps) return(D)
nm <- length(m)
Sbar <- vector("list", nm); names(Sbar) <- NAMES
K <- vector("list", nm); names(K) <- NAMES
for (i in 1:nm) {
m_i <- m[i]
D_i <- D[[i]]
blocks <- SbarBlocks(xt, d, m_i)
if (is.array(blocks)) {
Sbar.LTB <- .Call("C_SbarLTB", blocks, PACKAGE = "gps")
Sbar.Matrix <- LTB2Csparse(Sbar.LTB, symmetric = TRUE)
L.LTB <- LPBTRF(Sbar.LTB, overwrite = TRUE)
L.Matrix <- LTB2Csparse(L.LTB)
K_i <- Matrix::crossprod(L.Matrix, D_i)
} else {
k1 <- length(blocks)
Sbar.Matrix <- methods::new("ddiMatrix", diag = "N", x = blocks, Dim = c(k1, k1))
L.diag <- sqrt(blocks)
newx <- L.diag[D_i@i + 1L] * D_i@x
K_i <- methods::new("dgCMatrix", i = D_i@i, p = D_i@p, Dim = D_i@Dim, x = newx)
}
Sbar[[i]] <- Sbar.Matrix
K[[i]] <- K_i
}
sandwich <- list(D = D, Sbar = Sbar)
structure(K, sandwich = sandwich)
}
PriorCoef1 <- function (n, D) {
q <- D@Dim[1L]
S <- as_matrix(Matrix::crossprod(D))
ei <- base::eigen(S, symmetric = TRUE)
qs <- seq_len(q)
d <- ei$values[qs]
U <- ei$vectors[, qs]
di <- 1 / d
e <- rnorm(n * q)
if (n > 1) e <- ChangeDim(e, c(q, n))
b <- U %*% (sqrt(di) * e)
if (n == 1) b <- c(b)
b
}
PriorCoef2 <- function (n, D) {
Dim <- D@Dim
q <- Dim[1L]
p <- Dim[2L]
m <- p - q
e <- rnorm(n * q)
if (n > 1) e <- ChangeDim(e, c(q, n))
Dt <- Csparse2LTB(D)
y <- SolveLTB(transA = TRUE, Dt, e, overwrite = TRUE)
R <- as_matrix(D[, seq.int(q + 1, p), drop = FALSE])
X <- SolveLTB(transA = TRUE, Dt, R, overwrite = TRUE)
XtX <- base::crossprod(X)
Xty <- base::crossprod(X, y)
U <- chol.default(XtX + diag(m))
f <- forwardsolve(U, Xty, upper.tri = TRUE, transpose = TRUE)
b <- backsolve(U, f)
r <- y - X %*% b
if (n > 1) rbind(r, b) else c(r, b)
}
PriorCoef3 <- function (n, D) {
Dim <- D@Dim
q <- Dim[1L]
p <- Dim[2L]
QR <- Dt2ThinQR(D, form.Q = FALSE)
Rt <- QR$Rt
e <- rnorm(n * q)
if (n > 1) e <- ChangeDim(e, c(q, n))
x <- SolveLTB(transA = FALSE, Rt, e, overwrite = TRUE)
x <- SolveLTB(transA = TRUE, Rt, x, overwrite = TRUE)
b <- Matrix::crossprod(D, x)
if (n == 1) b@x else as_matrix(b)
}
PriorCoef <- function (n, D) {
b <- PriorCoef3(n, D)
sig <- sqrt(MPinv0(D, only.diag = TRUE))
scale.fctr <- 1 / sig[1]
b <- VecScal(scale.fctr, b, overwrite = TRUE)
sig <- VecScal(scale.fctr, sig, overwrite = TRUE)
list(coef = b, sigma = sig)
}
MPinvUUt <- function (D, method = "qr") {
if (method == "qr") {
QR <- Dt2ThinQR(D, form.Q = TRUE)
Qt <- QR$Qt
Rt <- QR$Rt
SolveLTB(transA = TRUE, Rt, Qt, overwrite = TRUE)
} else if (method == "eigen") {
q <- D@Dim[1L]
S <- as_matrix(Matrix::crossprod(D))
ei <- base::eigen(S, symmetric = TRUE)
qs <- seq_len(q)
d <- ei$values[qs]
Q <- ei$vectors[, qs]
di <- 1 / d
sqrt(di) * t.default(Q)
} else {
stop("'method' must be either \"qr\" or \"eigen\"!", call. = FALSE)
}
}
MPinv0 <- function (D, only.diag = FALSE) {
Ut <- MPinvUUt(D, method = "qr")
if (only.diag) {
base::colSums(Ut * Ut)
} else {
base::crossprod(Ut)
}
}
MPinv <- function (D, only.diag = FALSE) {
V <- MPinv0(D, only.diag)
VecScal(1 / V[1L], V, overwrite = TRUE)
}
Prior2Cond <- function (k) {
rho <- seq.int(0.05, 1, 0.05)
N <- length(rho)
cond.X <- numeric(N)
for (i in 1:N) {
Dt <- matrix(c(rho[i], -(1 + rho[i]), 1), nrow = 3, ncol = k + 2)
R <- matrix(0, nrow = k + 2, ncol = 2L)
R[c(k + 1, k + 2, 2 * k + 4)] <- c(1, -(1 + rho[i]), 1)
X <- SolveLTB(transA = TRUE, Dt, R, overwrite = TRUE)
XtX <- base::crossprod(X)
ka <- try(kappa(X, exact = TRUE), silent = TRUE)
cond.X[i] <- if (inherits(ka, "try-error")) Inf else ka
}
cond.XtX <- cond.X ^ 2
op <- par(mfrow = c(1, 2), mgp = c(1.5, 0.5, 0), mar = c(2.5, 2.5, 1.5, 0.5))
on.exit(op)
plot.default(rho, log10(cond.X), type = "b", xlab = expression(rho),
ylab = "log10(condition number)", main = "X")
plot.default(rho, log10(cond.XtX), type = "b", xlab = expression(rho),
ylab = "log10(condition number)", main = "X'X")
list(rho = rho, X = cond.X, XtX = cond.XtX)
}
GetPBC <- function (xt, d, compact = TRUE, transpose = FALSE, sparse = TRUE) {
d <- as.integer(d)
if (d < 2L) stop("d >= 2 required!")
degree <- d - 1L
K <- length(xt)
if (length(xt) < degree + 2L * d) {
stop("length(xt) >= 3 * d - 1 required!")
}
p <- K - d
a <- xt[d]
b <- xt[K - degree]
nDeriv <- seq.int(0L, length.out = degree)
knots.a <- xt[seq.int(from = 1L, length.out = 2L * d)]
knots.b <- xt[seq.int(to = K, length.out = 2L * d)]
Ca <- splines::splineDesign(knots.a, rep.int(a, degree), d, nDeriv)
Cb <- splines::splineDesign(knots.b, rep.int(b, degree), d, nDeriv)
Ca <- Ca[, seq.int(1L, degree)]
Cb <- Cb[, seq.int(2L, d)]
C.compact <- cbind(Ca, -Cb, deparse.level = 0L)
Ct.compact <- t.default(C.compact)
if (compact) {
if (transpose) Ct.compact else C.compact
} else {
if (sparse) {
if (transpose) {
Ct.i <- rep.int(c(nDeriv, seq.int(to = p - 1L, length.out = degree)), degree)
Ct.p <- seq.int(0L, by = 2L * degree, length.out = d)
ChangeDim(Ct.compact, NULL)
Ct.x <- Ct.compact
methods::new("dgCMatrix", i = Ct.i, p = Ct.p, Dim = c(p, degree), x = Ct.x)
} else {
blocksize <- degree * degree
C.i <- rep.int(nDeriv, 2L * degree)
C.p <- c(seq.int(0, by = degree, length.out = d),
rep.int(blocksize, p - 2L * degree),
seq.int(blocksize + degree, by = degree, length.out = degree))
ChangeDim(C.compact, NULL)
C.x <- C.compact
methods::new("dgCMatrix", i = C.i, p = C.p, Dim = c(degree, p), x = C.x)
}
} else {
ind <- c(seq.int(1L, degree), seq.int(to = p, length.out = degree))
if (transpose) {
Ct <- matrix(0, p, degree)
Ct[ind, ] <- Ct.compact
Ct
} else {
C <- matrix(0, degree, p)
C[, ind] <- C.compact
C
}
}
}
}
NullPBC <- function (xt, d, compact = TRUE) {
d <- as.integer(d)
degree <- d - 1L
Ct <- GetPBC(xt, d, compact = TRUE, transpose = TRUE)
QR <- qr.default(Ct)
Q <- qr.Q(QR, complete = TRUE)
Q.compact <- Q[, seq.int(d, 2L * degree), drop = FALSE]
if (compact) {
Q.compact
} else {
p <- length(xt) - d
pp <- p - degree
N <- pp - degree
Q.compact.i <- c(seq.int(0L, length.out = degree),
seq.int(to = p - 1L, length.out = degree))
Q.i <- c(seq.int(degree, length.out = N), rep.int(Q.compact.i, degree))
Q.p <- c(seq.int(0, length.out = N),
seq.int(N, by = 2L * degree, length.out = d))
Q.x <- c(rep.int(1, N), Q.compact)
methods::new("dgCMatrix", i = Q.i, p = Q.p, Dim = c(p, pp), x = Q.x)
}
}
pbsDesign <- function (x, xd, d, nDeriv = 0, sparse = FALSE, wrap = TRUE) {
d <- as.integer(d)
if (d < 2L) stop("d >= 2 required!")
degree <- d - 1L
if (nDeriv < 0 || nDeriv > degree) stop("0 <= nDeriv <= d - 1 required!")
xd <- as.double(xd)
if (!IsMonoInc(xd)) stop("'xd' is not strictly ascending!")
k2 <- length(xd)
if (k2 <= d) stop("length(xd) >= d + 1 required!")
a <- xd[1L]
b <- xd[k2]
if (min(x) < a || max(x) > b) stop("domain does not contain all x-values!")
if (wrap) {
pbsDesign1(x, xd, d, nDeriv, sparse)
} else {
pbsDesign2(x, xd, d, nDeriv, sparse)
}
}
pbsDesign1 <- function (x, xd, d, nDeriv = 0, sparse = FALSE) {
degree <- d - 1L
k2 <- length(xd)
p <- k2 - 1L
a <- xd[1L]
b <- xd[k2]
period <- b - a
raux <- xd[2:d] + period
xt <- c(xd, raux)
b.aux <- raux[degree]
ind.x.wrapped <- which(x < xd[d])
if (length(ind.x.wrapped) == 0L) {
PB <- splines::splineDesign(xt, x, d, nDeriv, sparse = sparse)
} else {
nx <- length(x)
x.wrapped <- x[ind.x.wrapped] + period
x.wrapped[x.wrapped < b] <- b
x.wrapped[x.wrapped > b.aux] <- b.aux
x.padded <- c(x, x.wrapped)
xt0 <- c(rep.int(a, degree), xt, rep.int(b.aux, degree))
B <- splines::splineDesign(xt0, x.padded, d, nDeriv, sparse = TRUE)
B.p <- B@p[seq.int(d, length(B@p) - degree)]
ind <- seq.int(B.p[1L] + 1L, B.p[p - degree + 1L])
unchanged.i <- B@i[ind]
unchanged.x <- B@x[ind]
B.p.vital <- B.p[seq.int(to = p + 1L, length.out = d)]
ind <- seq.int(B.p.vital[1L] + 1L, B.p.vital[d])
unordered.i <- B@i[ind]
unordered.x <- B@x[ind]
ind <- (unordered.i >= nx)
unordered.i[ind] <- ind.x.wrapped[unordered.i[ind] + 1L - nx] - 1L
unordered.j <- rep.int(1:degree, diff.default(B.p.vital))
ind <- order(unordered.j, unordered.i, method = "radix")
reordered.i <- unordered.i[ind]
reordered.x <- unordered.x[ind]
PB.p <- B.p - B.p[1L]
PB.i <- c(unchanged.i, reordered.i)
PB.x <- c(unchanged.x, reordered.x)
if (sparse) {
PB <- methods::new("dgCMatrix", i = PB.i, p = PB.p, Dim = c(nx, p), x = PB.x)
} else {
PB.j <- rep.int(1:p, diff.default(PB.p))
PB <- matrix(0, nx, p)
PB[cbind(PB.i + 1L, PB.j)] <- PB.x
}
}
PB
}
pbsDesign2 <- function (x, xd, d, nDeriv = 0, sparse = FALSE) {
degree <- d - 1L
k2 <- length(xd)
xt <- c(rep.int(xd[1L], degree), xd, rep.int(xd[k2], degree))
B <- splines::splineDesign(xt, x, d, nDeriv, sparse = TRUE)
p <- ncol(B)
ind <- seq.int(B@p[d] + 1L, B@p[p + 1L - degree])
PB.i <- B@i[ind]
PB.x <- B@x[ind]
PB.p <- B@p[seq.int(d, p + 1L)]
PB.p <- PB.p - PB.p[1L]
IND <- seq.int(to = k2, length.out = degree)
PB.p[IND] <- PB.p[k2 - degree]
bnd_cols <- c(seq.int(1L, degree), seq.int(k2, p))
nnz_per_col <- B@p[bnd_cols + 1L] - B@p[bnd_cols]
active_col <- nnz_per_col > 0L
bnd_cols <- bnd_cols[active_col]
num_active_cols <- length(bnd_cols)
if (num_active_cols > 0L) {
Q <- NullPBC(xt, d)
Q <- Q[active_col, , drop = FALSE]
i_lst <- x_lst <- vector("list", num_active_cols)
for (i in 1:num_active_cols) {
j <- bnd_cols[i]
ind <- seq.int(B@p[j] + 1L, B@p[j + 1L])
i_lst[[i]] <- B@i[ind]
x_lst[[i]] <- B@x[ind]
}
active_row <- sort.int(unique.default(unlist(i_lst)), method = "radix")
num_active_rows <- length(active_row)
block.i <- unlist(lapply(i_lst, match, active_row))
block.j <- rep.int(1:num_active_cols, nnz_per_col[active_col])
block.x <- unlist(x_lst)
B.boundary <- matrix(0, num_active_rows, num_active_cols)
B.boundary[cbind(block.i, block.j)] <- block.x
B.transformed <- B.boundary %*% Q
PB.i <- c(PB.i, rep.int(active_row, degree))
PB.x <- c(PB.x, B.transformed)
PB.p[IND] <- PB.p[IND] + num_active_rows * (1:degree)
}
nr <- length(x)
nc <- k2 - 1L
if (sparse) {
PB <- methods::new("dgCMatrix", i = PB.i, p = PB.p, Dim = c(nr, nc), x = PB.x)
} else {
PB.j <- rep.int(1:nc, diff.default(PB.p))
PB <- matrix(0, nr, nc)
PB[cbind(PB.i + 1L, PB.j)] <- PB.x
}
PB
}
SparsePD <- function (xd, d, wrap = TRUE) {
d <- as.integer(d)
if (d < 2L) stop("d >= 2 required!")
degree <- d - 1L
xd <- as.double(xd)
if (!IsMonoInc(xd)) stop("'xd' is not strictly ascending!")
k2 <- length(xd)
if (k2 <= d) stop("length(xd) >= d + 1 required!")
if (wrap) {
SparsePD1(xd, d)
} else {
SparsePD2(xd, d)
}
}
SparsePD1 <- function (xd, d) {
degree <- d - 1L
k2 <- length(xd)
p <- k2 - 1L
a <- xd[1L]
b <- xd[k2]
period <- b - a
laux <- xd[(k2 - degree):p] - period
raux <- xd[2:d] + period
xt <- c(laux, xd, raux)
D <- SparseD(xt, d)
PD <- vector("list", degree)
for (m in 1:degree) {
PD[[m]] <- D[[m]][, 1:p]
PD[[m]][, 1:degree] <- PD[[m]][, 1:degree] + D[[m]][, p + 1:degree]
}
PD
}
SparsePD2 <- function (xd, d) {
degree <- d - 1L
k2 <- length(xd)
xt <- c(rep.int(xd[1L], degree), xd, rep.int(xd[k2], degree))
D <- SparseD(xt, d)
Q <- NullPBC(xt, d, compact = FALSE)
PD <- vector("list", degree)
for (m in 1:degree) PD[[m]] <- D[[m]] %*% Q
m <- 1L
while (m < degree) {
Qm <- NullPBC(xt[seq.int(1L + m, length(xt) - m)], d - m, compact = FALSE)
PD[[m]] <- Matrix::crossprod(Qm, PD[[m]])
m <- m + 1L
}
PD
}
DemoKernel <- function (bw = 1) {
sig <- 0.25 * bw / qnorm(0.75)
xg <- seq.int(-4 * sig, 4 * sig, length.out = 100)
yg <- dnorm(xg, sd = sig)
op <- par(mgp = c(1.5, 0.5, 0), mar = c(2.5, 2.5, 1.5, 0.5), xpd = TRUE)
on.exit(par(op))
plot.default(xg, yg, type = "l", xlab = "u", ylab = "K(u)",
lwd = 2, yaxs = "i", bty = "l")
title("kernel functions used in stats::ksmooth")
segments(-0.5 * bw, 1 / bw, 0.5 * bw, 1 / bw, lwd = 2, lty = 3)
segments(-0.5 * bw, 0, -0.5 * bw, 1 / bw, lty = 2, col = 8)
segments(0.5 * bw, 0, 0.5 * bw, 1 / bw, lty = 2, col = 8)
legend("topleft", legend = c("box", "gaussian"), lwd = 2, lty = c(3, 1))
}
ksVec <- function (x, r = 0.1) {
n <- length(x)
i <- seq_len(n)
bw <- n * r * qnorm(0.75) / 0.75
stats::ksmooth(i, x, "normal", bw, x.points = i)$y
}
buildB <- function (x, domain = NULL, k = 10) {
if (is.null(domain)) {
a <- min(x)
b <- max(x)
domain <- c(a, b)
} else {
a <- domain[1L]
b <- domain[2L]
if (min(x) < a || max(x) > b) {
stop("'domain' does not contain all x-values!", call. = FALSE)
}
}
xd <- seq.int(a, b, length.out = k + 2)
h <- xd[2L] - xd[1L]
laux <- seq.int(to = a - h, by = h, length.out = 3L)
raux <- seq.int(from = b + h, by = h, length.out = 3L)
xt <- c(laux, xd, raux)
B <- splines::splineDesign(xt, x, ord = 4L, sparse = TRUE)
list(B = B, xt = xt, domain = domain)
}
rsplRW1 <- function (x, domain = NULL, k = 26, r = 0.2, n = 1, plot = FALSE) {
spl <- buildB(x, domain, k)
B <- spl$B
D <- SparseD(spl$xt, d = 4, m = 1, gps = TRUE)[[1]]
Coef <- PriorCoef(n, D)
b0 <- Coef$coef
i <- seq_len(k + 4)
m <- mean.default(Coef$sigma)
if (n > 1) {
bm <- b0 + rep(runif(n, -m, m), each = k + 4)
b <- matrix(0, nrow = k + 4, ncol = n)
for (j in 1:n) b[, j] <- ksVec(bm[, j], r)
y <- as_matrix(B %*% b)
} else {
bm <- b0 + runif(1, -m, m)
b <- ksVec(bm, r)
y <- (B %*% b)@x
}
if (plot) {
ylim <- range(b0, bm)
op <- par(mfrow = c(2, 2), mgp = c(1.5, 0.5, 0),
mar = c(2.5, 2.5, 1.5, 0.5), oma = c(0, 0, 1.5, 0))
on.exit(par(op))
matplot(i, b0, type = "l", col = 1:7, ylim = ylim,
xlab = "index", ylab = "B-spline coefficients")
title("(1) prior B-spline coefficients")
matplot(i, bm, "l", col = 1:7, ylim = ylim,
xlab = "index", ylab = "B-spline coefficients")
title("(2) mean-adjusted coefficients")
matplot(i, b, "l", col = 1:7, ylim = ylim,
xlab = "index", ylab = "B-spline coefficients")
title("(3) smoothed coefficients")
matplot(x, y, "l", col = 1:7, ylim = ylim, ylab = "spline(s)")
title("(4) random cubic spline(s)")
title("rw(1) prior for B-spline coefficients", outer = TRUE)
}
list(y = y, b = b, xt = spl$xt, domain = spl$domain)
}
rsplRW2 <- function (x, domain = NULL, k = 26, r = 0.15, n = 1, plot = FALSE) {
spl <- buildB(x, domain, k)
B <- spl$B
D <- SparseD(spl$xt, d = 4, m = 2, gps = TRUE)[[1]]
Coef <- PriorCoef(n, D)
b0 <- Coef$coef
i <- seq_len(k + 4)
ic <- 0.5 * (k + 5)
s <- 2 / (k + 3)
m <- mean.default(Coef$sig)
if (n > 1) {
l0 <- base::tcrossprod(i - ic, runif(n, -s, s))
bl0 <- b0 + l0
bl <- bl0 + rep(runif(n, -m, m), each = k + 4)
b <- matrix(0, nrow = k + 4, ncol = n)
for (j in 1:n) b[, j] <- ksVec(bl[, j], r)
y <- as_matrix(B %*% b)
} else {
l0 <- runif(1, -s, s) * (i - ic)
bl0 <- b0 + l0
bl <- bl0 + runif(1, -m, m)
b <- ksVec(bl, r)
y <- (B %*% b)@x
}
if (plot) {
ylim <- range(b0, bl0, bl)
op <- par(mfrow = c(2, 3), mgp = c(1.5, 0.5, 0),
mar = c(2.5, 2.5, 1.5, 0.5), oma = c(0, 0, 1.5, 0))
on.exit(par(op))
matplot(i, b0, type = "l", col = 1:7, ylim = ylim,
xlab = "index", ylab = "B-spline coefficients")
title("(1) prior B-spline coefficients")
matplot(i, l0, type = "l", col = 1:7, ylim = ylim,
xlab = "index", ylab = "B-spline coefficients")
title("(2) zero-mean random trend(s)")
matplot(i, bl0, "l", col = 1:7, ylim = ylim,
xlab = "index", ylab = "B-spline coefficients")
title("(3) trend-adjusted coefficients")
matplot(i, bl, "l", col = 1:7, ylim = ylim,
xlab = "index", ylab = "B-spline coefficients")
title("(4) mean-adjusted coefficients")
matplot(i, b, "l", col = 1:7, ylim = ylim,
xlab = "index", ylab = "B-spline coefficients")
title("(5) smoothed coefficients")
matplot(x, y, "l", col = 1:7, ylim = ylim, ylab = "spline(s)")
title("(6) random cubic spline(s)")
title("rw(2) prior for B-spline coefficients", outer = TRUE)
}
list(y = y, b = b, xt = spl$xt, domain = spl$domain)
}
NullARI <- function (rho, k = 26) {
Dt <- matrix(0, nrow = k + 4, ncol = k + 2)
diagInd <- seq.int(1, by = k + 5, length.out = k + 2)
Dt[diagInd] <- rho
Dt[diagInd + 1L] <- -(1 + rho)
Dt[diagInd + 2L] <- 1
QR <- qr.default(Dt)
E <- matrix(0, nrow = k + 4, ncol = 2)
E[c(k + 3, 2 * k + 8)] <- 1
Q <- qr.qy(QR, E)
Z1 <- rep.int(1 / sqrt(k + 4), k + 4)
g <- solve.default(Q[1:2, ], Z1[1:2])
Z2 <- c(Q %*% c(-g[2L], g[1L]))
Z2 <- VecScal(sign(Z2[1L]), Z2, overwrite = TRUE)
cbind(Z1, Z2, deparse.level = 0L)
}
rsplARI <- function (x, domain = NULL, rho = 0.7, k = 26,
r = 0.175, n = 1, plot = FALSE) {
spl <- buildB(x, domain, k)
B <- spl$B
Dt <- matrix(c(rho, -(1 + rho), 1), nrow = 3, ncol = k + 2)
D <- Matrix::t(LTB2Csparse(Dt, k = k + 4))
Coef <- PriorCoef(n, D)
b0 <- Coef$coef
Z2 <- NullARI(rho, k)[, 2]
i <- seq_len(k + 4)
s <- 1 / Z2[1L]
m <- mean.default(Coef$sigma)
if (n > 1) {
l0 <- base::tcrossprod(Z2, runif(n, -s, s))
bl0 <- b0 + l0
bl <- bl0 + rep(runif(n, -m, m), each = k + 4)
b <- matrix(0, nrow = k + 4, ncol = n)
for (j in 1:n) b[, j] <- ksVec(bl[, j], r)
y <- as_matrix(B %*% b)
} else {
l0 <- runif(1, -s, s) * Z2
bl0 <- b0 + l0
bl <- bl0 + runif(1, -m, m)
b <- ksVec(bl, r)
y <- (B %*% b)@x
}
if (plot) {
ylim <- range(b0, bl0, bl)
op <- par(mfrow = c(2, 3), mgp = c(1.5, 0.5, 0),
mar = c(2.5, 2.5, 1.5, 0.5), oma = c(0, 0, 1.5, 0))
on.exit(par(op))
matplot(i, b0, type = "l", col = 1:7, ylim = ylim,
xlab = "index", ylab = "B-spline coefficients")
title("(1) prior B-spline coefficients")
matplot(i, l0, type = "l", col = 1:7, ylim = ylim,
xlab = "index", ylab = "B-spline coefficients")
title("(2) zero-mean random trend(s)")
matplot(i, bl0, "l", col = 1:7, ylim = ylim,
xlab = "index", ylab = "B-spline coefficients")
title("(3) trend-adjusted coefficients")
matplot(i, bl, "l", col = 1:7, ylim = ylim,
xlab = "index", ylab = "B-spline coefficients")
title("(4) mean-adjusted coefficients")
matplot(i, b, "l", col = 1:7, ylim = ylim,
xlab = "index", ylab = "B-spline coefficients")
title("(5) smoothed coefficients")
matplot(x, y, "l", col = 1:7, ylim = ylim, ylab = "spline(s)")
title("(6) random cubic spline(s)")
title(sprintf("ARI(1, 1) prior for B-spline coefficients, with rho = %.2f", rho),
outer = TRUE)
}
list(y = y, b = b, xt = spl$xt, domain = spl$domain)
}
rsplm <- function (x, domain = NULL, m = 1.7, n = 1, plot = FALSE) {
if (m < 1 || m > 2) stop("m should be in [1, 2]!", call. = FALSE)
k <- 26
if (m == 1) {
rsplRW1(x, domain, k, 0.20, n, plot)
} else if (m == 2) {
rsplRW2(x, domain, k, 0.15, n, plot)
} else {
rsplARI(x, domain, m - 1, k, 0.175, n, plot)
}
}
rsplNLE <- function (x, domain = NULL, m = 1.7, n = 1000, plot = TRUE) {
y <- rsplm(x, domain, m, n)$y
k <- numeric(n)
for (j in 1:n) k[j] <- nle(y[, j])
f <- tabulate(k)
K <- 1:length(f)
N <- sum(f)
if (N < n) {
f <- c(n - N, f)
K <- c(0, K)
}
names(f) <- as.character(K)
p <- (1 / n) * f
if (plot) {
op <- par(mgp = c(1.5, 0.5, 0), mar = c(2.5, 2.5, 1.5, 0.5), xpd = TRUE)
xp <- barplot.default(p, xlab = "number of local extrema", ylab = "percentage")
text.default(xp, p, sprintf("%.2f", p), pos = 3, font = 2)
title(sprintf("m = %.2f", m))
}
p
}
rsplCMP <- function (x, domain = NULL, m = c(1, 1.7, 2), n = 1000) {
if (!is.matrix(x)) {
mchar <- sprintf("%.2f", m)
N <- length(m)
p.list <- vector("list", N)
for (l in 1:N) {
p.list[[l]] <- rsplNLE(x, domain, m[l], n, plot = FALSE) * 100
}
kl <- lapply(p.list, function (x) as.integer(names(x)))
k <- sort.int(unique.default(unlist(kl, use.names = FALSE)))
j <- lapply(kl, match, table = k)
i <- rep.int(1:N, lengths(j, use.names = FALSE))
p <- matrix(0, nrow = N, ncol = length(k), dimnames = list(mchar, as.character(k)))
p[cbind(i, unlist(j))] <- unlist(p.list)
} else {
p <- x
N <- nrow(p)
mchar <- rownames(p)
}
if (N > 3) {
warning("Can't draw more than 3 bars per group!", call. = FALSE)
} else {
my.col <- seq.int(2, by = 1, length.out = N)
op <- par(mgp = c(1.5, 0.5, 0), mar = c(2.5, 2.5, 1.5, 0.5), xpd = TRUE)
on.exit(par(op))
p <- p[, colSums(p < 0.5) < N]
xp <- barplot.default(p, beside = TRUE, col = my.col,
xlab = "number of local extrema", ylab = "percentage")
text.default(xp, p, sprintf("%.2f", p), pos = 3, col = my.col, font = 2)
legend("topright", legend = mchar, fill = my.col, text.col = my.col)
title("rsplm() with different m values")
}
p
}
rspl <- function (x, domain = NULL, n = 1) {
rsplARI(x, domain, rho = 0.7, k = 26, r = 0.175, n, plot = FALSE)
}
as_matrix <- function (A) {
if (is.matrix(A)) return(A)
sparse <- inherits(A, "dsparseMatrix")
dense <- inherits(A, "ddenseMatrix")
if (!sparse && !dense) {
stop("'A' is not a \"dsparseMatrix\" or \"ddenseMatrix\"!", call. = FALSE)
}
nnz <- length(A@x)
nr <- A@Dim[1]
nc <- A@Dim[2]
if (nnz == nr * nc) {
denseA <- matrix(A@x, nr, nc)
} else if (inherits(A, "dMatrix") && inherits(A, "CsparseMatrix")) {
i <- A@i
j <- rep.int(seq.int(0L, nc - 1L), diff.default(A@p))
denseA <- matrix(0, nr, nc)
denseA[j * nr + (i + 1L)] <- A@x
if (inherits(A, "dsCMatrix")) denseA[i * nc + (j + 1L)] <- A@x
} else {
stop("as_matrix() does not support this matrix A yet!", call. = FALSE)
}
denseA
}
VecDot <- function (x, y = NULL) {
if (is.null(y)) y <- x
.Call("C_VecDot", x, y, PACKAGE = "gps")
}
VecScal <- function (a, x, overwrite = FALSE) {
.Call("C_VecScal", a, x, overwrite, PACKAGE = "gps")
}
ChangeDim <- function (x, Value) {
.Call("C_SetDim", x, Value, PACKAGE = "gps")
}
nle <- function (x) sum(Diff(sign(Diff(x))) != 0)
Try the gps package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.