createDesignPlot improved; document lwls2dV2

DESCRIPTION

@@ -4,17 +4,42 @@ Title: PACE package for Functional Data Analysis and Empirical Dynamics
 URL: https://github.com/hadjipantelis/tPACE
 Version: 0.0.0.9000
 Date: 2015-08-14
-Author: Xiongtao Dai, 
-    Pantelis Z. Hadjipantelis, 
-    Hao Ji, 
-    Hans-Georg Mueller, 
+Author: Xiongtao Dai,
+    Pantelis Z. Hadjipantelis,
+    Hao Ji,
+    Hans-Georg Mueller,
     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.0)
-License: BSD_3_clause 
+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.0)
+License: BSD_3_clause
 LazyData: false
-Imports: Rcpp (>= 0.11.5), RcppEigen, gtools, Hmisc, plot3D, MASS, pracma, numDeriv
+Imports:
+    Rcpp (>= 0.11.5),
+    RcppEigen,
+    gtools,
+    Hmisc,
+    plot3D,
+    MASS,
+    pracma,
+    numDeriv
 LinkingTo: Rcpp, RcppEigen
-Suggests: testthat, rgl, aplpack
+Suggests:
+    testthat,
+    rgl,
+    aplpack
 NeedsCompilation: yes
+RoxygenNote: 5.0.1

NAMESPACE

@@ -1,22 +1,28 @@
-# Generated by roxygen2 (4.1.1): do not edit by hand
+# Generated by roxygen2: do not edit by hand
 
 S3method(fitted,FPCA)
 S3method(fitted,FregObj)
 S3method(print,FPCA)
 export(CrCovYX)
+export(CrCovYXfast)
 export(CrCovYZ)
 export(FPCA)
 export(FPCAder)
 export(FPCAregFunc)
+export(FPCAregFuncExp)
 export(FPCAregScalar)
+export(FPCAregScalarExp)
 export(FVPA)
 export(Rlwls1d)
+export(createBetaPlots)
 export(createCovPlot)
 export(createDesignPlot)
 export(createDiagnosticsPlot)
 export(createFuncBoxPlot)
+export(createModeOfVarPlot)
 export(createOutliersPlot)
 export(createScreePlot)
+export(lwls2dV2)
 export(makePACEinputs)
 export(sparsify)
 export(wiener)
@@ -32,5 +38,4 @@ importFrom(pracma,midpoint)
 importFrom(pracma,mod)
 importFrom(pracma,ones)
 importFrom(pracma,uniq)
-importFrom(rARPACK,eigs)
 useDynLib(tPACE)

R/RcppExports.R

@@ -21,10 +21,6 @@ interp2lin <- function(xin, yin, zin, xou, you) {
     .Call('tPACE_interp2lin', PACKAGE = 'tPACE', xin, yin, zin, xou, you)
 }
 
-RmullwlskCC <- function(bw, kernel_type, tPairs, cxxn, win, xgrid, ygrid, bwCheck) {
-    .Call('tPACE_RmullwlskCC', PACKAGE = 'tPACE', bw, kernel_type, tPairs, cxxn, win, xgrid, ygrid, bwCheck)
-}
-
 Rmullwlsk <- function(bw, kernel_type, tPairs, cxxn, win, xgrid, ygrid, bwCheck, transp = TRUE) {
     .Call('tPACE_Rmullwlsk', PACKAGE = 'tPACE', bw, kernel_type, tPairs, cxxn, win, xgrid, ygrid, bwCheck, transp)
 }
@@ -33,6 +29,10 @@ Rmullwlsk_old <- function(bw, kernel_type, tPairs, cxxn, win, xgrid, ygrid, bwCh
     .Call('tPACE_Rmullwlsk_old', PACKAGE = 'tPACE', bw, kernel_type, tPairs, cxxn, win, xgrid, ygrid, bwCheck, transp)
 }
 
+RmullwlskCC <- function(bw, kernel_type, tPairs, cxxn, win, xgrid, ygrid, bwCheck) {
+    .Call('tPACE_RmullwlskCC', PACKAGE = 'tPACE', bw, kernel_type, tPairs, cxxn, win, xgrid, ygrid, bwCheck)
+}
+
 Rrotatedmullwlsk <- function(bw, kernel_type, tPairs, cxxn, win, xygrid, npoly, bwCheck) {
     .Call('tPACE_Rrotatedmullwlsk', PACKAGE = 'tPACE', bw, kernel_type, tPairs, cxxn, win, xygrid, npoly, bwCheck)
 }

R/Rlwls1d.R

@@ -18,8 +18,11 @@
 
 Rlwls1d <- function( bw, kernel_type, win, xin, yin, xout, npoly = 1L, nder = 0L){
 
+  if (is.unsorted(xout))
+    stop('`xout` needs to be sorted in increasing order')
+
   # Deal with NA/NaN measurement values
-  NAinY = is.na(yin);
+  NAinY = is.na(xin) | is.na(yin) | is.na(win)
   if(any(NAinY)){
       win = win[!NAinY]
       xin = xin[!NAinY]

R/createDesignPlot.R

@@ -21,7 +21,7 @@
 #' createDesignPlot(sampWiener$tList, sort(unique(unlist(sampWiener$tList))))
 #' @export
 
-createDesignPlot = function(t, obsGrid = NULL, isColorPlot=FALSE, noDiagonal=TRUE, yname = NULL){
+createDesignPlot = function(t, obsGrid = NULL, isColorPlot=TRUE, noDiagonal=TRUE, yname = NULL, ...){
   if( is.null(obsGrid)){
     obsGrid = sort(unique(unlist(t)))
   }
@@ -41,16 +41,14 @@ createDesignPlot = function(t, obsGrid = NULL, isColorPlot=FALSE, noDiagonal=TRU
   #   xlab = 'Observed time grid', ylab = 'Observed time grid',
   #   main = titleString, col = 'white')
 
+  oldpty <- par()[['pty']]
+  par(pty="s")
   if(isColorPlot == TRUE){
-    par(pty="s")
-  	createColorPlot(res, obsGrid,titleString)
-  	par(pty="m")
+  	createColorPlot(res, obsGrid,titleString, ...)
   } else {
-    par(pty="s")
   	createBlackPlot(res, obsGrid,titleString)
-  	par(pty="m")
   }
-
+  par(pty=oldpty)
 }
 
 createBlackPlot = function(res, obsGrid,titleString ){
@@ -74,39 +72,47 @@ createBlackPlot = function(res, obsGrid,titleString ){
 
 }
 
-createColorPlot = function(res, obsGrid,titleString ){
-#  for(i in 1:length(obsGrid)){
-#    tmp = res[i,]
-#    idx = which(tmp > 0)
-#    if(length(idx) > 0){
-#      for(j in 1:length(idx)){
-#        points(obsGrid[i], obsGrid[idx[j]], col = searchCol(tmp[idx[j]]),
-#          pch = 19)
-#      }
-#    }
-#  }
-#  legend('right', c('1', '2', '3~5', '>=6'), pch = 14, 
-#    col = c('red', 'purple', 'green', 'blue'), title = 'Count')
-  res[res >4 ] = 4;
-  resVals = sort(unique(as.vector(res)));
-  colPalette = c('white', 'black', 'blue', 'green', 'red')
-  resColPalt = colPalette[resVals+1]
-
+createColorPlot = function(res, obsGrid,titleString, ... ){
+  res[res > 4] = 4;
+  # resVals = sort(unique(as.vector(res)));
+  # colPalette = c('white', 'black', 'blue', 'green', 'red')
+  # resColPalt = colPalette[resVals+1]
   # image(res, col= resColPalt, axes=FALSE, xlab = 'Observed support points', ylab = 'Observed support points', main = titleString)
 
-  u1 = as.vector(res)
-  u2 = as.vector(t(res))
-  t1 = rep(obsGrid, times = length(obsGrid) )
+  notZero <- res != 0
+  nnres <- res[notZero]
+  ddd <- list(...)
+
+  colVec <- c(`1`='black', `2`='blue', `3`='green', `4`='red')
+  if (!is.null(ddd[['col']]))
+    colVec[] <- ddd[['col']]
+
+  pchVec <- rep(19, length(colVec))
+  names(pchVec) <- names(colVec)
+  if (!is.null(ddd[['pch']]))
+    pchVec[] <- ddd[['pch']]
+
+  cexVec <- seq(from=0.3, by=0.1, length.out=length(colVec))
+  names(cexVec) <- names(colVec)
+  if (!is.null(ddd[['cex']]))
+    cexVec[] <- ddd[['cex']]
+
+  if (!is.null(ddd[['xlab']]))
+    xlab <- ddd[['xlab']]
+  else
+    xlab <- 'Observed time grid'
+
+  if (!is.null(ddd[['ylab']]))
+    ylab <- ddd[['ylab']]
+  else
+    ylab <- 'Observed time grid'
+
+  t1 = rep(obsGrid, times = length(obsGrid))
   t2 = rep(obsGrid, each = length(obsGrid)) 
-  plot(t1, t2, col= 'black' , t= 'n', xlab = 'Observed time grid', ylab = 'Observed time grid', main = titleString, pch = 19 )
+  plot(t1[notZero], t2[notZero], col= colVec[nnres], xlab=xlab, ylab=ylab, main = titleString, pch = pchVec[nnres], cex=cexVec[nnres] )
 
-  points(t1[u1 == 1], t2[u2 ==1], col= 'black', pch = 19, cex =0.3)
-  points(t1[u1 == 2], t2[u2 ==2], col= 'blue',  pch = 19, cex =0.4)
-  points(t1[u1 == 3], t2[u2 ==3], col= 'green', pch = 19, cex =0.5)
-  points(t1[u1 == 4], t2[u2 ==4], col= 'red',   pch = 19, cex =0.6)
-#  axis(1, obsGrid[round(seq(1,length(obsGrid), length.out=11))], obsGrid[round(seq(1,length(obsGrid), length.out=11))],col.axis="black")
-#  axis(2, obsGrid[round(seq(1,length(obsGrid), length.out=11))], obsGrid[round(seq(1,length(obsGrid), length.out=11))],col.axis="black")
-  legend('right', c('1', '2', '3', '4+'), pch = 19,  col = c('black','blue','green','red'), title = 'Count',bg='white' )
+  if (!identical(unique(nnres), 1))
+    legend('right', c('1', '2', '3', '4+'), pch = pchVec, col=colVec, pt.cex=cexVec, title = 'Count',bg='white' )
 }
 
 

R/designPlotCount.R

@@ -20,9 +20,8 @@ designPlotCount = function(t, obsGrid, noDiagonal, isColorPlot){
   N = length(obsGrid) # number of distinct observed time pts
   res = matrix(0, nrow = N, ncol = N)
 
-  for(i in 1:length(t)){
-    cur = t[[i]]
-    curidx = searchID(cur, obsGrid)
+  for(cur in t){
+    curidx = match(cur, obsGrid)
     if(isColorPlot == FALSE){
       res[curidx, curidx] = 1
     } else {
@@ -37,11 +36,11 @@ designPlotCount = function(t, obsGrid, noDiagonal, isColorPlot){
   return(res)
 }
 
-searchID = function(cur, obsGrid){
-  ni = length(cur)
-  id = rep(0, ni)
-  for(i in 1:ni){
-    id[i] = which(obsGrid == cur[i])
-  }
-  return(id)
-}
+# searchID = function(cur, obsGrid){
+  # ni = length(cur)
+  # id = rep(0, ni)
+  # for(i in 1:ni){
+    # id[i] = which(obsGrid == cur[i])
+  # }
+  # return(id)
+# }

R/lwls2dV2.R

@@ -1,16 +1,18 @@
-# Uses Pantelis' cpp code.
-# 2 dimensional local weighted least squares smoother. Only local linear smoother is implemented (no higher order, no derivative). 
-# bw: bandwidth, a scalar or a vector of length 2.
-# kern: kernel used: 'gauss', 'rect', 'gausvar', 'epan' (default), 'quar'
-# xin: an n by 2 dataframe or matrix of x-coordinate.
-# yin: a vector of y-coordinate.
-# win: a vector of weights on the observations. The number of count as in (maybe) raw covariance should be integrated into win.
-# xout1: a p1-vector of first output coordinate grid.
-# xout2: a p2-vector of second output coordinate grid. If both xout1 and xout2 are unspecified then the output gridpoints are set to the input gridpoints.
-# xout: alternative to xout1 and xout2. A matrix of p by 2 specifying the output points.
-
-# Returns a p1 by p2 matrix of fitted values if xout is not specified. Otherwise a vector of length p. 
+#' Two dimensional local linear kernel smoother.
+#'
+#' Two dimensional local weighted least squares smoother. Only local linear smoother for estimating the original curve is available (no higher order, no derivative). 
+#' @param bw A scalar or a vector of length 2 specifying the bandwidth.
+#' @param kern Kernel used: 'gauss', 'rect', 'gausvar', 'epan' (default), 'quar'.
+#' @param xin An n by 2 dataframe or matrix of x-coordinate.
+#' @param yin A vector of y-coordinate.
+#' @param win A vector of weights on the observations. 
+#' @param xout1: a p1-vector of first output coordinate grid. Defaults to the input gridpoints if left unspecified.
+#' @param xout2: a p2-vector of second output coordinate grid. Defaults to the input gridpoints if left unspecified.
+#' @param xout: alternative to xout1 and xout2. A matrix of p by 2 specifying the output points (may be inefficient if the size of \code{xout} is small).
+#' @return a p1 by p2 matrix of fitted values if xout is not specified. Otherwise a vector of length p corresponding to the rows of xout. 
+#' @export
 
+# Uses Pantelis' cpp code.
 lwls2dV2 <- function(bw, kern='epan', xin, yin, win=NULL, xout1=NULL, xout2=NULL, xout=NULL, subset=NULL, crosscov = FALSE, userNumCores = NULL ) {
   if (length(bw) == 1){
     bw <- c(bw, bw)
@@ -30,7 +32,7 @@ lwls2dV2 <- function(bw, kern='epan', xin, yin, win=NULL, xout1=NULL, xout2=NULL
   }
 
   if (!is.null(xout1) && !is.null(xout2) && !is.null(xout)) {
-    stop('Either xout1/xout2 or xout should be specified.')
+    stop('Either xout1/xout2 or xout should be specified, but not both.')
   }
 
   if (is.null(xout1))