#' Predict FPC scores and curves for a new sample given an FPCA object
#'
#' Return a list containing the matrix with the first k FPC scores according to conditional expectation or numerical integration, the matrix of predicted trajectories and the prediction work grid.
#'
#' @param object An FPCA object.
#' @param newLy A list of \emph{n} vectors containing the observed values for each individual.
#' @param newLt A list of \emph{n} vectors containing the observation time points for each individual corresponding to y.
#' @param sigma2 The user-defined measurement error variance. A positive scalar. (default: rho if applicable, otherwise sigma2 if applicable, otherwise 0 if no error. )
#' @param K The scalar defining the number of clusters to define; (default: from user-specified FPCA Object).
#' @param xiMethod The integration method used to calculate the functional principal component scores
#' (standard numerical integration 'IN' or conditional expectation 'CE'); default: 'CE'.
#' @param ... Not used.
#'
#' @return A list containing the following fields:
#' \item{scores}{A matrix of size \emph{n}-by-\emph{k} which comprise of the predicted functional principal component scores.}
#' \item{predCurves}{A matrix of size \emph{n}-by-\emph{l} where \emph{l} is the length of the work grid in \emph{object}.}
#' \item{predGrid}{A vector of length \emph{l} which is the output grid of the predicted curves. This is same is the workgrid of \emph{object}.}
#'
#' @examples
#' \donttest{
#' set.seed(1)
#' n <- 50
#' pts <- seq(0, 1, by=0.05)
#' # The first n samples are for training and the rest testing
#' sampWiener <- Wiener(2 * n, pts)
#' sparsity <- 2:5
#' train <- Sparsify(sampWiener[seq_len(n), , drop=FALSE], pts, sparsity)
#' test <- Sparsify(sampWiener[seq(n + 1, 2 * n), , drop=FALSE], pts, sparsity)
#' res <- FPCA(train$Ly, train$Lt)
#' pred <- predict(res, test$Ly, test$Lt, K=3)
#' plot(pred$predGrid, pred$predCurves[1,])
#' }
#'
#' @method predict FPCA
#' @export
predict.FPCA <- function(object, newLy, newLt, sigma2 = NULL, K = NULL, xiMethod = 'CE', ...){
fpcaObj = object;
if(toString(class(fpcaObj)) != "FPCA"){
stop('Please provide a valid FPCA object.')
}
if(! all.equal( sapply(newLt, length), sapply(newLy, length) ) ){
stop('The size of the vectors in newLt and newLy differ. They must be equal.')
}
# Standard data checks/massasing as in FPCA()
CheckData(newLy, newLt)
inputData <- HandleNumericsAndNAN(Lt = newLt, Ly = newLy);
newLy <- inputData$Ly;
newLt <- inputData$Lt;
if(is.null(sigma2)){
if(fpcaObj$optns$methodRho=="vanilla"){
sigma2=fpcaObj$sigma2
}else{
sigma2 <- ifelse( !is.null(fpcaObj$rho), fpcaObj$rho,
ifelse( !is.null(fpcaObj$sigma2), fpcaObj$sigma2,
ifelse( !fpcaObj$optns$error, 0, stop('sigma2 cannot be determined.'))))
}
} else {
if(!is.numeric(sigma2) ){
stop('sigma2 is not numeric; are you sure?')
}
}
if(is.null(K)){
K = fpcaObj$selectK
}
if(K > fpcaObj$selectK){
stop( paste0( collapse = '', 'You cannot get FPC scores for more components than what is already available. (', fpcaObj$selectK ,').' ))
}
if( !(xiMethod %in% c('CE','IN')) ){
stop( paste0( collapse = '', 'Unrecognised method to calculate the FPC scores.'))
}
if( (fpcaObj$optns$dataType == 'Sparse') && (xiMethod == 'IN')){
stop( 'Trapezoid Numerical intergration (IN) is invalid for sparse data.')
}
rangeOfOrigData <- range(unlist(fpcaObj$workGrid))
rangeOfNewData <- range(unlist(newLt))
if(rangeOfNewData[1] < rangeOfOrigData[1]){
stop("The new data's lower range is below the original data's lower range.")
} else{
if(rangeOfNewData[2] > rangeOfOrigData[2]){
stop("The new data's upper range is above the original data's upper range.")
}
}
MuObs = ConvertSupport(fromGrid = fpcaObj$workGrid, toGrid = fpcaObj$obsGrid, mu = fpcaObj$mu)
CovObs = ConvertSupport(fromGrid = fpcaObj$workGrid, toGrid = fpcaObj$obsGrid, Cov = fpcaObj$fittedCov)
PhiObs = ConvertSupport(fromGrid = fpcaObj$workGrid, toGrid = fpcaObj$obsGrid, phi = fpcaObj$phi)
# Get scores
if ( xiMethod == 'CE') {
scoresObj <- GetCEScores(y = newLy, t = newLt, optns = fpcaObj$optns, mu = MuObs, obsGrid = fpcaObj$obsGrid, sigma2 = sigma2,
fittedCov = CovObs, lambda = fpcaObj$lambda, phi = PhiObs)
finalXiEst <- t(do.call(cbind, scoresObj['xiEst', ]))[, seq_len(K), drop=FALSE]
} else if (xiMethod == 'IN') {
scoresObj <- mapply(function(yvec,tvec)
GetINScores(yvec, tvec,optns = fpcaObj$optns,obsGrid = fpcaObj$obsGrid,mu = MuObs,lambda =fpcaObj$lambda ,phi = PhiObs,sigma2 = sigma2),newLy,newLt)
finalXiEst <- t(do.call(cbind, scoresObj['xiEst', ]))[, seq_len(K), drop=FALSE]
}
#return(finalXiEst)
#Get predicted trajectories
trajectoryEst <- t(fpcaObj$mu+fpcaObj$phi[,1:K]%*%t(finalXiEst))
return(list(scores = finalXiEst, predCurves = trajectoryEst, predGrid = fpcaObj$workGrid))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.