project.basis <- function(y, argvals, basisobj, penalize=FALSE)
{
# Arguments for this function:
#
# 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.
# ARGVALS ... A vector of argument values. This must be of length
# length(Y) if Y is a vector or dim(Y)[1] otherwise.
# BASISOBJ ... A basis.fd object
# PENALIZE ... If TRUE, a penalty term is used to deal with a singular
# basis matrix. But this is not normally needed.
#
# Returns a coefficient vector or array. The first dimension is the number
# of basis functions and the other dimensions (if any) match
# the other dimensions of Y.
#
# Last modified: 29 October 2005
# Check BASISOBJ
if (!inherits(basisobj, "basisfd")) stop(
"BASISOBJ is not a basis object.")
#
# Calculate the basis and penalty matrices, using the default
# for the number of derivatives in the penalty.
#
basismat <- getbasismatrix(argvals, basisobj)
Bmat <- crossprod(basismat)
if (penalize) {
penmat <- eval.penalty(basisobj)
#
# Add a very small multiple of the identity to penmat
# and find a regularization parameter
#
penmat <- penmat + 1e-10 * max(penmat) * diag(dim(penmat)[1])
lambda <- (0.0001 * sum(diag(Bmat)))/sum(diag(penmat))
Cmat <- Bmat + lambda * penmat
} else {
Cmat <- Bmat
}
#
# Do the fitting by a simple solution of the
# equations taking into account smoothing
#
if (is.array(y) == FALSE) y <- as.array(y)
if(length(dim(y)) <= 2) {
Dmat <- crossprod(basismat, y)
coef <- symsolve(Cmat, Dmat)
} else {
nvar <- dim(y)[3]
coef <- array(0, c(basisobj$nbasis, dim(y)[2], nvar))
for(ivar in 1:nvar) {
Dmat <- crossprod(basismat, y[, , ivar])
coef[, , ivar] <- symsolve(Cmat, Dmat)
}
}
coef
}