fRegress.stderr <- function (y, y2cMap, SigmaE, returnMatrix = FALSE, ...) 
{
  
  #  FREGRESS.STDERR  computes standard error estimates for regression
  #       coefficient functions estimated by function FREGRESS.
  #
  #  Arguments:
  #
  #  Y            ... a list object produced by function FREGRESS with class name
  #                   'fRegress".  
  #                   This is indicated by Y in the arguments since R syntax
  #                   requires all of tghe fRegress family of functions to
  #                   use this notation.
  #  Y2CMAP       ... the matrix mapping from the vector of observed values
  #                   to the coefficients for the dependent variable.
  #                   This is output by function SMOOTH_BASIS.  If this is
  #                   supplied, confidence limits are computed, otherwise not.
  #  SIGMAE       ... Estimate of the covariances among the residuals.  This
  #                   can only be estimated after a preliminary analysis
  #                   with FREGRESS.
  #  RETURNMATRIX ... If False, a matrix in sparse storage model can be returned
  #               from a call to function BsplineS.  See this function for
  #               enabling this option.
  #
  #  Returns:
  #
  #  BETASTDERRLIST ... a list object, each list containing a fdPar object
  #                     for the standard error of a regression function.
  #  BVAR           ... the symmetric matrix of sampling variances and
  #                     covariances for the matrix of regression coefficients
  #                     for the regression functions.  These are stored
  #                     column-wise in defining BVARIANCE.
  #  C2BMAP         ... the matrix mapping from response variable coefficients
  #                     to coefficients for regression coefficients
  
  #  Last modified 16 December 2020
  
  #  retrieve objects from y
  
  yfdobj         <- y$yfdobj
  xfdlist        <- y$xfdlist
  betalist       <- y$betalist
  betaestlist    <- y$betaestlist
  yhatfdobj      <- y$yhatfdobj
  Cmat           <- y$Cmat
  Dmat           <- y$Dmat
  Cmatinv        <- y$Cmatinv
  wt             <- y$wt
  df             <- y$df
  betastderrlist <- y$betastderrlist
  YhatStderr     <- y$YhatStderr
  Bvar           <- y$Bvar
  c2bMap         <- y$c2bMap
  
  #  get number of independent variables
  
  p <- length(xfdlist)

  #  compute number of coefficients
  
  ncoef <- 0
  for (j in 1:p) {
    betaParfdj <- betalist[[j]]
    ncoefj <- betaParfdj$fd$basis$nbasis
    ncoef <- ncoef + ncoefj
  }
  
  if (inherits(yfdobj, "fdPar") || inherits(yfdobj, "fd")) {
    
    #  ----------------------------------------------------------------
    #           yfdobj is functional data object
    #  ----------------------------------------------------------------
    
    #  As of 2020, if yfdobj is an fdPar object, it is converted to an fd object.
    #  The added structure of the fdPar class is not used in any of the fRegress codes.
    #  The older versions of fda package used yfdPar as the name for the first member.
    
    if (inherits(yfdobj, "fdPar")) yfdobj <- yfdobj$fd
    
    #  get number of replications and basis information for YFDPAR
    
    N         <- dim(yfdobj$coefs)[2]
    ybasisobj <- yfdobj$basis
    rangeval  <- ybasisobj$rangeval
    ynbasis   <- ybasisobj$nbasis
    ninteg    <- max(501, 10 * ynbasis + 1)
    tinteg    <- seq(rangeval[1], rangeval[2], len = ninteg)
    deltat    <- tinteg[2] - tinteg[1]
    
    ybasismat <- eval.basis(tinteg, ybasisobj, 0, returnMatrix)
    
    #  compute BASISPRODMAT
    
    basisprodmat <- matrix(0, ncoef, ynbasis * N)
    
    mj2 <- 0
    for (j in 1:p) {
      betafdParj <- betalist[[j]]
      betabasisj <- betafdParj$fd$basis
      ncoefj <- betabasisj$nbasis
      bbasismatj <- eval.basis(tinteg, betabasisj, 0, returnMatrix)
      xfdj <- xfdlist[[j]]
      tempj <- eval.fd(tinteg, xfdj, 0, returnMatrix)
      mj1 <- mj2 + 1
      mj2 <- mj2 + ncoefj
      indexj <- mj1:mj2
      mk2 <- 0
      for (k in 1:ynbasis) {
        mk1 <- mk2 + 1
        mk2 <- mk2 + N
        indexk <- mk1:mk2
        tempk <- bbasismatj * ybasismat[, k]
        basisprodmat[indexj, indexk] <- deltat * crossprod(tempk, 
                                                           tempj)
      }
    }
    
    #  check dimensions of Y2CMAP
    
    y2cdim <- dim(y2cMap)
    if (y2cdim[1] != ynbasis || y2cdim[2] != dim(SigmaE)[1]) 
      stop("Dimensions of Y2CMAP not correct.")
    
    #  compute variances of regression coefficient function values
    
    Varc <- y2cMap %*% SigmaE %*% t(y2cMap)
    CVar <- kronecker(Varc, diag(rep(1, N)))
    C2BMap <- Cmatinv %*% basisprodmat
    Bvar <- C2BMap %*% CVar %*% t(C2BMap)
    nplot <- max(51, 10 * ynbasis + 1)
    tplot <- seq(rangeval[1], rangeval[2], len = nplot)
    betastderrlist <- vector("list", p)
    PsiMatlist <- vector("list", p)
    mj2 <- 0
    for (j in 1:p) {
      betafdParj <- betalist[[j]]
      betabasisj <- betafdParj$fd$basis
      ncoefj <- betabasisj$nbasis
      mj1 <- mj2 + 1
      mj2 <- mj2 + ncoefj
      indexj <- mj1:mj2
      bbasismat <- eval.basis(tplot, betabasisj, 0, returnMatrix)
      PsiMatlist <- bbasismat
      bvarj <- Bvar[indexj, indexj]
      bstderrj <- sqrt(diag(bbasismat %*% bvarj %*% t(bbasismat)))
      bstderrfdj <- smooth.basis(tplot, bstderrj, betabasisj)$fd
      betastderrlist[[j]] <- bstderrfdj
    }
    
    #  compute estimated variance-covariance matrix over plotting grid
    #  of fitted values
    
    YhatStderr <- matrix(0, nplot, N)
    B2YhatList <- vector("list", p)
    for (iplot in 1:nplot) {
      YhatVari <- matrix(0, N, N)
      tval <- tplot[iplot]
      for (j in 1:p) {
        Zmat <- eval.fd(tval, xfdlist[[j]])
        betabasisj <- betalist[[j]]$fd$basis
        PsiMatj <- eval.basis(tval, betabasisj)
        B2YhatMapij <- t(Zmat) %*% PsiMatj
        B2YhatList[[j]] <- B2YhatMapij
      }
      m2j <- 0
      for (j in 1:p) {
        m1j <- m2j + 1
        m2j <- m2j + betalist[[j]]$fd$basis$nbasis
        B2YhatMapij <- B2YhatList[[j]]
        m2k <- 0
        for (k in 1:p) {
          m1k <- m2k + 1
          m2k <- m2k + betalist[[k]]$fd$basis$nbasis
          B2YhatMapik <- B2YhatList[[k]]
          YhatVari <- YhatVari + 
            B2YhatMapij %*% Bvar[m1j:m2j,m1k:m2k] %*% t(B2YhatMapik)
        }
      }
      YhatStderr[iplot, ] <- matrix(sqrt(diag(YhatVari)), 1, N)
    }
  }
  
  else {
    
    #  ----------------------------------------------------------------
    #                   yfdobj is scalar or multivariate
    #  ----------------------------------------------------------------
    
    ymat <- as.matrix(yfdobj)
    N <- dim(ymat)[1]
    Zmat <- NULL
    for (j in 1:p) {
      xfdj <- xfdlist[[j]]
      if (inherits(xfdj, "fd")) {
        xcoef <- xfdj$coefs
        xbasis <- xfdj$basis
        betafdParj <- betalist[[j]]
        bbasis <- betafdParj$fd$basis
        Jpsithetaj <- inprod(xbasis, bbasis)
        Zmat <- cbind(Zmat, t(xcoef) %*% Jpsithetaj)
      }
      else if (inherits(xfdj, "numeric")) {
        Zmatj <- xfdj
        Zmat <- cbind(Zmat, Zmatj)
      }
    }
    
    #  compute linear mapping c2bMap takinging coefficients for
    #  response into coefficients for regression functions
    
    c2bMap <- Cmatinv %*% t(Zmat)
    y2bmap <- c2bMap
    Bvar <- y2bmap %*% as.matrix(SigmaE) %*% t(y2bmap)
    betastderrlist <- vector("list", p)
    mj2 <- 0
    for (j in 1:p) {
      betafdParj <- betalist[[j]]
      betabasisj <- betafdParj$fd$basis
      ncoefj <- betabasisj$nbasis
      mj1 <- mj2 + 1
      mj2 <- mj2 + ncoefj
      indexj <- mj1:mj2
      bvarj <- Bvar[indexj, indexj]
      xfdj <- xfdlist[[j]]
      if (inherits(xfdj, "fd")) {
        betarng <- betabasisj$rangeval
        ninteg <- max(c(501, 10 * ncoefj + 1))
        tinteg <- seq(betarng[1], betarng[2], len = ninteg)
        bbasismat <- eval.basis(tinteg, betabasisj, 0, 
                                returnMatrix)
        bstderrj <- sqrt(diag(bbasismat %*% bvarj %*% 
                                t(bbasismat)))
        bstderrfdj <- smooth.basis(tinteg, bstderrj, 
                                   betabasisj)$fd
      }
      else {
        bsterrj <- sqrt(diag(bvarj))
        onebasis <- create.constant.basis(betabasisj$rangeval)
        bstderrfdj <- fd(t(bstderrj), onebasis)
      }
      betastderrlist[[j]] <- bstderrfdj
    }
    
    #  compute estimated variance-covariance matrix for fitted values
    
    B2YhatList <- vector("list", p)
    YhatVari <- matrix(0, N, N)
    for (j in 1:p) {
      betabasisj <- betalist[[j]]$fd$basis
      Xfdj <- xfdlist[[j]]
      B2YhatMapij <- inprod(Xfdj, betabasisj)
      B2YhatList[[j]] <- B2YhatMapij
    }
    m2j <- 0
    for (j in 1:p) {
      m1j <- m2j + 1
      m2j <- m2j + betalist[[j]]$fd$basis$nbasis
      B2YhatMapij <- B2YhatList[[j]]
      m2k <- 0
      for (k in 1:p) {
        m1k <- m2k + 1
        m2k <- m2k + betalist[[k]]$fd$basis$nbasis
        B2YhatMapik <- B2YhatList[[k]]
        YhatVari <- YhatVari + B2YhatMapij %*% Bvar[m1j:m2j, 
                                                    m1k:m2k] %*% t(B2YhatMapik)
      }
    }
    YhatStderr <- matrix(sqrt(diag(YhatVari)), N, 1)
  }
  
  #  return output object of class fRegress
  
  fRegressList <- list(yfdobj = y$yfdobj, xfdlist = y$xfdlist, 
                       betalist = y$betalist, betaestlist = y$betaestlist, 
                       yhatfdobj = y$yhatfdobj, Cmat = y$Cmat, Dmat = y$Dmat, 
                       Cmatinv = y$Cmatinv, wt = y$wt, 
                       df = y$df, y2cMap = y2cMap, SigmaE = SigmaE, 
                       betastderrlist = betastderrlist, YhatStderr = YhatStderr, 
                       Bvar = Bvar, c2bMap = c2bMap)
  
  class(fRegressList) = "fRegress"
  return(fRegressList)
}