R/fRegress.CV.R

In drtagkim/mcgillfdar: Functional Data Analysis

fRegress.CV <- function(y, xfdlist, betalist, wt=NULL, CVobs=1:N, ...)
{

# FREGRESS.CV computes cross-validated error sum of squares
# for scalar or functional responses. NOTE: ordinary and
# generalized cross validation scores are now returned by fRegress
# when scalar responses are used.

# last modified 29 October 2009 by Jim Ramsay

argList  <- fRegressArgCheck(y, xfdlist, betalist, wt)
yfdPar   <- argList$yfdPar
xfdlist  <- argList$xfdlist
betalist <- argList$betalist
wt       <- argList$wt 

p <- length(xfdlist)
N <- dim(xfdlist[[1]]$coef)[2]

M <- length(CVobs)

if (inherits(yfdPar, "numeric"))  {

    yvec   <- yfdPar
    SSE.CV <- 0
    errfd  <- c()
    for (m in 1:M) {
      i        <- CVobs[m]  
      xfdlisti <- vector("list",p)
      for (j in 1:p) {
        xfdj          <- xfdlist[[j]]
        if (inherits(xfdj, "numeric")) {
          betafdParj <- betalist[[j]]
          betafdj    <- betafdParj$fd
          basisj     <- betafdj$basis
          betarangej <- basisj$rangeval
          conbasisj  <- create.constant.basis(betarangej)
          xfdj       <- fd(matrix(xfdj,1,N), conbasisj)
        }
        basisj <- xfdj$basis
        coefj  <- xfdj$coefs
        if (dim(coefj)[1] == 1) coefj <- matrix(coefj[-i],1,N-1)
        else                    coefj <- as.matrix(coefj[,-i])
        xfdlisti[[j]] <- fd(coefj,basisj)
      }
      yveci         <- yvec[-i]
      wti           <- wt[-i]
      fRegressListi <- fRegress(yveci, xfdlisti, betalist, wti)
      betaestlisti  <- fRegressListi$betaestlist
      yhati <- 0
      for (j in 1:p) {
        betafdParj <- betaestlisti[[j]]
        betafdj    <- betafdParj$fd
        xfdj       <- xfdlist[[j]]
        bbasisj    <- betafdj$basis
        rangej     <- bbasisj$rangeval
        nfine      <- max(101, bbasisj$nbasis*10+1)
        tfine      <- seq(rangej[1], rangej[2], len=nfine)
        delta      <- tfine[2]-tfine[1]
        betavec    <- eval.fd(tfine, betafdj)
        xveci      <- eval.fd(tfine, xfdj[i])
        yhati      <- yhati + delta*(sum(xveci*betavec) -
                                    0.5*( xveci[1]    *betavec[1] +
                                          xveci[nfine]*betavec[nfine] ))
      }
      errfd[i] = yvec[i] - yhati;
      SSE.CV <- SSE.CV + errfd[i]^2
    }
 } else { 
    yfd      <- yfdPar$fd
    SSE.CV   <- 0
    errcoefs <- c()
    for(m in 1:N){
      i <-  CVobs[m]
      txfdlist <- xfdlist           
      for(k in 1:p){
        txfdlist[[k]] <- xfdlist[[k]][-i]
      }
      tres <- fRegress(yfd[-i],txfdlist,betalist,wt)
      yhat <- 0                       
      for(k in 1:p){
        yhat <- yhat + xfdlist[[k]][i]*tres$betaestlist[[k]]$fd
      }
      errfdi   <- yfd[i] - yhat
      SSE.CV   <- SSE.CV + inprod(errfdi,errfdi)
      errcoefs <- cbind(errcoefs,errfdi$coefs)
    }
    errfd <- fd(errcoefs,errfdi$basis)
    names(errfd$fdnames)[[3]] <- "Xval Errors"
}
return(list(SSE.CV=SSE.CV,errfd.cv=errfd))
}