# R/slp-gam.R In slp: Discrete Prolate Spheroidal (Slepian) Sequence Regression Smoothers

#### Documented in slp

```################################################################################
#
#  Extension for gam package allowing a projection Slepian basis to be used;
#  part of the slp package
#  (C) wburr, July 2014
#
#
################################################################################

slp <- function(x, W = NA, K = NA, deltat = 1, naive = FALSE,
intercept = FALSE, customSVD = TRUE, forceC = FALSE,
returnS = FALSE) {

# logical checks (note: W is assumed to be in the same units as deltat, i.e., days)
stopifnot(is.numeric(x), (is.na(W) | is.numeric(W) & W > 0 & W < 1/(2*deltat)), deltat %in% c(1, 6),
(is.numeric(K) & K <= length(x) | is.na(K)), is.logical(customSVD))

namesx <- names(x)
x <- as.vector(x)

# x should be a contiguous time array; if not, convert to one, back-convert at the end
if(is.na(max(x[-1L] - x[-length(x)])) | max(x[-1L] - x[-length(x)]) > deltat) {
warning("slp: Input time array is not contiguous.")
minT <- min(x, na.rm = TRUE); maxT <- max(x, na.rm = TRUE)
stopifnot(is.numeric(minT), is.numeric(maxT))
wx <- seq(minT, maxT, deltat)

if(deltat == 6) {
wx <- seq(minT, maxT, deltat)
if(wx[length(wx)] != maxT) { stop("slp: Input time array not properly time aligned to 6-day samples per deltat = 6.") }
}

cat(" Input Time Array: \n")
cat(paste(" *        samples: ", length(x), "\n", sep = ""))
cat(paste(" *      max - min: ", maxT - minT + 1, "\n", sep = ""))
cat(paste(" *         deltat: ", deltat, "\n", sep = ""))
} else {
wx <- x
}

N <- length(wx)

if(is.na(W) & is.na(K)) { stop("Must set one of K or W for family selection.") }
if(!is.na(K) & !is.na(W)) {
if(K > ceiling(2 * N * W)) { stop("Using more than 2NW basis vectors is not recommended (or supported).") }
} else {
if(is.na(K)) {
K <- as.integer(round(2 * N * W))
} else {   # Case of is.na(W), K is set
if(floor(K) != K) {  # K doesn't have to be integer _class_, but it needs to be an integer
K <- floor(K)
warning(paste("slp: K choice not integer. Truncated to K = ", K, sep = ""))
}
W <- as.numeric((K)) / as.numeric((2 * length(wx)))
}
}

# Logical check: if N, K and W match one of the saved objects, simply
# return that object; otherwise, run through the generation process
Wn <- round(W * 365.2425)    # convert to df/year
if(checkSaved(N, Wn, K) & !forceC & !returnS) {

data(list = as.character(paste0("basis_N_", N, "_W_", Wn, "_K_", K)),
envir = environment())

if(!intercept) { basis <- basis[, -1] }

# need to convert back to original time array -- the NAs in the original
if(any(is.na(x))) {
basis[which(!(wx %in% x[!is.na(x)])), ] <- rep(NA, ncol(basis))
}
} else {

# start by generating baseline Slepian basis vectors
v <- .dpss(n = N, nw = N * W, k = K)
if(naive) {        # Case of non-mean-adjusted Slepians, SLP
basis <- v
if(returnS) {
sMat <- v %*% t(v)
}
} else {           # Case of mean-adjusted Slepians, SLP2 or SLP3

# Equations follow from Thomson (2001), see help file for full citation
alpha <- t(rep(1, N)) %*% v %*% t(v) %*% rep(1, N)
R <- rep(1/sqrt(alpha), N)
U <- t(v) %*% R

#  Two options: intercept = TRUE / FALSE
sRaw <- v %*% (diag(K) - U %*% t(U)) %*% t(v)

if(!returnS) {
if(customSVD) {
basis <- .slpsvd(A = sRaw, N = N, K = K)\$u      # optimized LAPACK SVD for our edge case
} else {
basis <- svd(sRaw, nu = K, nv = 0)\$u            # generic R -> C -> LAPACK SVD
}

basis <- basis[, -K]   # mean adjustment implies a mean-passing subspace of dim K-1,
# conveniently arranged so that basis vector K is out-of-band
if(intercept) { basis <- cbind(R, basis) }

dimnames(basis) <- list(names(wx), 1L:ncol(basis))

# need to convert back to original time array -- the NAs in the original
if(any(is.na(x))) {
basis[which(!(wx %in% x[!is.na(x)])), ] <- rep(NA, ncol(basis))
}

a <- list(K = K, W = W, N = N, naive = naive)
attributes(basis) <- c(attributes(basis), a)
class(basis) <- c("slp", "basis", "matrix")
} else {
sMat <- sRaw
}
}
}
if(!returnS) {
basis
} else {
sMat
}
}
```

## Try the slp package in your browser

Any scripts or data that you put into this service are public.

slp documentation built on May 2, 2019, 2:39 a.m.