R/pda.fd.R In fda: Functional Data Analysis

Documented in pda.fd

```pda.fd  <-  function(xfdlist, bwtlist=NULL, awtlist=NULL, ufdlist=NULL,
nfine=501, returnMatrix=FALSE)
{
#  PDA_FD computes the basis function expansions of the
#  estimates of the coefficient functions a_k(t) and b_j(t)
#  in the possibly nonhomogeneous linear differential operator
#
#    Lx(t) =
#       b_0(t)x(t) + b_1(t)Dx(t) + ... + b_{M-1}D^{M-1}x(t) + D^M x(t)
#       - a_1(t)u_1(t) - ... - a_k(t)u_K(t)
#
#  of order M = DIFEORDER that minimizes in a least squares sense the residual
#  functions f(t) = Lx(t).
#
#  The J equations may be of different orders and have different numbers
#     of forcing functions.  In the rather complicated description of the
#     arguments below, we use M_j to stand for the order of the jth equation.
#
#  If (DIFEORDER = 0, PDALIST fits the varying coefficient or pointwise
#  linear model using the functions x(t) as dependent variables and
#  the forcing functions u(t) as indep}ent variables.  In this case,
#  there must be at least one forcing function.
#
#  The functions x(t) are in functional data object XFDOBJ.
#  The forcing functions u_k(t) are in functional data object UFDOBJ.
#  The coefficient functions for u_k(t) and x(t) are expanded in terms of the
#  basis functions specified in AWTLIST and BWTLIST, respectively.
#
#
#  Arguments:
#  XFDLIST     list vector of length J of functional data objects for the
#                 J functions whose derivatives define the DIFE.
#  UFDLIST     list of length J whose components are themselves lists.
#                  The jth component of this list contains a list vector of
#                  length K_j of functional data objects for the
#                 independent variables or u-variables forcing the equation
#                  for the jth variable.
#              In the special univariate case where there is only a single
#                  DIFE, UFDLIST can have components which are the
#                  functional data objects for the forcing functions, rather
#                  than being a list with a single component list containing
#                  these functions.
#  BWTLIST     list vector of length J, the jth component of which is a list
#                   vector of length J, each component which is itself a list
#                   vector of length M_k equal to the order of the kth
#                   differential equation in the containing list.
#                   Each component of these lists within lists within a list
#                   is a functional parameter object defining a weighting
#                   coefficient function.
#              that is, BWTLIST is three levels or layers of lists, the top
#                   level a single list of length J, the second level a
#                   set of J lists, each corresponding to an equation, and the
#                   third level containing lists of length M_k containing
#                   coefficient functions defining the contribution of the
#                   m_kth derivative of variable k to the equation for
#                   variable j.
#              that is, if the orders of all the equations were the same
#                   and R supported list arrays, their dimensions would be
#                   J, J and M.
#              However, in the special case of J = 1, where only M
#                   coefficients are required, BWTLIST may be a simple
#                   list vector of length M.
#  AWTLIST     list of length J whose components are themselves lists.
#                 The jth component of this list contains a list vector of
#                  length K_j of functional parameter objects for the
#                  coefficient functions a_{jk}(t) multipling the
#                  corresponding forcing function in UFDLIST.
#              In the special univariate case where there is only a single
#                  DIFE, AWTLIST can have components which are the
#                  functional parameter objects for the forcing functions,
#                  rather than being a list with a single component list
#                  containing these functional parameter objects.
#  NFINE       number of sampling points for numerical integration, set by
#                  default to 501, but adjusted as required to define a mesh
#                  that is sufficient fine as to give a satisfactory
#                  approximation to an integral.

#  The value in each component of lists (within lists) XFDLIST and UFDLIST is a
#      scalar FD object.
#  The value in each component of lists (within lists) AWTLIST AND BWTLIST is a
#      scalar FDPAR object.

#  Returns:
#  BWTLIST     list structure identical to that of the argument BWTLIST
#  RESFDLIST   FD object for residual functions.
#  AWTLIST     list structure identical to that of the argument AWTLIST
#  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.

#  check dimensions of the lists

# check XFDLIST

if (inherits(xfdlist, "fd")) xfdlist = list(xfdlist)

if (!inherits(xfdlist, "list")) stop(
"XFDLIST is neither a list or a FD object")

nvar <- length(xfdlist)

#  ----------------------------------------------------------------
#     For efficiency, there are two versions of this code:
#     one for a single variable, and another for multiple variables.
#  ----------------------------------------------------------------

if (nvar == 1) {

#  ----------------------------------------------------------------
#                   Single variable case
#  ----------------------------------------------------------------

difeorder <- length(bwtlist)
difeordp1 <- difeorder + 1

xfdobj <- xfdlist[[1]]
xbasis <- xfdobj\$basis
xcoef  <- xfdobj\$coefs
xrange <- xbasis\$rangeval

#  check the dimensions of UFDLIST and AWTLIST and get number of forcing
#    functions NFORCE

if (is.null(ufdlist) | is.null(awtlist)) {
nforce  <- 0
} else {
if (inherits(ufdlist[[1]], "list")) {
#  UFDLIST is a list with a single component that is a list.
#  convert to a list of length NFORCE.
nforce <- length(ufdlist[[1]])
temp <- vector("list", nforce)
for (iu in 1:nforce) temp[[iu]] <- ufdlist[[1]][[iu]]
ufdlist <- temp
} else {
nforce <- length(ufdlist)
}
if (inherits(awtlist[[1]], "list")) {
#  AWTLIST is a list with a single component that is a list.
#  convert to a list of length NFORCE.
if (length(awtlist[[1]]) != nforce)
stop("The length of AWTLIST is incorrect.")
temp <- vector("list", nforce)
for (iu in 1:nforce) temp[[iu]] <- awtlist[[1]][[iu]]
awtlist <- temp
} else {
if (length(awtlist) != nforce)
stop("The length of AWTLIST is incorrect.")
}
}

#  check to see if there is anything to estimate

if (difeorder == 0 && nforce == 0)
stop("There are no coefficient functions to estimate.")

ncurve  <- dim(xcoef)[2]
nbasmax <- xbasis\$nbasis

#  check UFDLIST and AWTLIST

if (nforce > 0) {
errorwrd <- FALSE
for (iu in 1:nforce) {
if (!inherits(ufdlist[[iu]], "fd")) {
print(paste("UFDLIST[[",iu,
"]] is not a functional data object.",sep=""))
errorwrd <- TRUE
} else {
ufdi   <- ufdlist[[iu]]
urange <- ufdi\$basis\$rangeval
#  check that urange is equal to xrange
if (any(urange != xrange)) {
print(paste(
"XRANGE and URANGE are not identical for UFDLIST[[",
iu,"]].",sep=""))
errorwrd <- TRUE
}
}
afdPari <- awtlist[[iu]]
afdi    <- afdPari\$fd
if (!inherits(afdi, "fd")) {
print(paste(
"AFDI is not a functional data object for AWTLIST[[",
iu,"]].",sep=""))
errorwrd <- TRUE
} else {
basisi <- afdi\$basis
if (any(basisi\$rangeval != urange)) {
print(paste("Ranges are incompatible for AWTLIST[[",
iu,"]].",sep=""))
errorwrd <- TRUE
}
nbasmax <- max(c(nbasmax,basisi\$nbasis))
}
}
if (errorwrd) stop("")
}

#  check BWTLIST

#  convert to a single-layer list if necessary

if (inherits(bwtlist[[1]],"list")) {
temp <- vector("list",difeorder)
for (j in 1:nvar) {
if (inherits(bwtlist[[1]][[j]], "list")) {
bwtlist[[1]][[j]] <- bwtlist[[1]][[j]][[1]]
}
temp[[j]] <- bwtlist[[1]][[j]]
}
bwtlist <- temp
}

#  check the components
errorwrd <- FALSE
for (j in 1:difeorder) {
if (!is.null(bwtlist[[j]])) {
bfdParj <- bwtlist[[j]]
if (!inherits(bfdParj,"fdPar")) {
print(paste(
"BWTLIST[[",j,"]] is not a functional parameter object.",sep=""))
errorwrd <- TRUE
}  else {
bfdj <- bfdParj\$fd
if (!inherits(bfdj, "fd")) {
print(paste(
"BFDJ in BWTLIST[[",j,"]] is not a functional data object.",
sep=""))
errorwrd <- TRUE
} else {
basisj <- bfdj\$basis
if (any(basisj\$rangeval != xrange)) print(paste(
"Ranges are incompatible for BWTLIST[[",j,"]].",sep=""))
}
}
nbasmax <- max(c(nbasmax,basisj\$nbasis))
}
}
if (errorwrd) stop("")

#  Set up sampling values to be used in numerical integration
#    and set up matrix of basis values.  The number of sampling
#  NFINE is here set to a usually workable value if too small.

if (nfine < 5*nbasmax) nfine <- 5*nbasmax

deltax <- (xrange[2]-xrange[1])/(nfine-1)
tx     <- seq(xrange[1],xrange[2],deltax)

#  set up  YARRAY to hold values of x functions and their derivatives

yarray <- array(0,c(nfine,ncurve,difeordp1))
for (j in 1:difeordp1) yarray[,,j] <- eval.fd(tx, xfdobj, j-1, returnMatrix)

#  set up  UARRAY to hold values of u functions

if (nforce > 0) {
uarray <- array(0,c(nfine,ncurve,nforce))
for (iu in 1:nforce)
uarray[,,iu] <- eval.fd(tx, ufdlist[[iu]], returnMatrix=returnMatrix)
}

#  set up array YPROD to hold mean of products of values in YARRAY

yprod <- array(0,c(nfine,difeordp1,difeordp1))
for (j1 in 1:difeordp1) for (j2 in 1:j1) {
if (ncurve == 1) yprodval <- yarray[,1,j1]*yarray[,1,j2]
else             yprodval <- apply(yarray[,,j1]*yarray[,,j2],1,mean)
yprod[,j1,j2] <- yprodval
yprod[,j2,j1] <- yprodval
}

#  set up array YUPROD to hold mean of u-variables u times
#    x functions and their derivatives

if (nforce > 0) {
yuprod <- array(0,c(nfine, nforce, difeordp1))
for (iu in 1:nforce) {
for (j1 in 1:difeordp1) {
if (ncurve == 1) {
yuprodval <- yarray[,1,j1]*uarray[,1,iu]
} else {
yuprodval <- apply(yarray[,,j1]*uarray[,,iu],1,mean)
}
yuprod[,iu,j1] <- yuprodval
}
}
}

#  set up array UPROD to hold mean of products of u-variables u

if (nforce > 0) {
uprod <- array(0,c(nfine, nforce, nforce))
for (iu in 1:nforce) for (ju in 1:iu) {
if (ncurve == 1) uprodval <- uarray[,1,iu]*uarray[,1,ju]
else             uprodval <- apply(uarray[,,iu]*uarray[,,ju],1,mean)
uprod[,iu,ju] <- uprodval
uprod[,ju,iu] <- uprodval
}
}

#  set up an index array and some arrays of 1's

onesn <- rep(1,nfine)

#  set up array to hold coefficients for basis expansions

if (nforce > 0) {
aarray <- matrix(0,nfine,nforce)
} else {
aarray <- NULL
}

barray <- matrix(0,nfine,difeorder)

#  --------------  beginning of loop through variables  -------------------

#  get number of coefficients to be estimated for this equation

# loop through u-variables

neqns  <- 0

if (nforce > 0) {
for (iu in 1:nforce) {
if (!is.null(awtlist[[iu]])) {
afdPari <- awtlist[[iu]]
if (afdPari\$estimate)
neqns <- neqns + afdPari\$fd\$basis\$nbasis
}
}
}

# loop through x functions and their derivatives

for (j1 in 1:difeorder) {
if (!is.null(bwtlist[[j1]])) {
bfdParj <- bwtlist[[j1]]
if (bfdParj\$estimate)
neqns <- neqns + bfdParj\$fd\$basis\$nbasis
}
}

if (neqns < 1) stop(
"Number of equations to solve is not positive.")

#  set up coefficient array and right side array for linear equation

cmat   <- matrix(0,neqns, neqns)
dmat   <- matrix(0,neqns, 1)

#  evaluate default weight functions for this variable

if (nforce > 0) {
for (iu in 1:nforce) {
if (!is.null(awtlist[[iu]])) {
afdPari     <- awtlist[[iu]]
aarray[,iu] <- eval.fd(tx, afdPari\$fd,returnMatrix=returnMatrix)
}
}
}

for (j1 in 1:difeorder) {
if (!is.null(bwtlist[[j1]])) {
bfdParj     <- bwtlist[[j1]]
barray[,j1] <- eval.fd(tx, bfdParj\$fd, returnMatrix=returnMatrix)
}
}

#  loop through equations,
#    corresponding to rows for CMAT and DMAT

#  loop through equations for u-variables

mi12 <- 0
if (nforce > 0) {
for (iu1 in 1:nforce) {
if (!is.null(awtlist[[iu1]])) {
afdPari1   <- awtlist[[iu1]]
if (afdPari1\$estimate) {
abasisi1    <- afdPari1\$fd\$basis
abasismati1 <- getbasismatrix(tx, abasisi1,
returnMatrix=returnMatrix)
mi11 <- mi12 + 1
mi12 <- mi12 + abasisi1\$nbasis
indexi1 <- mi11:mi12
#  DMAT entry for u-variable
weighti1 <- -yuprod[,iu1,difeordp1]
dmat[indexi1] <-
trapzmat(abasismati1,onesn,deltax,weighti1)
# add terms corresponding to x-derivate weights
# that are not estimated
for (j1 in 1:difeorder) {
bfdParij <- bwtlist[[j1]]
if (!bfdParij\$estimate) {
weightij <- -yuprod[,iu1,j1]
dmat[indexi1] <- dmat[indexi1] +
trapzmat(abasismati1, barray[,j1],
deltax, weightij)
}
}
#  loop through weight functions to be estimated,
#    corresponding to columns for CMAT
#  begin with u-variables
mi22 <- 0
for (iu2 in 1:nforce) {
if (!is.null(awtlist[[iu2]])) {
afdPari2   <- awtlist[[iu2]]
if (afdPari2\$estimate) {
abasisi2    <- afdPari2\$fd\$basis
abasismati2 <- getbasismatrix(tx, abasisi2,
returnMatrix=returnMatrix)
weighti2    <- uprod[,iu1,iu2]
Cprod  <- trapzmat(abasismati1, abasismati2,
deltax, weighti2)
mi21 <- mi22 + 1
mi22 <- mi22 + abasisi2\$nbasis
indexi2 <- mi21:mi22
#  coefficient matrix CMAT entry
cmat[indexi1,indexi2] <- Cprod
}
}
}
#  remaining columns:
#    loop through u-variable -- x-derivative pairs
mij22 <- mi22
for (j2 in 1:difeorder) {
if (!is.null(bwtlist[[j2]])) {
bfdParj2     <- bwtlist[[j2]]
if (bfdParj2\$estimate) {
bbasisij2    <- bfdParj2\$fd\$basis
bbasismatij2 <- getbasismatrix(tx, bbasisij2,
returnMatrix=returnMatrix)
weightij12   <- -yuprod[,iu1,j2]
Cprod <- trapzmat(abasismati1,bbasismatij2,
deltax,weightij12)
mij21 <- mij22 + 1
mij22 <- mij22 + bbasisij2\$nbasis
indexij2  <- mij21:mij22
cmat[indexi1,indexij2] <- Cprod
}
}
}
#  add roughness penalty matrix to diagonal entries
lambdai1 <- afdPari1\$lambda
if (lambdai1 > 0) {
Lfdobj <- afdPari1\$Lfd
penmat <- lambdai1*eval.penalty(abasisi1, Lfdobj)
cmat[indexi1,indexi1] <- cmat[indexi1,indexi1] + penmat
}
}
}
}
}

#  loop through equations for x-derivatives

mij12 <- mi12
for (j1 in 1:difeorder) {
if (!is.null(bwtlist[[j1]])) {
bfdParj1 <- bwtlist[[j1]]
if (bfdParj1\$estimate) {
bbasisij1    <- bfdParj1\$fd\$basis
bbasismatij1 <- getbasismatrix(tx,bbasisij1,
returnMatrix=returnMatrix)
mij11 <- mij12 + 1
mij12 <- mij12 + bbasisij1\$nbasis
indexij1 <- mij11:mij12
#  DMAT entry for u-variable -- x-derivative pair
weightij1 <- yprod[,j1,difeordp1]
dmat[indexij1] <-
trapzmat(bbasismatij1,onesn,deltax,weightij1)
#  add terms corresponding to forcing functions
#  with unestimated coefficients
if (nforce > 0) {
for (iu in 1:nforce) {
if (!is.null(awtlist[[iu]])) {
afdPari <- awtlist[[iu]]
if (!afdPari\$estimate) {
weightijk <- -yuprod[,iu,j1]
dmat[indexij1] <- dmat[indexij1] +
trapzmat(bbasisij1, aarray[,iu],deltax, weightijk)
}
}
}
}
}
}
#  first columns of CMAT: u-variable entries
mi22 <- 0
if (nforce > 0) {
for (iu2 in 1:nforce) {
if (!is.null(awtlist[[iu2]])) {
afdPari2 <- awtlist[[iu2]]
if (afdPari2\$estimate) {
abasisi2    <- afdPari2\$fd\$basis
abasismati2 <- getbasismatrix(tx, abasisi2,
returnMatrix=returnMatrix)
weighti2    <- -yuprod[,iu2,j1]
Cprod <- trapzmat(bbasismatij1,abasismati2,deltax,weighti2)
mi21 <- mi22 + 1
mi22 <- mi22 + abasisi2\$nbasis
indexi2 <- mi21:mi22
cmat[indexij1,indexi2] <- Cprod
}
}
}
}
#  remaining columns: x-derivative pairs
mij22 <- mi22
for (j2 in 1:difeorder) {
if (!is.null(bwtlist[[j2]])) {
bfdParj2  <- bwtlist[[j2]]
bbasisij2 <- bfdParj2\$fd\$basis
if (bfdParj2\$estimate) {
bbasismatij2 <- getbasismatrix(tx, bbasisij2,
returnMatrix=returnMatrix)
weightij22   <- yprod[,j1,j2]
Cprod <- trapzmat(bbasismatij1,bbasismatij2,deltax,weightij22)
mij21 <- mij22 + 1
mij22 <- mij22 + bbasisij2\$nbasis
indexij2 <- mij21:mij22
cmat[indexij1,indexij2] <- Cprod
}
}
}
# add roughness penalty matrix to diagonal entries
lambdaj1 <- bfdParj1\$lambda
if (lambdaj1 > 0) {
Lfdobj <- bfdParj1\$Lfd
penmat <- lambdaj1*eval.penalty(bbasisij1, Lfdobj)
cmat[indexij1,indexij1] <- cmat[indexij1,indexij1] + penmat
}
}

#  --------------  end of loop through variables  -------------------

# solve for coefficients of basis expansions

dvec <- -symsolve(cmat,dmat)

#  set up u-function weight functions

mi2 <- 0
if (nforce > 0) {
for (iu in 1:nforce) {
if (!is.null(awtlist[[iu]])) {
afdPari <- awtlist[[iu]]
if (afdPari\$estimate) {
mi1 <- mi2 + 1
mi2 <- mi2 + afdPari\$fd\$basis\$nbasis
indexi <- mi1:mi2
afdPari\$fd\$coefs <- as.matrix(dvec[indexi])
awtlist[[iu]] <- afdPari
}
}
}
}

#  set up X-function derivative weight functions

mij2 <- mi2
for (j in 1:difeorder) {
if (!is.null(bwtlist[[j]])) {
bfdParj <- bwtlist[[j]]
if (bfdParj\$estimate) {
mij1 <- mij2 + 1
mij2 <- mij2 + bfdParj\$fd\$basis\$nbasis
indexij <- mij1:mij2
bfdParj\$fd\$coefs <- as.matrix(dvec[indexij])
bwtlist[[j]] <- bfdParj
}
}
}

#  set up residual list RESFDLIST

#  initialize with highest order derivative for this variable
resmat  <- eval.fd(tx, xfdobj, difeorder, returnMatrix)
#  add contributions from weighted u-functions
if (nforce > 0) {
onesncurve <- rep(1,ncurve)
for (iu in 1:nforce) {
if (!is.null(awtlist[[iu]])) {
afdPari  <- awtlist[[iu]]
aveci    <- as.vector(eval.fd(tx, afdPari\$fd,
returnMatrix=returnMatrix))
umati    <- eval.fd(tx, ufdlist[[iu]], returnMatrix=returnMatrix)
aumati   <- outer(aveci,onesncurve)*umati
resmat   <- resmat - aumati
}
}
}
#  add contributions from weighted x-function derivatives
for (j in 1:difeorder) {
if (!is.null(bwtlist[[j]])) {
bfdParj <- bwtlist[[j]]
bmatij  <- as.vector(eval.fd(tx, bfdParj\$fd,
returnMatrix=returnMatrix))
xmatij  <- eval.fd(tx, xfdobj, j-1, returnMatrix)
resmat  <- resmat + bmatij*xmatij
}
}
#  set up the functional data object
resbasis <- xbasis
resfd    <- smooth.basis(tx, resmat, resbasis)\$fd
resfdnames      <- xfdobj\$fdnames
resfdnames[[2]] <- "Residual function"
resfdnames[[3]] <- "Residual function value"
resfd\$fdnames   <- resfdnames
resfdlist       <- list(resfd)

#  ----------------------------------------------------------------
#                   End of single variable case
#  ----------------------------------------------------------------

} else {

#  ----------------------------------------------------------------
#                   Multiple variable case
#  ----------------------------------------------------------------

#  check the dimensions of UFDLIST and AWTLIST

if (is.null(ufdlist) || is.null(awtlist)) {
awtlist <- NULL
} else {
if (length(ufdlist) != nvar)
stop(paste("The length of UFDLIST",
" does not match that of XFDLIST."))
errorwrd = FALSE
for (j in 1:nvar) {
if (!is.null(ufdlist[[j]])) {
nforce <- length(ufdlist[[j]])
if (length(awtlist[[j]]) != nforce) {
print(paste("The length of AWTLIST[[",j,
"]] is incorrect.",sep=""))
errorwrd = TRUE
}
}
}
if (errorwrd) stop("")
}

#  check the dimensions of BWTLIST

if (length(bwtlist) != nvar) stop("Length of BWTLIST is incorrect.")
errorwrd = FALSE
for (ivar in 1:nvar) {
if (length(bwtlist[[ivar]]) != nvar) {
print(paste("The length of BWTLIST[[",ivar,
"]] is incorrect.",sep=""))
errorwrd = TRUE
}
}
if (errorwrd) stop("")

#  check XFDLIST and extract NCURVE and XRANGE

xfd1       <- xfdlist[[1]]
xcoef1     <- xfd1\$coefs
xbasis1    <- xfd1\$basis
xrange1    <- xbasis1\$rangeval
ncurve     <- dim(xcoef1)[2]
resfdnames <- xfd1\$fdnames

errorwrd = FALSE
for (ivar in 1:nvar) {
xfdi    <- xfdlist[[ivar]]
xcoefi  <- xfdi\$coefs
xbasisi <- xfdi\$basis
xrangei <- xbasisi\$rangeval
ncurvei <- dim(xcoefi)[2]
if (!inherits(xfdi, "fd")) {
print(paste("XFDLIST[[",ivar,
"]] is not a functional data object.",sep=""))
errorwrd = TRUE
} else {
if (any(xrangei != xrange1)) {
print("Ranges are incompatible for XFDLIST.")
errorwrd = TRUE
}
if (ncurvei != ncurve) {
print("Number of curves is incompatible for XFDLIST.")
errorwrd = TRUE
}
}
}
if (errorwrd) stop("")

nbasmax <- xbasis1\$nbasis
#  This will be the maximum number of basis functions

#  check compatibility of UFDLIST and AWTLIST

if (!(is.null(ufdlist) || is.null(awtlist))) {
urange <- ufdlist[[1]]\$basis\$rangeval
errorwrd <- FALSE
for (ivar in 1:nvar) {
if (!is.null(ufdlist[[ivar]])) {
for (iu in 1:length(ufdlist[[ivar]])) {
ufdiviu <- ufdlist[[ivar]][[iu]]
if (!inherits(ufdiviu, "fd")) {
print(paste("UFDLIST[[",ivar,",",iu,
"]] is not a functional data object.",
sep=""))
errorwrd <- TRUE
}
if (any(ufdiviu\$basis\$rangeval != urange)) {
print("Ranges are incompatible for UFDLIST.")
errorwrd <- TRUE
}
awtfdPari <- awtlist[[ivar]][[iu]]
if (!inherits(awtfdPari, "fdPar")) {
print(paste("AWTFDPAR[[",ivar,"]][[",iu,
"]] is not a functional parameter object.",sep=""))
errorwrd <- TRUE
}
afdi   <- awtfdPari\$fd
basisi <- afdi\$basis
if (any(basisi\$rangeval != urange)) {
print("Ranges are incompatible for AWTLIST.")
errorwrd <- TRUE
}
nbasmax <- max(c(nbasmax,basisi\$nbasis))
}
if (errorwrd) stop("")
}
}
}

#  check BWTLIST

errorwrd <- FALSE
for (ivar1 in 1:nvar) {
for (ivar2 in 1:nvar) {
difeorder <- length(bwtlist[[ivar1]][[ivar2]])
for (j in 1:difeorder) {
if (!is.null(bwtlist[[ivar1]][[ivar2]][[j]])) {
bfdPari1i2j <- bwtlist[[ivar1]][[ivar2]][[j]]
if (!inherits(bfdPari1i2j, "fdPar")) {
print(paste("BWTLIST[[",ivar1, ",",ivar2, ",",j,
"]] is not a functional parameter object.",sep=""))
errorwrd = TRUE
}
basisi1i2j <- bfdPari1i2j\$fd\$basis
if (any(basisi1i2j\$rangeval != xrange1)) {
print(paste("Ranges are incompatible for BWTLIST[[",
ivar1,"]][[",ivar2,"]][[",
j,"]]",sep=""))
errorwrd <- TRUE
}
nbasmax <- max(c(nbasmax,basisi1i2j\$nbasis))
}
}
}
}
if (errorwrd) stop("")

#  set up sampling values to be used in numerical integration
#    and set up matrix of basis values.  The number of sampling
#  NFINE is here set to a usually workable value if too small.

if (nfine < 5*nbasmax) nfine <- 5*nbasmax

deltax <- (xrange1[2]-xrange1[1])/(nfine-1)
tx     <- seq(xrange1[1],xrange1[2],deltax)

#  set up  YARRAY to hold values of x functions and their derivatives

yarray <- vector("list", 0)
for (ivar in 1:nvar) {
difeorder <- length(bwtlist[[ivar]][[ivar]])
difeordp1 <- difeorder + 1
yarray[[ivar]] <- array(0,c(nfine,ncurve,difeordp1))
for (j in 1:difeordp1){
yj <- eval.fd(tx, xfdlist[[ivar]], j-1, returnMatrix=returnMatrix)
yarray[[ivar]][,,j] <- as.matrix(yj)
}
}

#  set up  UARRAY to hold values of u functions

if (!is.null(ufdlist)) {
uarray <- vector("list", nvar)
for (ivar in 1:nvar) {
if (is.null(ufdlist[[ivar]])) {
uarray[[ivar]] <- NULL
} else {
nforce <- length(ufdlist[[ivar]])
uarray[[ivar]] <- vector("list", nforce)
for (iu in 1:nforce)
uarray[[ivar]][[iu]] <- matrix(0,nfine,ncurve)
}
}
for (ivar in 1:nvar) {
if (!is.null(ufdlist[[ivar]])) {
nforce <- length(ufdlist[[ivar]])
for (iu in 1:nforce)
uarray[[ivar]][[iu]] <- eval.fd(tx, ufdlist[[ivar]][[iu]],
returnMatrix=returnMatrix)
}
}
}

#  set up array YPROD to hold mean of products of values in YARRAY

yprod <- vector("list", nvar)
for (i1 in 1:nvar) yprod[[i1]] <- vector("list", nvar)
for (i1 in 1:nvar) {
difeord1p1 <- length(bwtlist[[i1]][[i1]]) + 1
for (i2 in 1:nvar) {
difeord2p1 <- length(bwtlist[[i2]][[i2]]) + 1
yprod[[i1]][[i2]] <- array(0,c(nfine,difeord2p1,difeord2p1))
}
}

for (i1 in 1:nvar) {
difeord1p1 <- length(bwtlist[[i1]][[i1]]) + 1
for (j1 in 1:difeordp1) {
for (i2 in 1:nvar) {
difeord2p1 <- length(bwtlist[[i2]][[i2]]) + 1
for (j2 in 1:difeord2p1) {
if (ncurve == 1) {
yprodval <-       yarray[[i1]][,1,j1]*yarray[[i2]][,1,j2]
} else {
yprodval <- apply(yarray[[i1]][,,j1]*yarray[[i2]][,,j2],1,mean)
}
yprod[[i1]][[i2]][,j1,j2] <- yprodval
}
}
}
}

#  set up array YUPROD to hold mean of u-variables u times
#    x functions and their derivatives

if (!is.null(ufdlist)) {
yuprod <- vector("list", nvar)
for (i1 in 1:nvar) {
if (!is.null(ufdlist[[i1]])) {
nforce <- length(ufdlist[[i1]])
if (nforce > 0) {
yuprod[[i1]] <- vector("list", nforce)
for (iu in 1:nforce) {
difeordp1 <- length(bwtlist[[i1]][[i1]]) + 1
yuprod[[i1]][[iu]] <- matrix(0,nfine,difeordp1)
}
}
}
}
onesncurve <- rep(1,ncurve)
for (i1 in 1:nvar) {
if (!is.null(ufdlist[[i1]])) {
nforce <- length(ufdlist[[i1]])
if (nforce > 0) {
difeordp1 <- length(bwtlist[[i1]][[i1]]) + 1
for (iu in 1:nforce) {
for (j1 in 1:difeordp1) {
if (ncurve == 1) {
yuprodval <- yarray[[i1]][,1,j1]*uarray[[i1]][[iu]]
} else {
yuprodval <- apply(yarray[[i1]][,,j1]*
outer(uarray[[i1]][[iu]],onesncurve),1,mean)
}
yuprod[[i1]][[iu]][,j1] <- yuprodval
}
}
}
}
}
}

#  set up array UPROD to hold mean of products of u-variables u

if (!is.null(ufdlist)) {
uprod <- vector("list", nvar)
for (ivar in 1:nvar) {
nforce <- length(ufdlist[[ivar]])
if (nforce > 0) {
uprod[[ivar]] <- array(0,c(nfine, nforce, nforce))
for (iu in 1:nforce) for (ju in 1:iu) {
uprodval <- uarray[[ivar]][[iu]]*uarray[[ivar]][[ju]]
uprod[[ivar]][,iu,ju] <- uprodval
uprod[[ivar]][,ju,iu] <- uprodval
}
}
}
}

#  set up an index array and some arrays of 1"s

onesn <- rep(1,nfine)

#  set up array to hold coefficients for basis expansions

#  --------------  beginning of loop through variables  -------------------

for (ivar in 1:nvar) {

#  get number of coefficients to be estimated for this equation

neqns  <- 0

# loop through u-variables  if required

if (is.null(ufdlist) || is.null(ufdlist[[ivar]])) {
nforce <- 0
} else {
nforce <- length(ufdlist[[ivar]])
if (nforce > 0) {
for (iu in 1:nforce) {
afdPari <- awtlist[[ivar]][[iu]]
if (afdPari\$estimate)
neqns <- neqns + afdPari\$fd\$basis\$nbasis
}
}
}

# loop through x functions and their derivatives

for (i2 in 1:nvar) {
difeorder <- length(bwtlist[[ivar]][[i2]])
for (j2 in 1:difeorder) {
if (!is.null(bwtlist[[ivar]][[i2]][[j2]])) {
bfdParij <- bwtlist[[ivar]][[i2]][[j2]]
if (bfdParij\$estimate) neqns <- neqns + bfdParij\$fd\$basis\$nbasis
}
}
}
if (neqns < 1)  stop("Number of equations to solve is not positive.")

#  set up coefficient array and right side array for linear equation

cmat   <- matrix(0,neqns, neqns)
dmat   <- matrix(0,neqns, 1)

#  evaluate default weight functions for this variable

if (nforce > 0) {
aarray <- matrix(0,nfine,nforce)
for (iu in 1:nforce) {
if (!is.null(awtlist[[ivar]][[iu]])) {
afdPari <- awtlist[[ivar]][[iu]]
aarray[,iu] <- eval.fd(tx, afdPari\$fd, returnMatrix=returnMatrix)
}
}
}

barray <- vector("list", nvar)
for (i in 1:nvar) {
difeorder <- length(bwtlist[[ivar]][[i]])
barray[[i]] <- array(0,c(nfine,nvar,difeorder))
for (j in 1:difeorder) {
if (!is.null(bwtlist[[ivar]][[i]][[j]])) {
bfdParij <- bwtlist[[ivar]][[i]][[j]]
bfdij    <- bfdParij\$fd
barray[[i]][,,j] <- as.matrix(eval.fd(tx, bfdij),
returnMatrix=returnMatrix)
}
}
}

#  loop through equations,
#    corresponding to rows for CMAT and DMAT

#  loop through equations for u-variables

mi12 <- 0
if (nforce > 0) {
for (iu1 in 1:nforce) {
if (!is.null(awtlist[[ivar]][[iu1]])) {
afdPari1   <- awtlist[[ivar]][[iu1]]
if (afdPari1\$estimate) {
abasisi1    <- afdPari1\$fd\$basis
abasismati1 <- getbasismatrix(tx, abasisi1, returnMatrix=returnMatrix)
mi11 <- mi12 + 1
mi12 <- mi12 + abasisi1\$nbasis
indexi1 <- mi11:mi12
#  DMAT entry for u-variable
weighti1 <- -yuprod[[ivar]][[iu1]][,difeordp1]
dmat[indexi1] <- trapzmat(abasismati1,onesn,deltax,weighti1)
#  add terms corresponding to x-derivative weights
#  that are not estimated
for (i in 1:nvar) {
difeorder <- length(bwtlist[[ivar]][[i]])
for (j in 1:difeorder) {
bfdParij <- bwtlist[[ivar]][[i]][[j]]
if (!bfdParij\$estimate) {
weightij <- -yuprod[[ivar]][[iu1]][,j]
dmat[indexi1] <- dmat[indexi1] +
trapzmat(abasismati1, barray[[ivar]][,,j],
deltax, weightij)
}
}
}
#  loop through weight functions to be estimated,
#    corresponding to columns for CMAT
#  begin with u-variables
mi22 <- 0
for (iu2 in 1:nforce) {
if (!is.null(awtlist[[ivar]][[iu2]])) {
afdPari2   <- awtlist[[ivar]][[iu2]]
if (afdPari2\$estimate) {
abasisi2    <- afdPari2\$fd\$basis
abasismati2 <- getbasismatrix(tx, abasisi2,
returnMatrix=returnMatrix)
weighti2    <- uprod[[ivar]][,iu1,iu2]
Cprod       <- trapzmat(abasismati1, abasismati2,
deltax, weighti2)
mi21 <- mi22 + 1
mi22 <- mi22 + abasisi2\$nbasis
indexi2 <- mi21:mi22
#  coefficient matrix CMAT entry
cmat[indexi1,indexi2] <- Cprod
}
}
}
#  remaining columns:
#    loop through u-variable -- x-derivative pairs
mij22 <- mi22
for (i2 in 1:nvar) {
if (!is.null(bwtlist[[ivar]][[i2]])) {
difeorder <- length(bwtlist[[ivar]][[i2]])
for (j2 in 1:difeorder) {
bfdParij2   <- bwtlist[[ivar]][[i2]][[j2]]
if (bfdParij2\$estimate) {
bbasisij2    <- bfdParij2\$fd\$basis
bbasismatij2 <- getbasismatrix(tx, bbasisij2,
returnMatrix=returnMatrix)
weightij12   <- -yuprod[[i2]][[iu1]][,j2]
Cprod        <- trapzmat(abasismati1,bbasismatij2,
deltax,weightij12)
mij21 <- mij22 + 1
mij22 <- mij22 + bbasisij2\$nbasis
indexij2  <- mij21:mij22
cmat[indexi1,indexij2] <- Cprod
}
}
}
}
#  add roughness penalty matrix to diagonal entries
lambdai1 <- afdPari1\$lambda
if (lambdai1 > 0) {
Lfdobj <- afdPari1\$Lfd
penmat <- lambdai1*eval.penalty(abasisi1,Lfdobj)
cmat[indexi1,indexi1] <- cmat[indexi1,indexi1] + penmat
}
}
}
}
}

#  end of loop through equations for u-variables

#  loop through equations for x-derivatives

mij12 <- mi12
for (i1 in 1:nvar) {
difeorder1 <- length(bwtlist[[ivar]][[i1]])
for (j1 in 1:difeorder1) {
if (!is.null(bwtlist[[ivar]][[i1]][[j1]])) {
bfdParij1 <- bwtlist[[ivar]][[i1]][[j1]]
if (bfdParij1\$estimate) {
bbasisij1    <- bfdParij1\$fd\$basis
bbasismatij1 <- getbasismatrix(tx, bbasisij1, returnMatrix=returnMatrix)
mij11 <- mij12 + 1
mij12 <- mij12 + bbasisij1\$nbasis
indexij1 <- mij11:mij12
#  DMAT entry for u-variable -- x-derivative pair
weightij1 <- yprod[[i1]][[ivar]][,j1,difeordp1]
dmat[indexij1] <- trapzmat(bbasismatij1,onesn,
deltax,weightij1)
#  add terms corresponding to forcing functions
#  with unestimated coefficients
if (nforce > 0) {
for (iu in 1:nforce) {
if (!is.null(awtlist[[ivar]][[iu]])) {
afdPari <- awtlist[[ivar]][[iu]]
if (!afdPari\$estimate) {
weightijk <- yprod[,ivar,iu,j1]
dmat[indexij1] <- dmat[indexij1] +
trapzmat(bbasismatij1,aarray[,iu],deltax,weightijk)
}
}
}
}
#  first columns of CMAT: u-variable entries
mi22 <- 0
if (nforce > 0) {
for (iu2 in 1:nforce) {
if (!is.null(awtlist[[ivar]][[iu2]])) {
afdPari2  <- awtlist[[ivar]][[iu2]]
if (afdPari2\$estimate) {
abasisi2    <- afdPari2\$fd\$basis
abasismati2 <- getbasismatrix(tx, abasisi2, returnMatrix=returnMatrix)
weighti2    <- -yuprod[[i1]][[iu2]][,j1]
Cprod <- trapzmat(bbasismatij1,abasismati2,deltax,weighti2)
mi21 <- mi22 + 1
mi22 <- mi22 + abasisi2\$nbasis
indexi2 <- mi21:mi22
cmat[indexij1,indexi2] <- Cprod
}
}
}
}
#  remaining columns: x-derivative pairs
mij22 <- mi22
for (i2 in 1:nvar) {
difeorder2 <- length(bwtlist[[ivar]][[i2]])
for (j2 in 1:difeorder2) {
if (!is.null(bwtlist[[ivar]][[i2]][[j2]])) {
bfdParij2 <- bwtlist[[ivar]][[i2]][[j2]]
bbasisij2    <- bfdParij2\$fd\$basis
bbasismatij2 <- getbasismatrix(tx, bbasisij2, returnMatrix=returnMatrix)
weightij22   <- yprod[[i1]][[i2]][,j1,j2]
Cprod <- trapzmat(bbasismatij1,bbasismatij2,deltax,weightij22)
if (bfdParij2\$estimate) {
mij21 <- mij22 + 1
mij22 <- mij22 + bbasisij2\$nbasis
indexij2 <- mij21:mij22
cmat[indexij1,indexij2] <- Cprod
}
}
}
}
#  add roughness penalty terms to diagonal entries
lambdaij1 <- bfdParij1\$lambda
if (lambdaij1 > 0) {
Lfdobj <- bfdParij1\$Lfd
penmat <- lambdaij1*eval.penalty(bbasisij1,Lfdobj)
cmat[indexij1,indexij1] <- cmat[indexij1,indexij1] +penmat
}
}
}
}
}

#  end of loop through x derivatives

dvec <- -symsolve(cmat,dmat)

#  set up u-function weight functions

mi2 <- 0
if (nforce > 0) {
for (iu in 1:nforce) {
if (!is.null(awtlist[[ivar]][[iu]])) {
afdPari <- awtlist[[ivar]][[iu]]
if (afdPari\$estimate) {
mi1 <- mi2 + 1
mi2 <- mi2 + afdPari\$fd\$basis\$nbasis
indexi <- mi1:mi2
afdPari\$fd\$coefs <- as.matrix(dvec[indexi])
awtlist[[ivar]][[iu]] <- afdPari
}
}
}
}

#  set up X-function derivative weight functions

mij2 <- mi2
for (i1 in 1:nvar) {
difeorder <- length(bwtlist[[ivar]][[i1]])
for (j1 in 1:difeorder) {
if (!is.null(bwtlist[[ivar]][[i1]][[j1]])) {
bfdParij <- bwtlist[[ivar]][[i1]][[j1]]
if (bfdParij\$estimate) {
mij1 <- mij2 + 1
mij2 <- mij2 + bfdParij\$fd\$basis\$nbasis
indexij <- mij1:mij2
bfdParij\$fd\$coefs <- as.matrix(dvec[indexij])
bwtlist[[ivar]][[i1]][[j1]] <- bfdParij
}
}
}
}
}

#  --------------  end of loop through variables  -------------------

#  set up residual list RESFDLIST

resfdlist <- vector("list", nvar)

for (ivar in 1:nvar) {
difeorder <- length(bwtlist[[ivar]][[ivar]])
xfdi      <- xfdlist[[ivar]]
resbasis  <- xfdi\$basis
#  initialize with highest order derivative for this variable
resmat    <- eval.fd(tx, xfdi, difeorder, returnMatrix)
#  add contributions from weighted u-functions
onesncurve <- rep(1,ncurve)
if (!is.null(ufdlist)) {
nforce <- length(ufdlist[[ivar]])
if (nforce > 0) {
for (iu in 1:nforce) {
if (!is.null(awtlist[[ivar]][[iu]])) {
afdPari  <- awtlist[[ivar]][[iu]]
amati    <- as.vector(eval.fd(tx, afdPari\$fd,
returnMatrix=returnMatrix))
umati    <- eval.fd(tx, ufdlist[[ivar]][[iu]],
returnMatrix=returnMatrix)
if (ncurve == 1) {
aumati <- amati*umati
} else {
aumati <- outer(amati,onesncurve)*umati
}
resmat   <- resmat - aumati
}
}
}
}
#  add contributions from weighted x-function derivatives
for (i1 in 1:nvar) {
difeorder <- length(bwtlist[[ivar]][[i1]])
for (j1 in 1:difeorder) {
if (!is.null(bwtlist[[ivar]][[i1]][[j1]])) {
bfdParij <- bwtlist[[ivar]][[i1]][[j1]]
bfdij    <- bfdParij\$fd
bvecij   <- as.vector(eval.fd(tx, bfdij, returnMatrix=returnMatrix))
if (ncurve == 1) {
bmatij <- bvecij
}  else  {
bmatij <- outer(bvecij,onesncurve)
}
xmatij <- eval.fd(tx, xfdlist[[i1]], j1-1, returnMatrix)
resmat <- resmat + bmatij*xmatij
}
}
}
#  set up the functional data object
resfdi            <- smooth.basis(tx, resmat, resbasis)\$fd
resfdnames        <- xfdi\$fdnames
resfdnames[[2]]   <- "Residual function"
resfdnames[[3]]   <- "Residual function value"
resfdlist[[ivar]] <- resfdi
}

#  ----------------------------------------------------------------
#                   End of multiple variable case
#  ----------------------------------------------------------------

}

pdaList <- list(bwtlist=bwtlist, resfdlist=resfdlist, awtlist=awtlist)
class(pdaList) <- 'pda.fd'
pdaList
}

# ---------------------------------------------------------------------------

trapzmat <- function(X,Y,delta=1,wt=rep(1,n)) {
#TRAPZMAT integrates the products of two matrices of values
#   using the trapezoidal rule, assuming equal spacing
#  X is the first  matrix of values
#  Y is the second matrix of values
#  DELTA is the spacing between argument values (one by default)
#  WT is a vector of weights (ones by default)
#
#  XtWY is a matrix of integral estimates, number of rows equal to
#  number of col of X, number of cols equal to number of cols of Y

X <- as.matrix(X)
Y <- as.matrix(Y)

n <- dim(X)[1]

if (dim(Y)[1] != n) {
stop("X and Y do not have same number of rows.")
}

if (length(wt) != n) {
stop("X and WT do not have same number of rows.")
}

if (delta <= 0) {
stop("DELTA is not a positive value.")
}

wt[c(1,n)] <- wt[c(1,n)]/2
wt <- wt*delta

X <- X*outer(wt,rep(1,dim(X)[2]))
XtWY <- crossprod(X,Y)
return(XtWY)
}
```

Try the fda package in your browser

Any scripts or data that you put into this service are public.

fda documentation built on May 2, 2019, 5:12 p.m.