#' Temporal Functional Principal Component
#'
#' This function is used on step II of the STDF algorithm. It does not need to be used
#' directly, since it's mostly an wrapper to \code{fda} package function calls. It
#' behaves similarly to the \code{pca.fd} function, but rescales the
#' eigenfunctions/eigenvalues in a way that's needed by the \code{stdf} function.
#' Notice \code{tfpca} do not handle NA's in the data values. It might be useful do to some
#' kind of nearest neightbor imputation if possible. The \code{fda} function \code{data2fd} does
#' a least-squares fit of the basis function, which can be helpful in preprocessing the data.
#' The tuning parameter zeta must be passed manually.
#'
#' @param MatY The data organized in a matrix, with each time series in a column. Missing
#' values are currently not handled properly.
#' @param L The number of functional principal components to be returned.
#' @param t.fit Values of the time-domain at which a value is fitted.
#' @param zeta Smoothness parameter, must be provided. Defaults to 0 (no smoothing).
#' @param tbas Number of spline basis functions used in the smooth.
#'
#' @export
#' @return List of three items
#' \item{values}{numeric vector of L Eigenvalue estimates}
#' \item{vectors}{matrix of L columns with Eigenfunction estimates observed at t.fit}
#' \item{harmfd}{Internal use. Prediction of the Eigenfunction at new values of t,
#' say under vector t.pred, can be obtained with the next three
#' parameters. The syntax is
#' fda::eval.fd(t.pred, harmfd)*t(matrix(sqrt(nObs/etan), L, length(t.pred)))}
#' \item{nObs}{Internal use. }
#' \item{etan}{Internal use. }
#'
#' @examples
#' # This example uses fda's CanadianWeather dataset:
#' library(fda)
#' test <- tfpca(CanadianWeather$dailyAv[,,"Temperature.C"], 2, 1:365, 0, 20)
#' test2 <- tfpca(CanadianWeather$dailyAv[,,"Temperature.C"], 2, 1:365, 1e5, 20)
#' layout(matrix(1:4, byrow = TRUE, ncol = 2))
#' plot(test$vectors[,1], type="l", main="First Principal Component")
#' plot(test$vectors[,2], type="l", main="Second Principal Component")
#' plot(test2$vectors[,1], type="l", main="Smooth First Principal Component")
#' plot(test2$vectors[,2], type="l", main="Smooth Second Principal Component")
#'
#' @references
#' Ramsay, J. O. (2006) \emph{Functional Data Analysis}. New York: Springer.
#'
#' @seealso \code{\link{pca.fd}}, \code{\link{stdf}}
#' @keywords Spatial Statistics
#' @keywords Functional Data Analysis
tfpca <- function(MatY, L, t.fit, zeta = 0, tbas = 20){
nObs <- nrow(MatY)
if(nObs < tbas) stop("Number of basis > number of observations, please adjust \"tbas\" parameter\n")
rng <- range(t.fit)
timebasis <- fda::create.bspline.basis(rng, nbasis = tbas) # Smoothness controlled with # basis
timefdPar <- fda::fdPar(timebasis, lambda = zeta) # zeta = 0 => no additional smoothing
smtempfd <- fda::smooth.basis(t.fit, MatY, timebasis)$fd # Preliminary smoothing
pcafd <- fda::pca.fd(smtempfd, nharm = L, timefdPar)
harmfd <- pcafd[[1]]
fdmat <- fda::eval.fd(t.fit, harmfd)
etan <- apply(fdmat^2,2,sum) # Approximation of integral?
values <- pcafd$values[1:L]*etan/nObs # Truncates at L, multiply by |phi|, divide by n
vectors <- fdmat*t(matrix(sqrt(nObs/etan), L, nObs)) # Vectors at scale of data?
return(list(values = values,
vectors = vectors,
harmfd = harmfd,
nObs = nObs,
etan = etan))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.