slp.mgcv: Generate a Basis Matrix for Discrete Prolate Spheroidal (Slepian) Sequences

Description

Generate the basis matrix for a particular N, W Slepian sequence family member, with the additional property that the smoother captures/passes constants without distortion. Can be quite slow in execution due to the latter property.

Based on smooth.construct.cr.smooth.spec for implementation with mgcv.

Usage

"smooth.construct"(object,data,knots)

Arguments

object

a smooth specification object, usually generated by a model term s(..., bs = 'slp', ..., xt = list(...)), and for this type, requiring an additional xt = list() object containing parameters. For examples, see below.

data

a list containing just the data required by this term, with names corresponding to object[['term']]. Typically just a single time index array.

knots

a list containing any knots supplied for basis setup -- should be NULL, as slp basis objects are not generated from knots.

Value

C: a constraint matrix which restricts mgcv from modifying the (already) orthogonal and mean-passing/capturing basis vectors.
K: the user-specified number of basis vectors, or the computed K from user-supplied W.
W: the user-specified bandwidth W, or the computed W from user-supplied K.
fullBasis: the full-span computed, normalized basis set, before contiguity is taken into account. Used by predict when given an object of type slp.smooth.
contiguous: a logical variable declaring whether or not the input time array was considered to be contiguous by the basis computation procedure.
wx: the ``corrected'' input time array; if contiguous == FALSE then this will be the same as data[object[['term']]].

Details

slp is based on .dpss, which generates a family of Discrete Prolate Spheroidal (Slepian) Sequences. These vectors are orthonormal, have alternating even/odd parity, and form the optimally concentrated basis set for the subspace of R^N corresponding to the bandwidth W. Full details are given in Slepian (1978). These basis functions have natural boundary conditions, and lack any form of knot structure. This version is returned for naive = TRUE.

The dpss basis vectors can be adapted to provide the additional useful property of capturing or passing constants perfectly. That is, the smoother matrix S formed from the returned rectangular matrix will either reproduce constants at near round-off precision, i.e., S %*% rep(1, N) = rep(1, N), for naive = FALSE with intercept = TRUE, or will pass constants, i.e., S %*% rep(1, N) = rep(0, N), for naive = FALSE with intercept = FALSE.

The primary use is in modeling formula to directly specify a Slepian time-based smoothing term in a model: see the examples.

For large N this routine can be very slow. If you are computing models with large N, we highly recommend pre-computing the basis object, then using it in your models without recomputation. The third example below demonstrates this approach.

References

Wood S.N. (2006) Generalized Additive Models: An Introduction with R. Chapman and Hall/CRC Press.

Hastie T.J. & Tibshirani, R.J. (1990) Generalized Additive Models. Chapman and Hall.

Thomson, D.J (1982) Spectrum estimation and harmonic analysis. Proceedings of the IEEE. Volume 70, number 9, pp. 1055-1096.

Slepian, David (1978) Prolate Spheroidal Wave Functions, Fourier Analysis, and Uncertainty V: the Discrete Case. Bell System Technical Journal. Volume 57, pp. 1371-1429.

Examples

Run this code

    # Examples using pkg:mgcv
    library("mgcv")
    library("slp")

    N <- 730
    W <- 8 / N
    K <- 16         # will actually use 15 df as intercept = FALSE
    x <- rnorm(n = N, sd = 1)
    y <- x + rnorm(n = N, sd = 2) + 5.0
    t <- seq(1, N)    

    # note: all three examples share identical results

    # example with in-call computation, using K (df)
    fit1 <- gam(y ~ x + s(t, bs = 'slp', xt = list(K = K)), family = gaussian) 

    # example with in-call computation, using W
    fit2 <- gam(y ~ x + s(t, bs = 'slp', xt = list(W = W)), family = gaussian)

Run the code above in your browser using DataLab