penutils {gps} | R Documentation |
Utility functions for working with wiggliness penalties
Description
Evaluate without using
.
Usage
DiffCoef(b, xt, d, m)
btSb(b, xt, d, m)
Arguments
b |
a vector of B-spline coefficients ( |
xt |
full knot sequence for ordinary B-splines ( |
d |
B-spline order ( |
m |
penalty order ( |
Details
Implicit Evaluation of the Penalty
Sometimes we want to evaluate the penalty for some
. The obvious way is to do the matrix-vector multiplication
then compute its
norm, however, implicit evaluation without using
is possible. For general P-splines, we can calculate
by taking order-
general differences between elements of
, and function
DiffCoef
does this. For O-splines, the evaluation can be more refined. Denote domain knots by , where
are interior knots and
,
are domain endpoints. The derivative penalty adds up local wiggliness measure on each interval:
. Function
btSb
calculates each integral in the summation and returns those additive components in a vector.
Value
DiffCoef
(for general P-splines only) returns as a vector.
btSb
(for O-splines only) returns a vector with element .
Examples
require(Matrix)
require(gps)
## 11 domain knots at equal quantiles of Beta(3, 3) distribution
xd <- qbeta(seq.int(0, 1, by = 0.1), 3, 3)
## full knots (with clamped boundary knots) for constructing cubic B-splines
xt <- c(0, 0, 0, xd, 1, 1, 1)
## compute 2nd order D matrix for O-splines
D.os <- SparseD(xt, d = 4, m = 2, gps = FALSE)
D2.os <- D.os$order.2
## get 2nd order D matrix for general P-splines
## we can of course compute it with D.gps <- SparseD(xt, d = 4, m = 2, gps = TRUE)
## but it is readily stored in the "sandwich" attribute of 'D.os'
D.gps <- attr(D.os, "sandwich")$D
D2.gps <- D.gps$order.2
## random B-spline coefficients
b <- rnorm(ncol(D2.gps))
## two ways to evaluate a difference penalty
diff.b1 <- DiffCoef(b, xt, d = 4, m = 2) ## implicit
diff.b2 <- as.numeric(D2.gps %*% b) ## explicit
range(diff.b1 - diff.b2) / max(abs(diff.b1))
## several ways to evaluate a derivative penalty
sum(btSb(b, xt, d = 4, m = 2)) ## recommended
sum(as.numeric(D2.os %*% b) ^ 2)
S2.os <- crossprod(D2.os); sum(b * as.numeric(S2.os %*% b))