https://github.com/cran/fda
Tip revision: ca5e2b4994971ec127b6a5ed2a08ce34abb2655c authored by J. O. Ramsay on 28 September 2021, 03:50:08 UTC
version 5.4.0
version 5.4.0
Tip revision: ca5e2b4
smooth.basis1.R
smooth.basis1 <- function (argvals=1:n, y, fdParobj,
wtvec=NULL, fdnames=NULL, covariates=NULL,
method="chol", dfscale=1, returnMatrix=FALSE)
{
# Arguments:
# ARGVALS A set of N argument values, set by default to equally spaced
# on the unit interval (0,1).
# Y an array containing values of curves
# If the array is a matrix, rows must correspond to argument
# values and columns to replications, and it will be assumed
# that there is only one variable per observation.
# If Y is a three-dimensional array, the first dimension
# corresponds to argument values, the second to replications,
# and the third to variables within replications.
# If Y is a vector, only one replicate and variable are assumed.
# FDPAROBJ A functional parameter or fdPar object. This object
# contains the specifications for the functional data
# object to be estimated by smoothing the data. See
# comment lines in function fdPar for details.
# This argument may also be either a FD object, or a
# BASIS object. In this case, the smoothing parameter
# LAMBDA is set to 0.
# WEIGHT A vector of N weights, set to one by default, that can
# be used to differentially weight observations, or
# a symmetric positive definite matrix of order N
# FDNAMES A cell of length 3 with names for
# 1. argument domain, such as "Time"
# 2. replications or cases
# 3. the function.
# COVARIATES A N by Q matrix Z of covariate values used to augment
# the smoothing function, where N is the number of
# data values to be smoothed and Q is the number of
# covariates. The process of augmenting a smoothing
# function in this way is often called "semi-parametric
# regression". The default is the null object NULL.
# METHOD The method for computing coefficients. The usual method
# computes cross-product matrices of the basis value matrix,
# adds the roughness penalty, and uses the Choleski decomposition
# of this to compute coefficients, analogous to using the
# normal equations in least squares fitting. But this approach,
# while fast, contributes unnecessary rounding error, and the qr
# decomposition of the augmented basis matrix is prefererable.
# But nothing comes for free, and the computational overhead of
# the qr approach can be a serious problem for large problems
# (n of 1000 or more). For this reason, the default is
# "method" = "chol", but if 'method' == 'qr', the qr
# decomposition is used.
# DFFACTOR A multiplier of df in GCV, set to one by default
#
# Returns a list containing:
# FDOBJ an object of class fd containing coefficients.
# DF a degrees of freedom measure.
# GCV a measure of lack of fit discounted for df.
# If the function is univariate, GCV is a vector
# containing the error sum of squares for each
# function, and if the function is multivariate,
# GCV is a NVAR by NCURVES matrix.
# COEF the coefficient matrix for the basis function
# expansion of the smoothing function
# SSE the error sums of squares.
# SSE is a vector or matrix of the same size as
# GCV.
# PENMAT the penalty matrix.
# Y2CMAP the matrix mapping the data to the coefficients.
#
# last modified 7 August 2020 by Jim Ramsay
# This version of smooth.basis, introduced in March 2011, permits ARGVALS
# to be a matrix, with the same dimensions as the first two dimensions of Y
# This allows the sampling points to vary from one record to another.
# This first section of code selects the version of smooth.basis to use
# depending on whether ARGVALS is a vector (case 1) or a matrix (case 2)
# The earlier version of smooth.basis is found at the end of the file where
# it is names smooth.basis1.
# ---------------------------------------------------------------------
# Check argments
# ---------------------------------------------------------------------
# check Y and set nrep, nvar and ndim
# set up matrix or array for coefficients of basis expansion,
# as well as names for replications and, if needed, variables
if (is.vector(y))y <- matrix(y,length(y),1)
dimy <- dim(y)
n <- dimy[1]
ycheck <- ycheck(y, n)
y <- ycheck$y
y0 <- y # preserve a copy of Y
nrep <- ycheck$ncurve
nvar <- ycheck$nvar
ndim <- ycheck$ndim
ydim <- dim(y)
# check ARGVALS
if (!is.numeric(argvals)) stop("'argvals' is not numeric.")
argvals <- as.vector(argvals)
# check fdParobj
fdParobj <- fdParcheck(fdParobj)
fdobj <- fdParobj$fd
lambda <- fdParobj$lambda
Lfdobj <- fdParobj$Lfd
penmat <- fdParobj$penmat
# check LAMBDA
if (lambda < 0) {
warning ("Value of 'lambda' was negative 0 used instead.")
lambda <- 0
}
# check WTVEC
wtlist <- wtcheck(n, wtvec)
wtvec <- wtlist$wtvec
onewt <- wtlist$onewt
matwt <- wtlist$matwt
# if (matwt) wtmat <- wtvec # else wtmat <- diag(as.vector(wtvec))
# set up names for first dimension of y
tnames <- dimnames(y)[[1]]
if (is.null(tnames)) tnames <- 1:n
# extract information from fdParobj
nderiv <- Lfdobj$nderiv
basisobj <- fdobj$basis
dropind <- basisobj$dropind
ndropind <- length(dropind)
nbasis <- basisobj$nbasis - ndropind
# get names for basis functions
names <- basisobj$names
if (ndropind > 0) {
names <- names[-dropind]
}
if (ndim == 2) {
coef <- matrix(0,nbasis,nrep)
ynames <- dimnames(y)[[2]]
vnames <- "value"
dimnames(coef) <- list(names, ynames)
}
if (ndim == 3) {
coef <- array(0,c(nbasis,nrep,nvar))
ynames <- dimnames(y)[[2]]
vnames <- dimnames(y)[[3]]
dimnames(coef) <- list(names, ynames, vnames)
}
# check COVARIATES and set value for q, the number of covariates
if (!is.null(covariates)) {
if (!is.numeric(covariates)) {
stop(paste("smooth_basis_LS:covariates",
"Optional argument COVARIATES is not numeric."))
}
if (dim(covariates)[1] != n) {
stop(paste("smooth_basis_LS:covariates",
"Optional argument COVARIATES has incorrect number of rows."))
}
q <- dim(covariates)[2]
} else {
q <- 0
beta. <- NULL
}
# set up names for first dimension of y
tnames <- dimnames(y)[[1]]
if (is.null(tnames)) tnames <- 1:n
# ----------------------------------------------------------------
# set up the linear equations for smoothing
# ----------------------------------------------------------------
# set up matrix of basis function values
basismat <- eval.basis(as.vector(argvals), basisobj, 0, returnMatrix)
if (method == "chol") {
# -----------------------------------------------------------------
# use the default choleski decomposition of the crossproduct of the
# basis value matrix plus the roughness penalty
# -----------------------------------------------------------------
if (n > nbasis + q || lambda > 0) {
# augment BASISMAT0 and BASISMAT by the covariate matrix
# if it is supplied
if (!is.null(covariates)) {
ind1 <- 1:n
ind2 <- (nbasis+1):(nbasis+q)
basismat <- as.matrix(basismat)
basismat <- cbind(basismat, matrix(0,dim(basismat) [1],q))
basismat[ind1,ind2] <- covariates
}
# Compute the product of the basis and weight matrix
if (matwt) {
wtfac <- chol(wtvec)
basisw <- wtvec %*% basismat
} else {
rtwtvec <- sqrt(wtvec)
rtwtmat <- matrix(rtwtvec,n,nrep)
basisw <- (wtvec %*% matrix(1,1,nbasis+q))*basismat
}
# the weighted crossproduct of the basis matrix
Bmat <- t(basisw) %*% basismat
Bmat0 <- Bmat
# set up right side of normal equations
if (ndim < 3) {
Dmat <- t(basisw) %*% y
} else {
Dmat <- array(0, c(nbasis+q, nrep, nvar))
for (ivar in 1:nvar) {
Dmat[,,ivar] <- crossprod(basisw,y[,,ivar])
}
}
if (lambda > 0) {
# smoothing required, add the contribution of the penalty term
if (is.null(penmat)) penmat <- eval.penalty(basisobj, Lfdobj)
Bnorm <- sqrt(sum(diag(t(Bmat0) %*% Bmat0)))
pennorm <- sqrt(sum(penmat*penmat))
condno <- pennorm/Bnorm
if (lambda*condno > 1e12) {
lambda <- 1e12/condno
warning(paste("lambda reduced to",lambda,
"to prevent overflow"))
}
if (!is.null(covariates)) {
penmat <- rbind(cbind(penmat, matrix(0,nbasis,q)),
cbind(matrix(0,q,nbasis), matrix(0,q,q)))
}
Bmat <- Bmat0 + lambda*penmat
} else {
penmat <- NULL
Bmat <- Bmat0
}
# compute inverse of Bmat
Bmat <- (Bmat+t(Bmat))/2
Lmat <- try(chol(Bmat), silent=TRUE)
if (inherits(Lmat, "try-error")) {
Beig <- eigen(Bmat, symmetric=TRUE)
BgoodEig <- (Beig$values>0)
Brank <- sum(BgoodEig)
if (Brank<dim(Bmat)[1])
warning("Matrix of basis function values has rank ",
Brank, " < dim(fdobj$basis)[2] = ",
length(BgoodEig), " ignoring null space")
goodVec <- Beig$vectors[, BgoodEig]
Bmatinv <- (goodVec %*% (Beig$values[BgoodEig] * t(goodVec)))
} else {
Lmatinv <- solve(Lmat)
Bmatinv <- Lmatinv %*% t(Lmatinv)
}
# compute coefficient matrix by solving normal equations
if (ndim < 3) {
coef <- Bmatinv %*% Dmat
if (!is.null(covariates)) {
beta. <- as.matrix(coef[(nbasis+1):(nbasis+q),])
coef <- as.matrix(coef[1:nbasis,])
} else {
beta. <- NULL
}
} else {
coef <- array(0, c(nbasis, nrep, nvar))
if (!is.null(covariates)) {
beta. <- array(0, c(q, nrep, nvar))
} else {
beta. <- NULL
}
for (ivar in 1:nvar) {
coefi <- Bmatinv %*% Dmat[,,ivar]
if (!is.null(covariates)) {
beta.[,,ivar] <- coefi[(nbasis+1):(nbasis+q),]
coef[,,ivar] <- coefi[1:nbasis,]
} else {
coef[,,ivar] <- coefi
}
}
}
} else {
if (n == nbasis + q) {
# code for n == nbasis, q == 0, and lambda == 0
if (ndim==2) {
coef <- solve(basismat, y)
} else {
for (ivar in 1:var)
coef[1:n, , ivar] <- solve(basismat, y[,,ivar])
}
penmat <- NULL
} else {
# n < nbasis+q: this is treated as an error
stop("The number of basis functions = ", nbasis+q, " exceeds ",
n, " = the number of points to be smoothed.")
}
}
} else {
# -------------------------------------------------------------
# computation of coefficients using the qr decomposition of the
# augmented basis value matrix
# -------------------------------------------------------------
if (n > nbasis || lambda > 0) {
# Multiply the basis matrix and the data pointwise by the square root
# of the weight vector if the weight vector is not all ones.
# If the weights are in a matrix, multiply the basis matrix by its
# Choleski factor.
if (!onewt) {
if (matwt) {
wtfac <- chol(wtvec)
basismat.aug <- wtfac %*% basismat
if (ndim < 3) {
y <- wtfac %*% y
} else {
for (ivar in 1:nvar) {
y[,,ivar] <- wtfac %*% y[,,ivar]
}
}
} else {
rtwtvec <- sqrt(wtvec)
basismat.aug <- matrix(rtwtvec,n,nbasis) * basismat
if (ndim < 3) {
y <- matrix(rtwtvec,n,nrep) * y
} else {
for (ivar in 1:nvar) {
y[,,ivar] <- matrix(rtwtvec,n,nrep) * y[,,ivar]
}
}
}
} else {
basismat.aug <- basismat
}
# set up additional rows of the least squares problem for the
# penalty term.
if (lambda > 0) {
if (is.null(penmat)) penmat <- eval.penalty(basisobj, Lfdobj)
eiglist <- eigen(penmat)
Dvec <- eiglist$values
Vmat <- eiglist$vectors
# Check that the lowest eigenvalue in the series that is to be
# kept is positive.
neiglow <- nbasis - nderiv
naug <- n + neiglow
if (Dvec[neiglow] <= 0) {
stop(paste("smooth_basis:eig",
"Eigenvalue(NBASIS-NDERIV) of penalty matrix ",
"is not positive check penalty matrix."))
}
# Compute the square root of the penalty matrix in the subspace
# spanned by the first N - NDERIV eigenvectors
indeig <- 1:neiglow
penfac <- Vmat[,indeig] %*% diag(sqrt(Dvec[indeig]))
# Augment basismat by sqrt(lambda)*t(penfac)
basismat.aug <- rbind(basismat.aug, sqrt(lambda)*t(penfac))
# Augment data vector by n - nderiv 0's
if (ndim < 3) {
y <- rbind(y, matrix(0,nbasis-nderiv,nrep))
} else {
y <- array(0,c(ydim[1]+nbasis-nderiv,ydim[2],ydim[3]))
y[1:ydim[1],,] <- y0
ind1 <- (1:(nbasis-nderiv)) + ydim[1]
for (ivar in 1:nvar) {
y[ind1,,ivar] <- matrix(0,nbasis-nderiv,nrep)
}
}
} else {
penmat <- NULL
}
# augment BASISMAT0 and BASISMAT by the covariate matrix
# if it is supplied
if (!is.null(covariates)) {
ind1 <- 1:n
ind2 <- (nbasis+1):(nbasis+q)
basismat.aug <- cbind(basismat.aug, matrix(0,naug,q))
if (!onewt) {
if (matwt) {
basismat.aug[ind1,ind2] <- wtfac %*% covariates
} else {
wtfac <- matrix(rtwtvec,n,q)
basismat.aug[ind1,ind2] <- wtfac*covariates
}
} else {
basismat.aug[ind1,ind2] <- covariates
}
penmat <- rbind(cbind(penmat, matrix(0,nbasis,q)),
cbind(matrix(0,q,nbasis), matrix(0,q,q)))
}
# solve the least squares problem using the QR decomposition with
# one iteration to improve accuracy
qr <- qr(basismat.aug)
if (ndim < 3) {
coef <- qr.coef(qr,y)
if (!is.null(covariates)) {
beta. <- coef[ind2,]
coef <- coef[1:nbasis,]
} else {
beta. <- NULL
}
} else {
coef <- array(0, c(nbasis, nrep, nvar))
if (!is.null(covariates)) {
beta. <- array(0, c(q,nrep,nvar))
} else {
beta. <- NULL
}
for (ivar in 1:nvar) {
coefi <- qr.coef(qr,y[,,ivar])
if (!is.null(covariates)) {
beta.[,,ivar] <- coefi[ind2,]
coef[,,ivar] <- coefi[1:nbasis,]
} else {
coef[,,ivar] <- coefi
}
}
}
} else {
if (n == nbasis + q) {
# code for n == nbasis and lambda == 0
if (ndim==2){
coef <- solve(basismat, y)
} else {
for (ivar in 1:var)
coef[,,ivar] <- solve(basismat, y[,,ivar])
}
penmat <- NULL
} else {
stop(paste("The number of basis functions = ", nbasis, " exceeds ",
n, " = the number of points to be smoothed. "))
}
}
}
# ----------------------------------------------------------------
# compute SSE, yhat, GCV and other fit summaries
# ----------------------------------------------------------------
# compute map from y to c
if (onewt) {
temp <- t(basismat) %*% basismat
if (lambda > 0) {
temp <- temp + lambda*penmat
}
L <- chol(temp)
MapFac <- solve(t(L),t(basismat))
y2cMap <- solve(L,MapFac)
} else {
if(matwt){
temp <- t(basismat) %*% wtvec %*% basismat
} else {
temp <- t(basismat) %*% (as.vector(wtvec)*basismat)
}
if (lambda > 0) {
temp <- temp + lambda*penmat
}
L <- chol((temp+t(temp))/2)
MapFac <- solve(t(L),t(basismat))
if(matwt){
y2cMap <- solve(L, MapFac%*%wtvec)
} else {
y2cMap <- solve(L,MapFac*rep(as.vector(wtvec), e=nrow(MapFac)))
}
}
# compute degrees of freedom of smooth
df. <- sum(diag(y2cMap %*% basismat))
# compute error sum of squares
if (ndim < 3) {
yhat <- basismat[,1:nbasis] %*% coef
SSE <- sum((y[1:n,] - yhat)^2)
if (is.null(ynames)) ynames <- dimnames(yhat)[[2]]
} else {
SSE <- 0
yhat <- array(0,c(n, nrep, nvar))
dimnames(yhat) <- list(dimnames(basismat)[[1]],
dimnames(coef)[[2]],
dimnames(coef)[[3]])
for (ivar in 1:nvar) {
yhat[,,ivar] <- basismat[,1:nbasis] %*% coef[,,ivar]
SSE <- SSE + sum((y[1:n,,ivar] - yhat[,,ivar])^2)
}
if (is.null(ynames))ynames <- dimnames(yhat)[[2]]
if (is.null(vnames))vnames <- dimnames(yhat)[[2]]
}
if (is.null(ynames)) ynames <- paste("rep", 1:nrep, sep="")
if (is.null(vnames)) vnames <- paste("value", 1:nvar, sep="")
# compute GCV index
if (!is.null(df.) && df. < n) {
if (ndim < 3) {
gcv <- rep(0,nrep)
for (i in 1:nrep) {
SSEi <- sum((y[1:n,i] - yhat[,i])^2)
gcv[i] <- (SSEi/n)/((n - df.)/n)^2
}
if (ndim > 1) names(gcv) <- ynames
} else {
gcv <- matrix(0,nrep,nvar)
for (ivar in 1:nvar) {
for (i in 1:nrep) {
SSEi <- sum((y[1:n,i,ivar] - yhat[,i,ivar])^2)
gcv[i,ivar] <- (SSEi/n)/((n - df.)/n)^2
}
}
dimnames(gcv) <- list(ynames, vnames)
}
} else {
gcv <- NULL
}
#------------------------------------------------------------------
# Set up the functional data objects for the smooths
# ------------------------------------------------------------------
# set up default fdnames
if (is.null(fdnames)) {
fdnames <- list(time=tnames, reps=ynames, values=vnames)
}
# set up the functional data object
if (ndim < 3) {
coef <- as.matrix(coef)
fdobj <- fd(coef[1:nbasis,], basisobj, fdnames)
} else {
fdobj <- fd(coef[1:nbasis,,], basisobj, fdnames)
}
# return penalty matrix to original state if there were covariates
if (!is.null(penmat) && !is.null(covariates))
penmat <- penmat[1:nbasis,1:nbasis]
# assemble the fdSmooth object returned by the function
smoothlist <- list(fd=fdobj, df=df., gcv=gcv, beta=beta.,
SSE=SSE, penmat=penmat, y2cMap=y2cMap,
argvals=argvals, y=y0)
class(smoothlist) <- "fdSmooth"
return(smoothlist)
}