penutils: Utility functions for working with wiggliness penalties

In gps: General P-Splines

Utility functions for working with wiggliness penalties

Description

Evaluate \|\bm{D\beta}\|^2 without using \bm{D}.

Usage

DiffCoef(b, xt, d, m)

btSb(b, xt, d, m)

Arguments

b	a vector of B-spline coefficients (length(b) == length(xt) - d).
xt	full knot sequence for ordinary B-splines (length(xt) >= 2 * d).
d	B-spline order (d \ge 2).
m	penalty order (1 \le m \le d - 1).

Details

Implicit Evaluation of the Penalty

Sometimes we want to evaluate the penalty \|\bm{D\beta}\|^2 for some \bm{\beta}. The obvious way is to do the matrix-vector multiplication \bm{D\beta} then compute its L_2 norm, however, implicit evaluation without using \bm{D} is possible. For general P-splines, we can calculate \bm{D}_{\textrm{gps}}\bm{\beta} by taking order-m general differences between elements of \bm{\beta}, and function DiffCoef does this. For O-splines, the evaluation can be more refined. Denote domain knots by s_0,\ s_1,\ s_2,\ \ldots,\ s_k,\ s_{k + 1}, where (s_j)_1^k are interior knots and s_0 = a, s_{k + 1} = b are domain endpoints. The derivative penalty adds up local wiggliness measure on each interval: \int_a^b f^{(m)}(x)^2\mathrm{d}x = \sum_{j = 0}^k\int_{s_j}^{s_{j + 1}} f^{(m)}(x)^2\mathrm{d}x. Function btSb calculates each integral in the summation and returns those additive components in a vector.

Value

DiffCoef (for general P-splines only) returns \bm{D}_{\textrm{gps}}\bm{\beta} as a vector.

btSb (for O-splines only) returns a vector with element \int_{s_j}^{s_{j + 1}} f^{(m)}(x)^2\mathrm{d}x.

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))

gps documentation built on Nov. 2, 2023, 6:08 p.m.