fRegress.stderr <- function(fRegressList, y2cMap, SigmaE) {

#  FREGRESS.STDERR  computes standard error estimates for regression
#       coefficient functions estimated by function FREGRESS.
#
#  Arguments:
#
#  FREGRESSLIST ... a list object produced by function FREGRESS.
#  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.
#
#  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 2007.11.28 by Spencer Graves
#  previosly modified 5 March 2005

#  get number of independent variables

xfdlist  <- fRegressList$xfdlist
yfdPar   <- fRegressList$yfdPar
betalist <- fRegressList$betalist
Cmatinv  <- fRegressList$Cmatinv

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(yfdPar, "fdPar") || inherits(yfdPar, "fd")) {
	
    #  ----------------------------------------------------------------
    #           YFDPAR is functional for a functional parameter
    #  ----------------------------------------------------------------

    if (inherits(yfdPar, "fd")) yfdPar <- fdPar(yfdPar)

    #  get number of replications and basis information for YFDPAR

    yfd       <- yfdPar$fd
    N         <- dim(yfd$coefs)[2]
    ybasisobj <- yfdPar$fd$basis
    rangeval  <- ybasisobj$rangeval
    ynbasis   <- ybasisobj$nbasis
    nfine     <- max(501,10*ynbasis+1)
    tfine     <- seq(rangeval[1], rangeval[2], len=nfine)
    deltat    <- tfine[2] - tfine[1]

    ybasismat <- eval.basis(tfine, ybasisobj)

    #  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(tfine, betabasisj)
        xfdj       <- xfdlist[[j]]
        tempj      <- eval.fd(tfine, xfdj)
        #  row indices of BASISPRODMAT to fill
        mj1    <- mj2 + 1
        mj2    <- mj2 + ncoefj
        indexj <- mj1:mj2
        #  inner products of beta basis and response basis
        #    weighted by covariate basis functions
        mk2 <- 0
        for (k in 1:ynbasis) {
            #  row indices of BASISPRODMAT to fill
            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
		
    c2bMap    <- Cmatinv %*% basisprodmat
    VarCoef   <- y2cMap %*% SigmaE %*% t(y2cMap)
    CVariance <- kronecker(VarCoef,diag(rep(1,N)))
    bvar      <- c2bMap %*% CVariance %*% t(c2bMap)
    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
        bbasismat  <- eval.basis(tfine, betabasisj)
        bvarj      <- bvar[indexj,indexj]
        bstderrj   <- sqrt(diag(bbasismat %*% bvarj %*% t(bbasismat)))
        bstderrfdj <- data2fd(bstderrj, tfine, betabasisj)
        betastderrlist[[j]] <- bstderrfdj
    }

} else {
	
    #  ----------------------------------------------------------------
    #                   YFDPAR is scalar or multivariate
    #  ----------------------------------------------------------------

    ymat <- as.matrix(yfdPar)
    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, "double")) {
            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
            nfine      <- max(c(501,10*ncoefj+1))
            tfine      <- seq(betarng[1], betarng[2], len=nfine)
            bbasismat  <- eval.basis(tfine, betabasisj)
            bstderrj   <- sqrt(diag(bbasismat %*% bvarj %*% t(bbasismat)))
            bstderrfdj <- data2fd(bstderrj, tfine, betabasisj)
        } else {
	         bsterrj    <- sqrt(diag(bvarj))
#	         bstderrfdj <- fd(t(bstderrj), onebasis)
	         bstderrfdj <- fd(t(bstderrj), ynbasis)
        }
        betastderrlist[[j]] <- bstderrfdj
    }
}

stderrList <-
       list(betastderrlist = betastderrlist,
            bvar           = bvar,
            c2bMap         = c2bMap)

return(stderrList)

}