Added full PACE-QUO functionality.

DESCRIPTION

@@ -11,7 +11,7 @@ Author: Xiongtao Dai,
     Jane-Ling Wang
 Maintainer: Pantelis Z. Hadjipantelis <[email protected]>
 Description: PACE is a versatile package that provides implementation of various methods of Functional Data Analysis (FDA) and Empirical Dynamics. The core of this package is Functional Principal Component Analysis (FPCA), a key technique for functional data analysis, for sparsely or densely sampled random trajectories and time courses, via the Principal Analysis by Conditional Estimation (PACE) algorithm. PACE is useful for the analysis of data that have been generated by a sample of underlying (but usually not fully observed) random trajectories. It does not rely on pre-smoothing of trajectories, which is problematic if functional data are sparsely sampled. PACE provides options for functional regression and correlation, for Longitudinal Data Analysis, the analysis of stochastic processes from samples of realized trajectories, and for the analysis of underlying dynamics.
-Depends: R (>= 3.1.1)
+Depends: R (>= 3.1.0)
 License: BSD_3_clause 
 LazyData: false
 Imports: Rcpp (>= 0.11.5), RcppEigen, rARPACK, gtools, Hmisc, caret, plot3D, MASS, pracma, numDeriv

R/FPCAder.R

@@ -3,6 +3,7 @@
 #' Derivative FPCA for dense or sparse functional data. 
 #' 
 #' @param fpcaObj A object of class FPCA returned by the function FPCA(). 
+#' @param variant A string specifying the methodology used ('simple', 'QUO'; default: 'simple')
 #'
 #' @examples
 #' set.seed(1)
@@ -17,31 +18,50 @@
 #' @export
 
 
-FPCAder <-  function (fpcaObj) {
+FPCAder <-  function (fpcaObj, variant = 'simple') {
   # Use FPCA object to get derivative information object 
 
   if (class(fpcaObj) != 'FPCA'){
     stop("FPCAder() requires an FPCA class object as basic input")
   }
-
-  fpcaObjDer <- list( 
-    phi = apply(fpcaObj$phi, 2, getDerivative, t= fpcaObj$workGrid),
-    mu = getDerivative(y = fpcaObj$mu, t = fpcaObj$obsGrid), 
-    obsGrid = fpcaObj$obsGrid)
-
-  class(fpcaObjDer) <- 'FPCAder'
-  return(fpcaObjDer) 
+  if( variant == 'simple') {
+    fpcaObjDer <- list( 
+      phi = apply(fpcaObj$phi, 2, getDerivative, t= fpcaObj$workGrid),
+      mu = getDerivative(y = fpcaObj$mu, t = fpcaObj$obsGrid), 
+      obsGrid = fpcaObj$obsGrid)
+
+    class(fpcaObjDer) <- 'FPCAder'
+    return(fpcaObjDer) 
+  } else if ( variant == 'QUO' ){
+    impSample <- fitted(fpcaObj);
+    impSampleDer <- t(apply( impSample,1,getDerivative, fpcaObj$workGrid));
+    N = dim(impSample)[1];
+    M = dim(impSample)[2];
+    L = makePACEinputs(IDs = rep(1:N,each=M), tVec=rep(fpcaObj$workGrid,N), as.vector(t(impSampleDer)))
+    prefpcaObjDer = FPCA(y= L$Ly, t= L$Lt)
+    fpcaObjDer = list( 
+      phi = prefpcaObjDer$phi, mu =  prefpcaObjDer$mu, obsGrid = fpcaObj$obsGrid)     
+
+    class(fpcaObjDer) <- 'FPCAder'
+   return(fpcaObjDer)
+  } else {
+    stop("Invalid FPCAder variant requested.")
+    return( NULL )
+  }
 } 
 
 getEnlargedGrid <- function(x){ 
   N <- length(x)
   return (  c( x[1] - 0.1 * diff(x[1:2]), x, x[N] + 0.1 * diff(x[(N-1):N])) )
 }
 
-getDerivative <- function(y,t){ 
+getDerivative <- function(y,t){
+  if( length(y) != length(t) ){
+    stop("getDerivative y/t lengths are unequal.")    
+  }
   newt = getEnlargedGrid(t)
-  newy = approxExtrap(x=t, y=y, xout= newt)$y 
-  return (numDeriv::grad( splinefun(newt, newy) , x = t ) )
+  newy = Hmisc::approxExtrap(x=t, y=y, xout= newt)$y 
+  return (numDeriv::grad( stats::splinefun(newt, newy) , x = t ) )
 }