Nothing
R/pda.fd.R
In fda: Functional Data Analysis
Defines functions
pda.fd
trapzmat
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.
# last modified 16 April 2014 by Jim Ramsay
# 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.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions pda.fd trapzmat
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.
# last modified 16 April 2014 by Jim Ramsay
# 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.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.