# R/basisfd.R In fda: Functional Data Analysis

#### Documented in basisfdsummary.basisfd

```#  Generator function of class basisfd

basisfd <- function(type, rangeval, nbasis, params, dropind=vector("list",0),
basisvalues=vector("list",0))
{
#  BASISFD  generator function of "basisfd" class.
#  Arguments:
#  TYPE    ...a string indicating the type of basisobj.
#             This may be one of:
#             "Bspline", "bspline", "Bsp", "bsp",
#             "con", "const", "constant"
#             "exp", "exponen", "exponential"
#             "fdVariance"
#             "FEM"
#             "Fourier", "fourier", "Fou", "fou",
#             "mon", "monom", "monomial",
#             "polyg", "polygon", "polygonal"
#             "power" "pow"
#  RANGEVAL...an array of length 2 containing the lower and upper
#             boundaries for (the rangeval of argument values
#             If basis is of FEM type, rangeval is not used
#  NBASIS ... the number of basis functions
#  PARAMS ... If the basis is "fourier", this is a single number indicating
#               the period.  That is, the basis functions are periodic on
#               the interval (0,PARAMS) or any translation of it.
#             If the basis is "bspline", the values are interior points at
#               which the piecewise polynomials join.
#               Note that the number of basis functions NBASIS is equal
#               to the order of the Bspline functions plus the number of
#               interior knots, that is the length of PARAMS.
#             This means that NBASIS must be at least 1 larger than the
#               length of PARAMS.
#  DROPIND...A set of indices in 1:NBASIS of basis functions to drop when
#              basis objects are arguments.  Default is vector("list",0)
#              Note that argument NBASIS is reduced by the number of
#              indices, and the derivative matrices in VALUES are also clipped.
#              points to be used in a fixed point quadrature.  The second
#              contains quadrature weights.  For example, for (Simpson"s
#              rule for (NQUAD = 7, the points are equally spaced and the
#              weights are delta.*[1, 4, 2, 4, 2, 4, 1]/3.  DELTA is the
#              spacing between quadrature points.  The default is
#              matrix("numeric",0,0).
#  VALUES ...A list, with entries containing the values of
#              the basis function derivatives starting with 0 and
#              going up to the highest derivative needed.  The values
#              up to the user to decide whether or not to multiply
#              the derivative values by the square roots of the
#              quadrature weights so as to make numerical integration
#              a simple matrix multiplication.
#              Values are checked against QUADVALS to ensure the correct
#              number of rows, and against NBASIS to ensure the correct
#              number of columns.
#              The default value of is VALUES is vector("list",0).
#              VALUES contains values of basis functions and derivatives at
#              These values are only generated as required, and only if slot
#  BASISVALUES...A vector of lists, allocated by code such as
#              vector("list",1).
#              This field is designed to avoid evaluation of a
#              basis system repeatedly at a set of argument values.
#              Each list within the vector corresponds to a specific set
#              of argument values, and must have at least two components,
#              which may be tagged as you wish.
#              The first component in an element of the list vector contains the
#              argument values.
#              The second component in an element of the list vector
#              contains a matrix of values of the basis functions evaluated
#              at the arguments in the first component.
#              The third and subsequent components, if present, contain
#              matrices of values their derivatives up to a maximum
#              derivative order.
#              Whenever function getbasismatrix is called, it checks
#              the first list in each row to see, first, if the number of
#              argument values corresponds to the size of the first dimension,
#              and if this test succeeds, checks that all of the argument
#              values match.  This takes time, of course, but is much
#              faster than re-evaluation of the basis system.  Even this
#              time can be avoided by direct retrieval of the desired
#              array.
#              For example, you might set up a vector of argument values
#              called "evalargs" along with a matrix of basis function
#              values for these argument values called "basismat".
#              You might want too use tags like "args" and "values",
#              respectively for these.  You would then assign them
#              to BASISVALUES with code such as
#                basisobj\$basisvalues <- vector("list",1)
#                basisobj\$basisvalues[[1]] <-
#                             list(args=evalargs, values=basismat)
#
#  Returns
#  BASISOBJ  ... a basisfd object with slots
#         type
#         rangeval
#         nbasis
#         params
#         dropind
#         values
#         basisvalues
#  Slot VALUES contains values of basis functions and derivatives at
#   These values are only generated as required, and only if slot
#
#  An alternative name for (this function is CREATE.BASIS, but PARAMS argument
#     must be supplied.
#  Specific types of bases may be set up more conveniently using functions
#  CREATE.BSPLINE.BASIS     ...  creates a b-spline basis
#  CREATE.CONSTANT.BASIS    ...  creates a constant basis
#  CREATE.EXPONENTIAL.BASIS ...  creates an exponential basis
#  CREATE.FDVARIANCE.BASIS  ...  creates an fdVariance basis
#  CREATE.FEM.BASIS         ...  creates an FEM basis
#  CREATE.FOURIER.BASIS     ...  creates a fourier basis
#  CREATE.MONOMIAL.BASIS    ...  creates a monomial basis
#  CREATE.POLYGON.BASIS     ...  creates a polygonal basis
#  CREATE.POLYNOMIAL.BASIS  ...  creates a polynomial basis
#  CREATE.POWER.BASIS       ...  creates a monomial basis

# value -> values 2012.12.27 by spencer graves

#  Set up default basis if there are no arguments:
#     order 2 monomial basis over [0,1]

if (nargs()==0) {
type        <- "bspline"
rangeval    <- c(0,1)
nbasis      <- 2
params      <- vector("list",0)
dropind     <- vector("list",0)
values      <- vector("list",0)
basisvalues <- vector("list",0)

basisobj  <- list(type=type,     rangeval=rangeval, nbasis=nbasis,
values=values, basisvalues=basisvalues)
oldClass(basisobj) <- "basisfd"
return(basisobj)
}

#  if first argument is a basis object, return

if (inherits(type,"basisfd")) {
basisobj <- type
return(basisobj)
}

#  check basistype

# type <- moreNames(type)

#  recognize type of basis by use of several variant spellings

if(type == "bspline" ||
type == "Bspline" ||
type == "spline"  ||
type == "Bsp"     ||
type == "bsp") {
type = "bspline"
}
else if(type == "con"      ||
type == "const"    ||
type == "constant") {
type = "const"
}
else if(type == "exp"    ||
type == "expon"  ||
type == "exponential") {
type = "expon"
}
else if(type == "fdVariance" ||
type == "fdVar") {
type = "fdVariance"
}
else if(type == "FEM") {
type = "FEM"
}
else if(type == "Fourier" ||
type == "fourier" ||
type == "Fou"     ||
type == "fou") {
type = "fourier"
}
else if(type == "mon" ||
type == "monom"  ||
type == "monomial") {
type = "monom"
}
else if(type == "polyg"    ||
type == "polygon"  ||
type == "polygonal") {
type = "polygonal"
}
else if(type == "polynomial"    ||
type == "polynom") {
type = "polynomial"
}
else if(type == "pow"    ||
type == "power") {
type = "power"
}
else {
type = "unknown"
}

if (type=="unknown"){
stop("'type' unrecognizable.")
}

#  check rangeval if the object is not of type FEM

if (!type == "FEM") {
rangeval = as.vector(rangeval)
if (!is.numeric(rangeval)) stop("Argument rangeval is not numeric.")
if (length(rangeval) != 2) stop("Argument rangeval is not of length 2.")
if (!(rangeval[2] > rangeval[1]))
stop("Argument rangeval is not strictly increasing.")
}

#  check nbasis

if (nbasis <= 0)             stop("Argument nbasis is not positive.")
if (round(nbasis) != nbasis) stop("Argument nbasis is not an integer.")

#  check if QUADVALS is present, and set to default if not

if ((nquad == 2) && (ncol > 2)){
}
if (ncol != 2) stop("'quadvals' does not have two columns.")
}

#  check VALUES is present, and set to a single empty list if not.
if(!(length(values) == 0 || missing(values) || is.null(values))) {
n <- dim(values)[1]
k <- dim(values)[2]
stop(paste("Number of rows in 'values' not equal to number of",
if (k != nbasis)
stop(paste("Number of columns in 'values' not equal to number of",
"basis functions."))
}
else values <- vector("list",0)

#  check BASISVALUES is present, and set to vector("list",0) if not.
#  If present, it must be a two-dimensional list created by a command like
#  listobj <- matrix("list", 2, 3)

if(!(length(basisvalues) == 0 || missing(basisvalues) || !is.null(basisvalues))) {
if (!is.list(basisvalues)) stop("BASISVALUES is not a list object.")
sizevec <- dim(basisvalues)
if (length(sizevec) != 2) stop("BASISVALUES is not 2-dimensional.")
for (i in 1:sizevec[1]) {
if (length(basisvalues[[i,1]]) != dim(basisvalues[[i,2]])[1]) stop(
paste("Number of argument values not equal number",
"of values."))
}
}
else basisvalues <- vector("list",0)

#  check if DROPIND is present, and set to default if not

if(missing(dropind)) dropind <- vector("list",0)

if (length(dropind) > 0) {
#  check DROPIND
ndrop = length(dropind)
if (ndrop >= nbasis) stop('Too many index values in DROPIND.')
dropind = sort(dropind)
if (ndrop > 1 && any(diff(dropind)) == 0)
stop('Multiple index values in DROPIND.')
for (i in 1:ndrop) {
if (dropind[i] < 1 || dropind[i] > nbasis)
stop('A DROPIND index value is out of range.')
}
#  drop columns from VALUES cells if present
nvalues = length(values)
if (nvalues > 0 && length(values[[1]] > 0)) {
for (ivalue in 1:nvalues) {
derivvals = values[[ivalue]]
derivvals = derivvals[,-dropind]
values[[ivalue]] = derivvals
}
}
}

#  select the appropriate type and process

if (type=="fourier"){
paramvec   <- rangeval[2] - rangeval[1]
period     <- params[1]
if (period <= 0)  stop("Period must be positive for (a Fourier basis")
params <- period
if ((2*floor(nbasis/2)) == nbasis)  nbasis <- nbasis + 1
} else if(type=="bspline"){
if (!missing(params)){
nparams  <- length(params)
if(nparams>0){
if (params[1] <= rangeval[1])
stop("Smallest value in BREAKS not within RANGEVAL")
if (params[nparams] >= rangeval[2])
stop("Largest value in BREAKS not within RANGEVAL")
}
}
} else if(type=="expon") {
if (length(params) != nbasis)
stop("No. of parameters not equal to no. of basis fns for (exponential basisobj\$")
}   else if(type=="fdVariance") {
if (length(params) != 2)
stop("No. of parameters not equal to 8 for (FEM basisobj\$")
} else if(type=="FEM") {
# if (length(params) != 8)
#   stop("No. of parameters not equal to 8 for (FEM basisobj\$")
} else if(type=="polynomial") {
if (length(params) != nbasis)
stop("No. of parameters not equal to no. of basis fns for (polynomialal basisobj\$")
} else if(type=="power") {
if (length(params) != nbasis)
stop("No. of parameters not equal to no. of basis fns for (power basisobj\$")
} else if(type=="const") {
params <- 0
} else if(type=="monom") {
if (length(params) != nbasis)
stop("No. of parameters not equal to no. of basis fns for (monomial basisobj\$")
} else if(type=="polygonal") {
if (length(params) != nbasis)
stop("No. of parameters not equal to no. of basis fns for polygonal basisobj\$")
} else stop("Unrecognizable basis")

#  Save call

obj.call <- match.call()

#  S4 definition

# basisobj <- new("basisfd", call=obj.call, type=type, rangeval=rangeval,
#                 nbasis=nbasis,  params=params, dropind=dropind,

#  S3 definition

basisobj <- list(call=obj.call, type=type, rangeval=rangeval, nbasis=nbasis,
values=values, basisvalues=basisvalues)
oldClass(basisobj) <- "basisfd"

basisobj

}

#  --------------------------------------------------------------------------
#                  print for basisfd class
#  --------------------------------------------------------------------------

print.basisfd <- function(x, ...)
{

basisobj <- x
cat("\nBasis object:\n")
if (!inherits(basisobj, "basisfd"))
stop("Argument not a functional data object")

#  print type

cat(paste("\n  Type:  ", basisobj\$type,"\n"))

#  print range

cat(paste("\n  Range: ", basisobj\$rangeval[1],
" to ",        basisobj\$rangeval[2],"\n"))

#  return if a constant basis

if (basisobj\$type == "const") {
return()
}

#  print number of basis functions

cat(paste("\n  Number of basis functions: ",
basisobj\$nbasis,     "\n"))

#  print parameters according to type of basis

if (basisobj\$type == "fourier")
cat(paste("\n  Period: ",basisobj\$params,"\n"))
if (basisobj\$type == "bspline") {
norder <- basisobj\$nbasis - length(basisobj\$params)
cat(paste("\n  Order of spline: ", norder, "\n"))
if (length(basisobj\$params) > 0) {
print("  Interior knots")
print(basisobj\$params)
} else {
print("  There are no interior knots.")
}
}
if (basisobj\$type == "polyg") {
print("  Argument values")
print(basisobj\$params)
}
if (basisobj\$type == "expon") {
print("  Rate coefficients")
print(basisobj\$params)
}
if (basisobj\$type == "monom") {
print("  Exponents")
print(basisobj\$params)
}
if (basisobj\$type == "power") {
print("  Exponents")
print(basisobj\$params)
}

#  display indices of basis functions to be dropped

if (length(basisobj\$dropind) > 0) {
print("  Indices of basis functions to be dropped")
print(basisobj\$dropind)
}

}

#  --------------------------------------------------------------------------
#                  summary for basisfd class
#  --------------------------------------------------------------------------

summary.basisfd <- function(object, ...)
{
basisobj <- object
cat("\nBasis object:\n")
if (!inherits(basisobj, "basisfd"))
stop("Argument not a functional data object")
cat(paste("\n  Type:  ", basisobj\$type,"\n"))
cat(paste("\n  Range: ", basisobj\$rangeval[1],
" to ",        basisobj\$rangeval[2],"\n"))
if (basisobj\$type == "const") {
return()
}
cat(paste("\n  Number of basis functions: ",
basisobj\$nbasis,     "\n"))
if (basisobj\$type == "fourier")
cat(paste("\n  Period: ",basisobj\$params,"\n"))
if (length(basisobj\$dropind) > 0) {
print(paste(length(basisobj\$dropind),
"indices of basis functions to be dropped"))
}
}

#  --------------------------------------------------------------------------
#  equality for basisfd class
#  --------------------------------------------------------------------------

"==.basisfd" <- function(basis1, basis2)
{

# EQ assesses whether two bases are equivalent.

type1   <- basis1\$type
range1  <- basis1\$rangeval
nbasis1 <- basis1\$nbasis
pars1   <- basis1\$params
drop1   <- basis1\$dropind

type2   <- basis2\$type
range2  <- basis2\$rangeval
nbasis2 <- basis2\$nbasis
pars2   <- basis2\$params
drop2   <- basis2\$dropind

basisequal <- TRUE

#  check types

if (!(type1 == type2)) {
basisequal <- FALSE
return(basisequal)
}

#  check ranges

if (range1[1] != range2[1] || range1[2] != range2[2]) {
basisequal <- FALSE
return(basisequal)
}

#  check numbers of basis functions

if (nbasis1 != nbasis2) {
basisequal <- FALSE
return(basisequal)
}

#  check parameter vectors

if (!(all(pars1 == pars2))) {
basisequal <- FALSE
return(basisequal)
}

#  check indices of basis function to drop

if (!(all(drop1 == drop2))) {
basisequal <- FALSE
return(basisequal)
}

return(basisequal)

}

#  --------------------------------------------------------------------------
#  pointwise multiplication method for basisfd class
#  --------------------------------------------------------------------------

"*.basisfd" <- function (basisobj1, basisobj2)
{
# TIMES for (two basis objects sets up a basis suitable for (
#  expanding the pointwise product of two functional data
#  objects with these respective bases.
# In the absence of a true product basis system in this code,
#  the rules followed are inevitably a compromise:
#  (1) if both bases are B-splines, the norder is the sum of the
#      two orders - 1, and the breaks are the union of the
#      two knot sequences, each knot multiplicity being the maximum
#      of the multiplicities of the value in the two break sequences.
#      Order, however, is not allowed to exceed 20.
#      That is, no knot in the product knot sequence will have a
#      multiplicity greater than the multiplicities of this value
#      in the two knot sequences.
#      The rationale this rule is that order of differentiability
#      of the product at each value will be controlled  by
#      whichever knot sequence has the greater multiplicity.
#      In the case where one of the splines is order 1, or a step
#      function, the problem is dealt with by replacing the
#      original knot values by multiple values at that location
#      to give a discontinuous derivative.
#  (2) if both bases are Fourier bases, AND the periods are the
#      the same, the product is a Fourier basis with number of
#      basis functions the sum of the two numbers of basis fns.
#  (3) if only one of the bases is B-spline, the product basis
#      is B-spline with the same knot sequence and order two
#      higher.
#  (4) in all other cases, the product is a B-spline basis with
#      number of basis functions equal to the sum of the two
#      numbers of bases and equally spaced knots.

#  Of course the ranges must also match.

#  check the ranges

range1 <- basisobj1\$rangeval
range2 <- basisobj2\$rangeval
if (range1[1] != range2[1] || range1[2] != range2[2])
stop("Ranges are not equal.")

#  get the types

type1 <- basisobj1\$type
type2 <- basisobj2\$type

#  deal with constant bases

if (type1 == "const" && type2 == "const") {
prodbasisobj <- create.constant.basis(range1)
return(prodbasisobj)
}

if (type1 == "const") {
prodbasisobj <- basisobj2
return(prodbasisobj)
}

if (type2 == "const") {
prodbasisobj <- basisobj1
return(prodbasisobj)
}

#  get the numbers of basis functions

nbasis1 <- basisobj1\$nbasis
nbasis2 <- basisobj2\$nbasis

#  work through the cases

if (type1 == "bspline" && type2 == "bspline") {
#  both are bases B-splines
#  get orders
interiorknots1 <- basisobj1\$params
interiorknots2 <- basisobj2\$params
#    uniqueknots    <- sort(union(interiorknots1, interiorknots2))
interiorknots1.2 <- union(interiorknots1, interiorknots2)
uniqueknots <- {
if(is.null(interiorknots1.2)) NULL else sort(interiorknots1.2)
}
nunique <- length(uniqueknots)
multunique <- rep(0,nunique)
for (i in seq(length=nunique)) {
mult1 <- {
if(length(interiorknots1)>0)
length(interiorknots1[interiorknots1==uniqueknots[i]])
else 0
}
mult2 <- {
if(length(interiorknots2)>0)
length(interiorknots2[interiorknots2==uniqueknots[i]])
else 0
}
multunique[i] <- max(mult1,mult2)
}
#
allknots <- rep(0,sum(multunique))
m2 <- 0
for (i in seq(length=nunique)) {
m1 <- m2 + 1
m2 <- m2 + multunique[i]
allknots[m1:m2] <- uniqueknots[i]
}
norder1 <- nbasis1 - length(interiorknots1)
norder2 <- nbasis2 - length(interiorknots2)
#  norder is not allowed to exceed 20
norder  <- min(c(norder1 + norder2 - 1,20))
allbreaks  <- c(range1[1], allknots, range1[2])
nbasis <- length(allbreaks) + norder - 2
prodbasisobj <-
create.bspline.basis(range1, nbasis, norder, allbreaks)
return(prodbasisobj)
}

if (type1 == "fourier" && type2 == "fourier") {
#  both bases Fourier
#  check whether periods match
#  if they do not, default to the basis below.
period1 <- basisobj1\$params
period2 <- basisobj2\$params
nbasis  <- nbasis1 + nbasis2-1
if (period1 == period2) {
prodbasisobj <- create.fourier.basis(range1, nbasis, period1)
return(prodbasisobj)
}
}

#  Default case when all else fails: the product basis is B-spline
#  When neither basis is a B-spline basis, the order
#  is the sum of numbers of bases, but no more than 8.
#  When one of the bases if B-spline and the other isn"t,
#  the order is the smaller of 8 or the order of the spline
#  plus 2.  Under no circumstances can the order exceed 20, however.
#  See BsplineS where this restriction is tested.

if (type1 == "bspline" || type2 == "bspline") {
norder <- 8
if (type1 == "bspline") {
interiorknots1 <- basisobj1\$params
norder1        <- nbasis1 - length(interiorknots1)
norder         <- min(c(norder1+2, norder))
}
if (type2 == "bspline") {
interiorknots2 <- basisobj2\$params
norder2        <- nbasis2 - length(interiorknots2)
norder         <- min(c(norder2+2, norder))
}
} else {
#  neither basis is B-spline
norder <- min(c(8, nbasis1+nbasis2))
}
#  set up the default B-spline product basis
nbasis <- max(c(nbasis1+nbasis2, norder+1))
prodbasisobj <- create.bspline.basis(range1, nbasis, norder)
return(prodbasisobj)
}

#  ---------------------------------------------------------
#        Subscripted reference to a basis object
#  ---------------------------------------------------------

"[.basisfd" <- function(basisobj, subs=TRUE)
{
#  select subsets of basis functions in a basis object

dropind = vector("numeric", 0)
nbasis <- basisobj\$nbasis
for (i in 1:nbasis) {
if (!any(subs==i)) dropind = c(dropind, i)
}
basisobj\$dropind <- dropind
return(basisobj)
}
```

## Try the fda package in your browser

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

fda documentation built on May 29, 2024, 11:26 a.m.