Decomposition for Amplitude and Phase Variation

Description

Registration is the process of aligning peaks, valleys and other features in a sample of curves. Once the registration has taken place, this function computes two mean squared error measures, one for amplitude variation, and the other for phase variation. It also computes a squared multiple correlation index of the amount of variation in the unregistered functions is due to phase.

Usage

AmpPhaseDecomp(xfd, yfd, hfd, rng=xrng)

Arguments

Details

The decomposition can yield negative values for MS.phas if the registration does not improve the alignment of the curves, or if used to compare two registration processes based on different principles, such as is the case for functions landmarkreg and register.fd.

Value

a named list with the following components:

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

landmarkreg, register.fd, smooth.morph

Examples

#See the analysis for the growth data in the examples.

Arithmatic on functional basis objects

Description

Arithmatic on functional basis objects

Usage

## S3 method for class 'basisfd'
basis1 == basis2

Arguments

Value

basisobj1 == basisobj2 returns a logical scalar.

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2006), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

basisfd, basisfd.product arithmetic.fd

Arithmetic on functional data ('fd') objects

Description

Arithmetic on functional data objects

Usage

## S3 method for class 'fd'
e1 + e2
## S3 method for class 'fd'
e1 - e2
## S3 method for class 'fd'
e1 * e2
plus.fd(e1, e2, basisobj=NULL)
minus.fd(e1, e2, basisobj=NULL)
times.fd(e1, e2, basisobj=NULL)

Arguments

Value

A function data object corresponding to the pointwise sum, difference or product of e1 and e2.

If both arguments are functional data objects, the bases are the same, and the coefficient matrices are the same dims, the indicated operation is applied to the coefficient matrices of the two objects. In other words, e1+e2 is obtained for this case by adding the coefficient matrices from e1 and e2.

If e1 or e2 is a numeric scalar, that scalar is applied to the coefficient matrix of the functional data object.

If either e1 or e2 is a numeric vector, it must be the same length as the number of replicated functional observations in the other argument.

When both arguments are functional data objects, they need not have the same bases. However, if they don't have the same number of replicates, then one of them must have a single replicate. In the second case, the singleton function is replicated to match the number of replicates of the other function. In either case, they must have the same number of functions. When both arguments are functional data objects, and the bases are not the same, the basis used for the sum is constructed to be of higher dimension than the basis for either factor according to rules described in function TIMES for two basis objects.

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

basisfd, basisfd.product, exponentiate.fd

Examples

oldpar <- par(no.readonly=TRUE)
##
## add a parabola to itself
##
bspl4 <- create.bspline.basis(nbasis=4)
parab4.5 <- fd(c(3, -1, -1, 3)/3, bspl4)

coef2 <- matrix(c(6, -2, -2, 6)/3, 4)
dimnames(coef2) <- list(NULL, 'reps 1')

all.equal(coef(parab4.5+parab4.5), coef2)


##
## Same example with interior knots at 1/3 and 1/2
##
bspl5.3 <- create.bspline.basis(breaks=c(0, 1/3, 1))
plot(bspl5.3)
x. <- seq(0, 1, .1)
para4.5.3 <- smooth.basis(x., 4*(x.-0.5)^2, bspl5.3)$fd
plot(para4.5.3)

bspl5.2 <- create.bspline.basis(breaks=c(0, 1/2, 1))
plot(bspl5.2)
para4.5.2 <- smooth.basis(x., 4*(x.-0.5)^2, bspl5.2)$fd
plot(para4.5.2)

#str(para4.5.3+para4.5.2)

coef2. <- matrix(0, 9, 1)
dimnames(coef2.) <- list(NULL, 'rep1')

all.equal(coef(para4.5.3-para4.5.2), coef2.)


##
## product
##
quart <- para4.5.3*para4.5.2


# norder(quart) = norder(para4.5.2)+norder(para4.5.3)-1 = 7

norder(quart) == (norder(para4.5.2)+norder(para4.5.3)-1)


# para4.5.2 with knot at 0.5 and para4.5.3 with knot at 1/3
# both have (2 end points + 1 interior knot) + norder-2
#     = 5 basis functions
# quart has (2 end points + 2 interior knots)+norder-2
#     = 9 basis functions
# coefficients look strange because the knots are at
# (1/3, 1/2) and not symmetrical


all.equal(as.numeric(coef(quart)),
0.1*c(90, 50, 14, -10, 6, -2, -2, 30, 90)/9)


plot(para4.5.3*para4.5.2) # quartic, not parabolic ...

##
## product with Fourier bases
##
f3 <- fd(c(0,0,1), create.fourier.basis())
f3^2 # number of basis functions = 7?

##
## fd+numeric
##
coef1 <- matrix(c(6, 2, 2, 6)/3, 4)
dimnames(coef1) <- list(NULL, 'reps 1')

all.equal(coef(parab4.5+1), coef1)



all.equal(1+parab4.5, parab4.5+1)


##
## fd-numeric
##
coefneg <- matrix(c(-3, 1, 1, -3)/3, 4)
dimnames(coefneg) <- list(NULL, 'reps 1')

all.equal(coef(-parab4.5), coefneg)


plot(parab4.5-1)

plot(1-parab4.5)

par(oldpar)

Reshape a vector or array to have 3 dimensions.

Description

Coerce a vector or array to have 3 dimensions, preserving dimnames if feasible. Throw an error if length(dim(x)) > 3.

Usage

as.array3(x) 

Arguments

Details

1. dimx <- dim(x); ndim <- length(dimx)

2. if(ndim==3)return(x).

3. if(ndim>3)stop.

4. x2 <- as.matrix(x)

5. dim(x2) <- c(dim(x2), 1)

6. xnames <- dimnames(x)

7. if(is.list(xnames))dimnames(x2) <- list(xnames[[1]], xnames[[2]], NULL)

Value

A 3-dimensional array with names matching x

Author(s)

Spencer Graves

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

dim, dimnames checkDims3

Examples

  ##
  ## vector -> array 
  ##
  as.array3(c(a=1, b=2)) 

  ##
  ## matrix -> array 
  ##
  
  as.array3(matrix(1:6, 2))
  as.array3(matrix(1:6, 2, dimnames=list(letters[1:2],
      LETTERS[3:5]))) 

  ##
  ## array -> array 
  ##
  
  as.array3(array(1:6, 1:3)) 

  ##
  ## 4-d array 
  ##
  # These lines throw an error because the dimensionality woud be 4
  # and as.array3 only allows dimensions 3 or less.
  # if(!CRAN()) {
  #   as.array3(array(1:24, 1:4)) 
  # }

Convert a spline object to class 'fd'

Description

Translate a spline object of another class into the Functional Data (class fd) format.

Usage

as.fd(x, ...)
## S3 method for class 'fdSmooth'
as.fd(x, ...)
## S3 method for class 'function'
as.fd(x, ...)
## S3 method for class 'smooth.spline'
as.fd(x, ...)

Arguments

Details

The behavior depends on the class and nature of x.

as.fd.fdSmooth

extract the fd component

as.fd.function

Create an fd object from a function of the form created by splinefun. This will translate method = 'fmn' and 'natural' but not 'periodic': 'fmn' splines are isomorphic to standard B-splines with coincident boundary knots, which is the basis produced by create.bspline.basis. 'natural' splines occupy a subspace of this space, with the restriction that the second derivative at the end points is zero (as noted in the Wikipedia spline article). 'periodic' splines do not use coincident boundary knots and are not currently supported in fda; instead, fda uses finite Fourier bases for periodic phenomena.

as.fd.smooth.spline

Create an fd object from a smooth.spline object.

Author(s)

Spencer Graves

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

spline entry in Wikipedia https://en.wikipedia.org/wiki/Spline_(mathematics)

See Also

fd splinefun

Examples

##
## as.fd.fdSmooth
##
girlGrowthSm <- with(growth, 
  smooth.basisPar(argvals=age, y=hgtf, lambda=0.1))
girlGrowth.fd <- as.fd(girlGrowthSm)

##
## as.fd.function(splinefun(...), ...)
##
x2 <- 1:7
y2 <- sin((x2-0.5)*pi)
fd_function <- splinefun(x2, y2)
fd.  <- as.fd(fd_function)
x.   <- seq(1, 7, .02)
fdx. <- fda::eval.fd(x., fd.)

# range(y2, fx., fdx.) generates an error 2012.04.22

rfdx <- range(fdx.)

oldpar <- par(no.readonly= TRUE)
plot(range(x2), range(y2, fdx., rfdx), type='n')
points(x2, y2)
lines(x., sin((x.-0.5)*pi), lty='dashed')
lines(x., fdx., col='blue')
lines(x., eval.fd(x., fd.), col='red', lwd=3, lty='dashed')
# splinefun and as.fd(splineful(...)) are close
# but quite different from the actual function
# apart from the actual 7 points fitted,
# which are fitted exactly
# ... and there is no information in the data
# to support a better fit!

# Translate also a natural spline
fn <- splinefun(x2, y2, method='natural')
fn. <- as.fd(fn)
lines(x., fn(x.), lty='dotted', col='blue')
lines(x., eval.fd(x., fn.), col='green', lty='dotted', lwd=3)

if(!CRAN()) {
# Will NOT translate a periodic spline
# fp <- splinefun(x, y, method='periodic')
# as.fd(fp)
# Error in as.fd.function(fp) :
#  x (fp)  uses periodic B-splines, and as.fd is programmed
#   to translate only B-splines with coincident boundary knots.
}

##
## as.fd.smooth.spline ... this doesn't work (24 January 2024)
##
#cars.spl <- with(cars, smooth.spline(speed, dist))
#cars.fd  <- as.fd(cars.spl)

#plot(dist~speed, cars)
#lines(cars.spl)
#sp. <- with(cars, seq(min(speed), max(speed), len=101))
#d. <- eval.fd(sp., cars.fd)
#lines(sp., d., lty=2, col='red', lwd=3)
par(oldpar)

as.POXIXct for number of seconds since the start of 1970.

Description

as.POSIXct.numeric requires origin to be specified. This assumes that it is the start of 1970.

Usage

as.POSIXct1970(x, tz="GMT", ...)

Arguments

Details

o1970 <- strptime('1970-01-01', ' o1970. <- as.POSIXct(o1970), as.POSIXct(x, origin=o1970.)

Value

Returns a vector of class POSIXct.

Author(s)

Spencer Graves

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

as.POSIXct, ISOdate, strptime

Examples

sec <- c(0, 1, 60, 3600, 24*3600, 31*24*3600, 365*24*3600)
Sec <- as.POSIXct1970(sec)

all.equal(sec, as.numeric(Sec))

Mark Intervals on a Plot Axis

Description

Adds an axis (axisintervals) or two axes (axesIntervals) to the current plot with tick marks delimiting interval described by labels

Usage

axisIntervals(side=1, atTick1=fda::monthBegin.5, atTick2=fda::monthEnd.5,
      atLabels=fda::monthMid, labels=month.abb, cex.axis=0.9, ...)
axesIntervals(side=1:2, atTick1=fda::monthBegin.5, atTick2=fda::monthEnd.5,
      atLabels=fda::monthMid, labels=month.abb, cex.axis=0.9, las=1, ...)

Arguments

Value

The value from the third (labels) call to 'axis'. This function is usually invoked for its side effect, which is to add an axis to an already existing plot.

axesIntervals calls axisIntervals(side[1], ...) then axis(side[2], ...).

Side Effects

An axis is added to the current plot.

Author(s)

Spencer Graves

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

axis, par monthBegin.5 monthEnd.5 monthMid month.abb monthLetters

Examples

oldpar <- par(no.readonly= TRUE)
daybasis65 <- create.fourier.basis(c(0, 365), 65)

daytempfd <- with(CanadianWeather, smooth.basis(
       day.5,  dailyAv[,,"Temperature.C"], 
       daybasis65, fdnames=list("Day", "Station", "Deg C"))$fd )

with(CanadianWeather, plotfit.fd(
      dailyAv[,,"Temperature.C"], argvals=day.5,
          daytempfd, index=1, titles=place, axes=FALSE) )
# Label the horizontal axis with the month names
axisIntervals(1)
axis(2)
# Depending on the physical size of the plot,
# axis labels may not all print.
# In that case, there are 2 options:
# (1) reduce 'cex.lab'.
# (2) Use different labels as illustrated by adding
#     such an axis to the top of this plot

with(CanadianWeather, plotfit.fd(
      dailyAv[,,"Temperature.C"], argvals=day.5,
          daytempfd, index=1, titles=place, axes=FALSE) )
# Label the horizontal axis with the month names
axesIntervals()

axisIntervals(3, labels=monthLetters, cex.lab=1.2, line=-0.5)
# 'line' argument here is passed to 'axis' via '...'
par(oldpar)

Define a Functional Basis Object

Description

This is the constructor function for objects of the basisfd class. Each function that sets up an object of this class must call this function. This includes functions create.bspline.basis, create.constant.basis, create.fourier.basis, and so forth that set up basis objects of a specific type. Ordinarily, user of the functional data analysis software will not need to call this function directly, but these notes are valuable to understanding what the "slots" or "members" of the basisfd class are.

Usage

basisfd(type, rangeval, nbasis, params,
        dropind=vector('list', 0),
        quadvals=vector('list', 0),
        values=vector("list", 0),
        basisvalues=vector('list', 0))

Arguments

Details

Previous versions of the 'fda' software used the name basis for this class, and the code in Matlab still does. However, this class name was already used elsewhere in the S languages, and there was a potential for a clash that might produce mysterious and perhaps disastrous consequences.

To check that an object is of this class, use function is.basis.

It is comparatively simple to add new basis types. The code in the following functions needs to be estended to allow for the new type: basisfd, getbasismatrix and getbasispenalty. In addition, a new "create" function should be written for the new type, as well as functions analogous to fourier and fourierpen for evaluating basis functions for basis penalty matrices.

The "create" function names are rather long, and users who mind all that typing might be advised to modify these to versions with shorter names, such as "splbas", "conbas", and etc. However, a principle of good programming practice is to keep the code readable, preferably by somebody other than the programmer.

Normally only developers of new basis types will actually need to use this function, so no examples are provided.

Value

an object of class basisfd, being a list with the following components:

References

Ramsay, James O., Hooker, Giles and Graves, Spencer (2009) Functional Data Analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

is.basis, is.eqbasis, plot.basisfd, getbasismatrix, getbasispenalty, create.bspline.basis, create.constant.basis, create.exponential.basis, create.fourier.basis, create.monomial.basis, create.polygonal.basis, create.power.basis

Product of two basisfd objects

Description

pointwise multiplication method for basisfd class

Usage

## S3 method for class 'basisfd'
basisobj1 * basisobj2

Arguments

Details

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. 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 eachy 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.

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

basisfd

Examples

  f1 <- create.fourier.basis()
  f1.2 <- f1*f1
  
  all.equal(f1.2, create.fourier.basis(nbasis=5))
  

Create a bivariate functional data object

Description

This function creates a bivariate functional data object, which consists of two bases for expanding a functional data object of two variables, s and t, and a set of coefficients defining this expansion. The bases are contained in "basisfd" objects.

Usage

bifd (coef=matrix(0,2,1), sbasisobj=create.bspline.basis(),
      tbasisobj=create.bspline.basis(), fdnames=defaultnames)

Arguments

Value

A bivariate functional data object = a list of class 'bifd' with the following components:

Author(s)

Spencer Graves

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

basisfd objAndNames

Examples

Bspl2 <- create.bspline.basis(nbasis=2, norder=1)
Bspl3 <- create.bspline.basis(nbasis=3, norder=2)

(bBspl2.3 <- bifd(array(1:6, dim=2:3), Bspl2, Bspl3))
str(bBspl2.3)

Define a Bivariate Functional Parameter Object

Description

Functional parameter objects are used as arguments to functions that estimate functional parameters, such as smoothing functions like smooth.basis. A bivariate functional parameter object supplies the analogous information required for smoothing bivariate data using a bivariate functional data object $x(s,t)$. The arguments are the same as those for fdPar objects, except that two linear differential operator objects and two smoothing parameters must be applied, each pair corresponding to one of the arguments $s$ and $t$ of the bivariate functional data object.

Usage

bifdPar(bifdobj, Lfdobjs=int2Lfd(2), Lfdobjt=int2Lfd(2), lambdas=0, lambdat=0,
      estimate=TRUE)

Arguments

Value

a bivariate functional parameter object (i.e., an object of class bifdPar), which is a list with the following components:

Source

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009) Functional Data Analysis in R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

linmod

Examples

#See the prediction of precipitation using temperature as
#the independent variable in the analysis of the daily weather
#data, and the analysis of the Swedish mortality data.

B-Spline Penalty Matrix

Description

Computes the matrix defining the roughness penalty for functions expressed in terms of a B-spline basis.

Usage

  bsplinepen(basisobj, Lfdobj=2, rng=basisobj$rangeval, returnMatrix=FALSE)

Arguments

Details

A roughness penalty for a function $x(t)$ is defined by integrating the square of either the derivative of $x(t)$ or, more generally, the result of applying a linear differential operator $L$ to it. The most common roughness penalty is the integral of the square of the second derivative, and this is the default. To apply this roughness penalty, the matrix of inner products of the basis functions (possibly after applying the linear differential operator to them) defining this function is necessary. This function just calls the roughness penalty evaluation function specific to the basis involved.

Value

a symmetric matrix of order equal to the number of basis functions defined by the B-spline basis object. Each element is the inner product of two B-spline basis functions after applying the derivative or linear differential operator defined by Lfdobj.

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

Examples

##
## bsplinepen with only one basis function
##
bspl1.1 <- create.bspline.basis(nbasis=1, norder=1)
pen1.1 <- bsplinepen(bspl1.1, 0)

##
## bspline pen for a cubic spline with knots at seq(0, 1, .1)
##
basisobj <- create.bspline.basis(c(0,1),13)
#  compute the 13 by 13 matrix of inner products of second derivatives
penmat <- bsplinepen(basisobj)

##
## with rng of class Date or POSIXct
##
# Date
invasion1 <- as.Date('1775-09-04')
invasion2 <- as.Date('1812-07-12')
earlyUS.Canada <- c(invasion1, invasion2)
BspInvade1 <- create.bspline.basis(earlyUS.Canada)
Binvadmat <- bsplinepen(BspInvade1)

# POSIXct
AmRev.ct <- as.POSIXct1970(c('1776-07-04', '1789-04-30'))
BspRev1.ct <- create.bspline.basis(AmRev.ct)
Brevmat <- bsplinepen(BspRev1.ct)

B-spline Basis Function Values

Description

Evaluates a set of B-spline basis functions, or a derivative of these functions, at a set of arguments.

Usage

bsplineS(x, breaks, norder=4, nderiv=0, returnMatrix=FALSE)

Arguments

Value

a matrix of function values. The number of rows equals the number of arguments, and the number of columns equals the number of basis functions.

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

Examples

# Minimal example:  A B-spline of order 1 (i.e., a step function)
# with 0 interior knots:
bS <- bsplineS(seq(0, 1, .2), 0:1, 1, 0)

# check

all.equal(bS, matrix(1, 6))


#  set up break values at 0.0, 0.2,..., 0.8, 1.0.
breaks <- seq(0,1,0.2)
#  set up a set of 11 argument values
x <- seq(0,1,0.1)
#  the order willl be 4, and the number of basis functions
#  is equal to the number of interior break values (4 here)
#  plus the order, for a total here of 8.
norder <- 4
#  compute the 11 by 8 matrix of basis function values
basismat <- bsplineS(x, breaks, norder)

Canadian average annual weather cycle

Description

Daily temperature and precipitation at 35 different locations in Canada averaged over 1960 to 1994.

Usage

CanadianWeather
daily

Format

'CanadianWeather' and 'daily' are lists containing essentially the same data. 'CanadianWeather' may be preferred for most purposes; 'daily' is included primarily for compatibility with scripts written before the other format became available and for compatibility with the Matlab 'fda' code.

CanadianWeather

A list with the following components:

dailyAv

a three dimensional array c(365, 35, 3) summarizing data collected at 35 different weather stations in Canada on the following:

[,,1] = [,, 'Temperature.C']: average daily temperature for each day of the year

[,,2] = [,, 'Precipitation.mm']: average daily rainfall for each day of the year rounded to 0.1 mm.

[,,3] = [,, 'log10precip']: base 10 logarithm of Precipitation.mm after first replacing 27 zeros by 0.05 mm (Ramsay and Silverman 2006, p. 248).

place

Names of the 35 different weather stations in Canada whose data are summarized in 'dailyAv'. These names vary between 6 and 11 characters in length. By contrast, daily[["place"]] which are all 11 characters, with names having fewer characters being extended with trailing blanks.

province

names of the Canadian province containing each place

coordinates

a numeric matrix giving 'N.latitude' and 'W.longitude' for each place.

region

Which of 4 climate zones contain each place: Atlantic, Pacific, Continental, Arctic.

monthlyTemp

A matrix of dimensions (12, 35) giving the average temperature in degrees celcius for each month of the year.

monthlyPrecip

A matrix of dimensions (12, 35) giving the average daily precipitation in millimeters for each month of the year.

geogindex

Order the weather stations from East to West to North

daily

A list with the following components:

place

Names of the 35 different weather stations in Canada whose data are summarized in 'dailyAv'. These names are all 11 characters, with shorter names being extended with trailing blanks. This is different from CanadianWeather[["place"]], where trailing blanks have been dropped.

tempav

a matrix of dimensions (365, 35) giving the average temperature in degrees celcius for each day of the year. This is essentially the same as CanadianWeather[["dailyAv"]][,, "Temperature.C"].

precipav

a matrix of dimensions (365, 35) giving the average temperature in degrees celcius for each day of the year. This is essentially the same as CanadianWeather[["dailyAv"]][,, "Precipitation.mm"].

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

monthAccessories MontrealTemp

Examples

##
## 1.  Plot (latitude & longitude) of stations by region
##
with(CanadianWeather, plot(-coordinates[, 2], coordinates[, 1], type='n',
                           xlab="West Longitude", ylab="North Latitude",
                           axes=FALSE) )
Wlon <- pretty(CanadianWeather$coordinates[, 2])
axis(1, -Wlon, Wlon)
axis(2)

rgns <- 1:4
names(rgns) <- c('Arctic', 'Atlantic', 'Continental', 'Pacific')
Rgns <- rgns[CanadianWeather$region]
with(CanadianWeather, points(-coordinates[, 2], coordinates[, 1],
                             col=Rgns, pch=Rgns) )
legend('topright', legend=names(rgns), col=rgns, pch=rgns)

##
## 2.  Plot dailyAv[, 'Temperature.C'] for 4 stations
##
data(CanadianWeather)
# Expand the left margin to allow space for place names
op <- par(mar=c(5, 4, 4, 5)+.1)
# Plot
stations <- c("Pr. Rupert", "Montreal", "Edmonton", "Resolute")
matplot(day.5, CanadianWeather$dailyAv[, stations, "Temperature.C"],
        type="l", axes=FALSE, xlab="", ylab="Mean Temperature (deg C)")
axis(2, las=1)
# Label the horizontal axis with the month names
axis(1, monthBegin.5, labels=FALSE)
axis(1, monthEnd.5, labels=FALSE)
axis(1, monthMid, monthLetters, tick=FALSE)
# Add the monthly averages
matpoints(monthMid, CanadianWeather$monthlyTemp[, stations])
# Add the names of the weather stations
mtext(stations, side=4,
      at=CanadianWeather$dailyAv[365, stations, "Temperature.C"],
     las=1)
# clean up
par(op)

Functional Canonical Correlation Analysis

Description

Carry out a functional canonical correlation analysis with regularization or roughness penalties on the estimated canonical variables.

Usage

cca.fd(fdobj1, fdobj2=fdobj1, ncan = 2,
       ccafdPar1=fdPar(basisobj1, 2, 1e-10),
       ccafdPar2=ccafdPar1, centerfns=TRUE)

Arguments

Value

an object of class cca.fd with the 5 slots:

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

plot.cca.fd, varmx.cca.fd, pca.fd

Examples

#  Canonical correlation analysis of knee-hip curves

gaittime  <- (1:20)/21
gaitrange <- c(0,1)
gaitbasis <- create.fourier.basis(gaitrange,21)
lambda    <- 10^(-11.5)
harmaccelLfd <- vec2Lfd(c(0, 0, (2*pi)^2, 0))

gaitfd    <- fda::fd(matrix(0,gaitbasis$nbasis,1), gaitbasis)
gaitfdPar <- fda::fdPar(gaitfd, harmaccelLfd, lambda)
gaitfd    <- fda::smooth.basis(gaittime, gait, gaitfdPar)$fd
ccafdPar  <- fda::fdPar(gaitfd, harmaccelLfd, 1e-8)
ccafd0    <- cca.fd(gaitfd[,1], gaitfd[,2], ncan=3, ccafdPar, ccafdPar)
#  display the canonical correlations
round(ccafd0$ccacorr[1:6],3)
#  compute a VARIMAX rotation of the canonical variables
ccafd <- varmx.cca.fd(ccafd0)
#  plot the canonical weight functions
oldpar <- par(no.readonly= TRUE)
plot.cca.fd(ccafd)
par(oldpar)

Center Functional Data

Description

Subtract the pointwise mean from each of the functions in a functional data object; that is, to center them on the mean function.

Usage

 center.fd(fdobj)

Arguments

Value

a functional data object whose mean is zero.

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

mean.fd, sum.fd, stddev.fd, std.fd

Examples

daytime    <- (1:365)-0.5
daybasis   <- create.fourier.basis(c(0,365), 365)
harmLcoef  <- c(0,(2*pi/365)^2,0)
harmLfd    <- vec2Lfd(harmLcoef, c(0,365))
templambda <- 0.01
dayfd      <- fda::fd(matrix(0, daybasis$nbasis, 1), daybasis)
tempfdPar  <- fda::fdPar(dayfd, harmLfd, templambda)

# do not run on CRAN because it takes too long.
tempfd     <- smooth.basis(daytime,
       CanadianWeather$dailyAv[,,"Temperature.C"], tempfdPar)$fd
tempctrfd  <- center.fd(tempfd)
oldpar <- par(no.readonly= TRUE)
plot(tempctrfd, xlab="Day", ylab="deg. C",
     main = "Centered temperature curves")
par(oldpar)

Compare dimensions and dimnames of arrays

Description

Compare selected dimensions and dimnames of arrays, coercing objects to 3-dimensional arrays and either give an error or force matching.

Usage

checkDim3(x, y=NULL, xdim=1, ydim=1, defaultNames='x',
         subset=c('xiny', 'yinx', 'neither'),
         xName=substring(deparse(substitute(x)), 1, 33),
         yName=substring(deparse(substitute(y)), 1, 33) )
checkDims3(x, y=NULL, xdim=2:3, ydim=2:3, defaultNames='x',
         subset=c('xiny', 'yinx', 'neither'),
         xName=substring(deparse(substitute(x)), 1, 33),
         yName=substring(deparse(substitute(y)), 1, 33) )

Arguments

Details

For checkDims3, confirm that xdim and ydim have the same length, and call checkDim3 for each pair.

For checkDim3, proceed as follows:

1. if((xdim>3) | (ydim>3)) throw an error.

2. ixperm <- list(1:3, c(2, 1, 3), c(3, 2, 1))[xdim]; iyperm <- list(1:3, c(2, 1, 3), c(3, 2, 1))[ydim];

3. x3 <- aperm(as.array3(x), ixperm); y3 <- aperm(as.array3(y), iyperm)

4. xNames <- dimnames(x3); yNames <- dimnames(y3)

5. Check subset. For example, for subset='xiny', use the following: if(is.null(xNames)){ if(dim(x3)[1]>dim(y3)[1]) stop else y. <- y3[1:dim(x3)[1],,] dimnames(x) <- list(yNames[[1]], NULL, NULL) } else { if(is.null(xNames[[1]])){ if(dim(x3)[1]>dim(y3)[1]) stop else y. <- y3[1:dim(x3)[1],,] dimnames(x3)[[1]] <- yNames[[1]] } else { if(any(!is.element(xNames[[1]], yNames[[1]])))stop else y. <- y3[xNames[[1]],,] } }

6. return(list(x=aperm(x3, ixperm), y=aperm(y., iyperm)))

Value

a list with components x and y.

Author(s)

Spencer Graves

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

as.array3 plotfit.fd

Examples

# Select the first two rows of y
stopifnot(all.equal(
checkDim3(1:2, 3:5),
list(x=array(1:2, c(2,1,1), list(c('x1','x2'), NULL, NULL)),
     y=array(3:4, c(2,1,1), list(c('x1','x2'), NULL, NULL)) )
))

# Select the first two rows of a matrix y
stopifnot(all.equal(
checkDim3(1:2, matrix(3:8, 3)),
list(x=array(1:2,         c(2,1,1), list(c('x1','x2'), NULL, NULL)),
     y=array(c(3:4, 6:7), c(2,2,1), list(c('x1','x2'), NULL, NULL)) )
))

# Select the first column of y
stopifnot(all.equal(
checkDim3(1:2, matrix(3:8, 3), 2, 2),
list(x=array(1:2,         c(2,1,1), list(NULL, 'x', NULL)),
     y=array(3:5, c(3,1,1), list(NULL, 'x', NULL)) )
))

# Select the first two rows and the first column of y
stopifnot(all.equal(
checkDims3(1:2, matrix(3:8, 3), 1:2, 1:2),
list(x=array(1:2, c(2,1,1), list(c('x1','x2'), 'x', NULL)),
     y=array(3:4, c(2,1,1), list(c('x1','x2'), 'x', NULL)) )
))

# Select the first 2 rows of y
x1 <- matrix(1:4, 2, dimnames=list(NULL, LETTERS[2:3]))
x1a <- x1. <- as.array3(x1)
dimnames(x1a)[[1]] <- c('x1', 'x2')
y1 <- matrix(11:19, 3, dimnames=list(NULL, LETTERS[1:3]))
y1a <- y1. <- as.array3(y1)
dimnames(y1a)[[1]] <- c('x1', 'x2', 'x3')

stopifnot(all.equal(
checkDim3(x1, y1),
list(x=x1a, y=y1a[1:2, , , drop=FALSE])
))

# Select columns 2 & 3 of y
stopifnot(all.equal(
checkDim3(x1, y1, 2, 2),
list(x=x1., y=y1.[, 2:3, , drop=FALSE ])
))

# Select the first 2 rows and  columns 2 & 3 of y
stopifnot(all.equal(
checkDims3(x1, y1, 1:2, 1:2),
list(x=x1a, y=y1a[1:2, 2:3, , drop=FALSE ])
))

# y = columns 2 and 3 of x
x23 <- matrix(1:6, 2, dimnames=list(letters[2:3], letters[1:3]))
x23. <- as.array3(x23)
stopifnot(all.equal(
checkDim3(x23, xdim=1, ydim=2),
list(x=x23., y=x23.[, 2:3,, drop=FALSE ])
))

# Transfer dimnames from y to x
x4a <- x4 <- matrix(1:4, 2)
y4 <- matrix(5:8, 2, dimnames=list(letters[1:2], letters[3:4]))
dimnames(x4a) <- dimnames(t(y4))
stopifnot(all.equal(
checkDims3(x4, y4, 1:2, 2:1),
list(x=as.array3(x4a), y=as.array3(y4))
))

# as used in plotfit.fd
daybasis65 <- create.fourier.basis(c(0, 365), 65)

daytempfd <- with(CanadianWeather, smooth.basis(
       day.5, dailyAv[,,"Temperature.C"], 
       daybasis65, fdnames=list("Day", "Station", "Deg C"))$fd )

defaultNms <- with(daytempfd, c(fdnames[2], fdnames[3], x='x'))
subset <- checkDims3(CanadianWeather$dailyAv[, , "Temperature.C"],
               daytempfd$coef, defaultNames=defaultNms)
# Problem:  dimnames(...)[[3]] = '1'
# Fix:
subset3 <- checkDims3(
        CanadianWeather$dailyAv[, , "Temperature.C", drop=FALSE],
               daytempfd$coef, defaultNames=defaultNms)

Does an argument satisfy required conditions?

Description

Check whether an argument is a logical vector of a certain length or a numeric vector in a certain range and issue an appropriate error or warning if not:

checkLogical throws an error or returns FALSE with a warning unless x is a logical vector of exactly the required length.

checkNumeric throws an error or returns FALSE with a warning unless x is either NULL or a numeric vector of at most length with x in the desired range.

checkLogicalInteger returns a logical vector of exactly length unless x is neither NULL nor logical of the required length nor numeric with x in the desired range.

Usage

checkLogical(x, length., warnOnly=FALSE)
checkNumeric(x, lower, upper, length., integer=TRUE, unique=TRUE,
             inclusion=c(TRUE,TRUE), warnOnly=FALSE)
checkLogicalInteger(x, length., warnOnly=FALSE)

Arguments

Details

1. xName <- deparse(substitute(x)) to use in any required error or warning.

2. if(is.null(x)) handle appropriately: Return FALSE for checkLogical, TRUE for checkNumeric and rep(TRUE, length.) for checkLogicalInteger.

3. Check class(x).

4. Check other conditions.

Value

checkLogical returns a logical vector of the required length., unless it issues an error message.

checkNumeric returns a numeric vector of at most length. with all elements between lower and upper, and optionally unique, unless it issues an error message.

checkLogicalInteger returns a logical vector of the required length., unless it issues an error message.

Author(s)

Spencer Graves

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

Examples

##
## checkLogical
##
checkLogical(NULL, length=3, warnOnly=TRUE)
checkLogical(c(FALSE, TRUE, TRUE), length=4, warnOnly=TRUE)
checkLogical(c(FALSE, TRUE, TRUE), length=3)

##
## checkNumeric
##
checkNumeric(NULL, lower=1, upper=3)
checkNumeric(1:3, 1, 3)
checkNumeric(1:3, 1, 3, inclusion=FALSE, warnOnly=TRUE)
checkNumeric(pi, 1, 4, integer=TRUE, warnOnly=TRUE)
checkNumeric(c(1, 1), 1, 4, warnOnly=TRUE)
checkNumeric(c(1, 1), 1, 4, unique=FALSE, warnOnly=TRUE)

##
## checkLogicalInteger
##
checkLogicalInteger(NULL, 3)
checkLogicalInteger(c(FALSE, TRUE), warnOnly=TRUE) 
checkLogicalInteger(1:2, 3) 
checkLogicalInteger(2, warnOnly=TRUE) 
checkLogicalInteger(c(2, 4), 3, warnOnly=TRUE)

##
## checkLogicalInteger names its calling function 
## rather than itself as the location of error detection
## if possible
##
tstFun <- function(x, length., warnOnly=FALSE){
   checkLogicalInteger(x, length., warnOnly) 
}
tstFun(NULL, 3)
tstFun(4, 3, warnOnly=TRUE)

tstFun2 <- function(x, length., warnOnly=FALSE){
   tstFun(x, length., warnOnly)
}
tstFun2(4, 3, warnOnly=TRUE)

Extract functional coefficients

Description

Obtain the coefficients component from a functional object (functional data, class fd, functional parameter, class fdPar, a functional smooth, class fdSmooth, or a Taylor spline representation, class Taylor.

Usage

## S3 method for class 'fd'
coef(object, ...)
## S3 method for class 'fdPar'
coef(object, ...)
## S3 method for class 'fdSmooth'
coef(object, ...)
## S3 method for class 'fd'
coefficients(object, ...)
## S3 method for class 'fdPar'
coefficients(object, ...)
## S3 method for class 'fdSmooth'
coefficients(object, ...)

Arguments

Details

Functional representations are evaluated by multiplying a basis function matrix times a coefficient vector, matrix or 3-dimensional array. (The basis function matrix contains the basis functions as columns evaluated at the evalarg values as rows.)

Value

A numeric vector or array of the coefficients.

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

coef fd fdPar smooth.basisPar smooth.basis

Examples

##
## coef.fd
##
bspl1.1 <- create.bspline.basis(norder=1, breaks=0:1)
fd.bspl1.1 <- fd(0, basisobj=bspl1.1)
coef(fd.bspl1.1)


##
## coef.basisPar 
##
rangeval <- c(-3,3)
#  set up some standard normal data
x <- rnorm(50)
#  make sure values within the range
x[x < -3] <- -2.99
x[x >  3] <-  2.99
#  set up basis for W(x)
basisobj <- create.bspline.basis(rangeval, 11)
#  set up initial value for Wfdobj
Wfd0 <- fd(matrix(0,11,1), basisobj)
WfdParobj <- fdPar(Wfd0)

coef(WfdParobj)


##
## coef.fdSmooth
##

girlGrowthSm <- with(growth, smooth.basisPar(argvals=age, y=hgtf, 
                                             lambda=0.1)$fd)
coef(girlGrowthSm)

Correlation matrix from functional data object(s)

Description

Compute a correlation matrix for one or two functional data objects.

Usage

cor.fd(evalarg1, fdobj1, evalarg2=evalarg1, fdobj2=fdobj1)

Arguments

Details

1. var1 <- var.fd(fdobj1) 2. evalVar1 <- eval.bifd(evalarg1, evalarg1, var1) 3. if(missing(fdobj2)) Convert evalVar1 to correlations 4. else: 4.1. var2 <- var.fd(fdobj2) 4.2. evalVar2 <- eval.bifd(evalarg2, evalarg2, var2) 4.3. var12 <- var.df(fdobj1, fdobj2) 4.4. evalVar12 <- eval.bifd(evalarg1, evalarg2, var12) 4.5. Convert evalVar12 to correlations

Value

A matrix or array:

With one or two functional data objects, fdobj1 and possibly fdobj2, the value is a matrix of dimensions length(evalarg1) by length(evalarg2) giving the correlations at those points of fdobj1 if missing(fdobj2) or of correlations between eval.fd(evalarg1, fdobj1) and eval.fd(evalarg2, fdobj2).

With a single multivariate data object with k variables, the value is a 4-dimensional array of dim = c(nPts, nPts, 1, choose(k+1, 2)), where nPts = length(evalarg1).

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

mean.fd, sd.fd, std.fd stdev.fd var.fd

Examples

daybasis3 <- create.fourier.basis(c(0, 365))
daybasis5 <- create.fourier.basis(c(0, 365), 5)
tempfd3 <- with(CanadianWeather, smooth.basis(
       day.5, dailyAv[,,"Temperature.C"], 
       daybasis3, fdnames=list("Day", "Station", "Deg C"))$fd )
precfd5 <- with(CanadianWeather, smooth.basis(
       day.5, dailyAv[,,"log10precip"], 
       daybasis5, fdnames=list("Day", "Station", "Deg C"))$fd )

# Correlation matrix for a single functional data object
(tempCor3 <- cor.fd(seq(0, 356, length=4), tempfd3))

# Cross correlation matrix between two functional data objects 
# Compare with structure described above under 'value':
(tempPrecCor3.5 <- cor.fd(seq(0, 365, length=4), tempfd3,
                          seq(0, 356, length=6), precfd5))

# The following produces contour and perspective plots

daybasis65 <- create.fourier.basis(rangeval=c(0, 365), nbasis=65)
daytempfd <- with(CanadianWeather, smooth.basis(
       day.5, dailyAv[,,"Temperature.C"], 
       daybasis65, fdnames=list("Day", "Station", "Deg C"))$fd )
dayprecfd <- with(CanadianWeather, smooth.basis(
       day.5, dailyAv[,,"log10precip"], 
       daybasis65, fdnames=list("Day", "Station", "log10(mm)"))$fd )

str(tempPrecCor <- cor.fd(weeks, daytempfd, weeks, dayprecfd))
# dim(tempPrecCor)= c(53, 53)

op <- par(mfrow=c(1,2), pty="s")
contour(weeks, weeks, tempPrecCor, 
        xlab="Average Daily Temperature",
        ylab="Average Daily log10(precipitation)",
        main=paste("Correlation function across locations\n",
          "for Canadian Anual Temperature Cycle"),
        cex.main=0.8, axes=FALSE)
axisIntervals(1, atTick1=seq(0, 365, length=5), atTick2=NA, 
            atLabels=seq(1/8, 1, 1/4)*365,
            labels=paste("Q", 1:4) )
axisIntervals(2, atTick1=seq(0, 365, length=5), atTick2=NA, 
            atLabels=seq(1/8, 1, 1/4)*365,
            labels=paste("Q", 1:4) )
persp(weeks, weeks, tempPrecCor,
      xlab="Days", ylab="Days", zlab="Correlation")
mtext("Temperature-Precipitation Correlations", line=-4, outer=TRUE)
par(op)

# Correlations and cross correlations
# in a bivariate functional data object
gaittime   <- (1:20)/21
gaitbasis5 <- create.fourier.basis(c(0,1),nbasis=5)
gaitfd5    <- smooth.basis(gaittime, gait, gaitbasis5)$fd

gait.t3 <- (0:2)/2
(gaitCor3.5 <- cor.fd(gait.t3, gaitfd5))
# Check the answers with manual computations
gait3.5 <- eval.fd(gait.t3, gaitfd5)
all.equal(cor(t(gait3.5[,,1])), gaitCor3.5[,,,1])
# TRUE
all.equal(cor(t(gait3.5[,,2])), gaitCor3.5[,,,3])
# TRUE
all.equal(cor(t(gait3.5[,,2]), t(gait3.5[,,1])),
               gaitCor3.5[,,,2])
# TRUE

# NOTE:
dimnames(gaitCor3.5)[[4]]
# [1] Hip-Hip
# [2] Knee-Hip 
# [3] Knee-Knee
# If [2] were "Hip-Knee", then
# gaitCor3.5[,,,2] would match 
# cor(t(gait3.5[,,1]), t(gait3.5[,,2]))
# *** It does NOT.  Instead, it matches:  
# cor(t(gait3.5[,,2]), t(gait3.5[,,1]))

Estimate of the covariance surface

Description

Function covPACE does a bivariate smoothing for estimating the covariance surface for data that has not yet been smoothed

Usage

  covPACE(data,rng , time, meanfd, basis, lambda, Lfdobj)

Arguments

Value

a list with these two named entries:

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

Test if running as CRAN

Description

This function allows package developers to run tests themselves that should not run on CRAN or with "R CMD check –as-cran" because of compute time constraints with CRAN tests.

Usage

  CRAN(CRAN_pattern, n_R_CHECK4CRAN)

Arguments

Details

The "Writing R Extensions" manual says that "R CMD check" can be customized "by setting environment variables _R_CHECK_*_:, as described in" the Tools section of the "R Internals" manual.

'R CMD check' was tested with R 3.0.1 under Fedora 18 Linux and with Rtools 3.0 from April 16, 2013 under Windows 7. With the '–as-cran' option, 7 matches were found; without it, only 3 were found. These numbers were unaffected by the presence or absence of the '–timings' parameter. On this basis, the default value of n_R_CHECK4CRAN was set at 5.

1. x. <- Sys.getenv()

2. Fix CRAN_pattern and n_R_CHECK4CRAN if missing.

3. Let i be the indices of x. whose names match all the patterns in the vector x.

4. Assume this is CRAN if length(i) >= n_R_CHECK4CRAN

Value

a logical scalar with attributes 'Sys.getenv' containing the results of Sys.getenv() and 'matches' contining i per step 3 above.

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

Sys.getenv

Examples

  # Look in fda-ex.R
  # for the results from R CMD check
  # Modified defaults were tested with environment variables
  # _CRAN_pattern_ = 'A' and _n_R_CHECK4CRAN_ = '10'
  # 'R CMD check' found 26 matches and CRAN() returned TRUE
  cran <- CRAN()
  str(cran)
  gete <- attr(cran, 'Sys.getenv')
  (ngete <- names(gete))
  iget <- grep('^_', names(gete))
  gete[iget]
  # if (!CRAN()) {
  #   if (CRAN()) {
  #     stop("CRAN")
  #   } else {
  #     stop("NOT CRAN")
  #   }
  # }

Create Basis Set for Functional Data Analysis

Description

Functional data analysis proceeds by selecting a finite basis set and fitting data to it. The current fda package supports fitting via least squares penalized with lambda times the integral over the (finite) support of the basis set of the squared deviations from a linear differential operator.

Details

The most commonly used basis in fda is probably B-splines. For periodic phenomena, Fourier bases are quite useful. A constant basis is provided to facilitation arithmetic with functional data objects. To restrict attention to solutions of certain differential equations, it may be useful to use a corresponding basis set such as exponential, monomial or power basis sets.

Power bases support the use of negative and fractional powers, while monomial bases are restricted only to nonnegative integer exponents.

The polygonal basis is essentially a B-spline of order 2, degree 1.

The following summarizes arguments used by some or all of the current create.basis functions:

rangeval

a vector of length 2 giving the lower and upper limits of the range of permissible values for the function argument.

For bspline bases, this can be inferred from range(breaks). For polygonal bases, this can be inferred from range(argvals). In all other cases, this defaults to 0:1.

nbasis

an integer giving the number of basis functions.

This is not used for two of the create.basis functions: For constant this is 1, so there is no need to specify it. For polygonal bases, it is length(argvals), and again there is no need to specify it.

For bspline bases, if nbasis is not specified, it defaults to (length(breaks) + norder - 2) if breaks is provided. Otherwise, nbasis defaults to 20 for bspline bases.

For exponential bases, if nbasis is not specified, it defaults to length(ratevec) if ratevec is provided. Otherwise, in fda_2.0.2, ratevec defaults to 1, which makes nbasis = 1; in fda_2.0.4, ratevec will default to 0:1, so nbasis will then default to 2.

For monomial and power bases, if nbasis is not specified, it defaults to length(exponents) if exponents is provided. Otherwise, nbasis defaults to 2 for monomial and power bases. (Temporary exception: In fda_2.0.2, the default nbasis for power bases is 1. This will be increased to 2 in fda_2.0.4.)

In addition to rangeval and nbasis, all but constant bases have one or two parameters unique to that basis type or shared with one other:

bspline

Argument norder = the order of the spline, which is one more than the degree of the polynomials used. This defaults to 4, which gives cubic splines.

Argument breaks = the locations of the break or join points; also called knots. This defaults to seq(rangeval[1], rangeval[2], nbasis-norder+2).

polygonal

Argument argvals = the locations of the break or join points; also called knots. This defaults to seq(rangeval[1], rangeval[2], nbasis).

fourier

Argument period defaults to diff(rangeval).

exponential

Argument ratevec. In fda_2.0.2, this defaulted to 1. In fda_2.0.3, it will default to 0:1.

monomial, power

Argument exponents. Default = 0:(nbasis-1). For monomial bases, exponents must be distinct nonnegative integers. For power bases, they must be distinct real numbers.

Beginning with fda_2.1.0, the last 6 arguments for all the create.basis functions will be as follows; some but not all are available in the previous versions of fda:

dropind

a vector of integers specifiying the basis functions to be dropped, if any.

quadvals

a matrix with two columns and a number of rows equal to the number of quadrature points for numerical evaluation of the penalty integral. The first column of quadvals contains the quadrature points, and the second column the quadrature weights. A minimum of 5 values are required for each inter-knot interval, and that is often enough. For Simpson's rule, these points are equally spaced, and the weights are proportional to 1, 4, 2, 4, ..., 2, 4, 1.

values

a list of matrices with one row for each row of quadvals and one column for each basis function. The elements of the list correspond to the basis functions and their derivatives evaluated at the quadrature points contained in the first column of quadvals.

basisvalues

A list 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 the following:

basisobj$basisvalues <- vector("list",1)

basisobj$basisvalues[[1]] <- list(args=evalargs, values=basismat)

names

either a character vector of the same length as the number of basis functions or a simple stem used to construct such a vector.

For bspline bases, this defaults to paste('bspl', norder, '.', 1:nbreaks, sep=”).

For other bases, there are crudely similar defaults.

axes

an optional list used by selected plot functions to create custom axes. If this axes argument is not NULL, functions plot.basisfd, plot.fd, plot.fdSmooth plotfit.fd, plotfit.fdSmooth, and plot.Lfd will create axes via do.call(x$axes[[1]], x$axes[-1]). The primary example of this is to create CanadianWeather plots using list("axesIntervals")

Author(s)

J. O. Ramsay and Spencer Graves

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

create.bspline.basis create.constant.basis create.exponential.basis create.fourier.basis create.monomial.basis create.polygonal.basis create.power.basis

Create a B-spline Basis

Description

Functional data objects are constructed by specifying a set of basis functions and a set of coefficients defining a linear combination of these basis functions. The B-spline basis is used for non-periodic functions. B-spline basis functions are polynomial segments jointed end-to-end at at argument values called knots, breaks or join points. The segments have specifiable smoothness across these breaks. B-spline basis functions have the advantages of very fast computation and great flexibility. A polygonal basis generated by create.polygonal.basis is essentially a B-spline basis of order 2, degree 1. Monomial and polynomial bases can be obtained as linear transformations of certain B-spline bases.

Usage

create.bspline.basis(rangeval=NULL, nbasis=NULL, norder=4,
      breaks=NULL, dropind=NULL, quadvals=NULL, values=NULL,
      basisvalues=NULL, names="bspl")

Arguments

Details

Spline functions are constructed by joining polynomials end-to-end at argument values called break points or knots. First, the interval is subdivided into a set of adjoining intervals separated the knots. Then a polynomial of order $m$ (degree $m-1$) is defined for each interval. To make the resulting piecewise polynomial smooth, two adjoining polynomials are constrained to have their values and all their derivatives up to order $m-2$ match at the point where they join.

Consider as an illustration the very common case where the order is 4 for all polynomials, so that degree of each polynomials is 3. That is, the polynomials are cubic. Then at each break point or knot, the values of adjacent polynomials must match, and so also for their first and second derivatives. Only their third derivatives will differ at the point of junction.

The number of degrees of freedom of a cubic spline function of this nature is calculated as follows. First, for the first interval, there are four degrees of freedom. Then, for each additional interval, the polynomial over that interval now has only one degree of freedom because of the requirement for matching values and derivatives. This means that the number of degrees of freedom is the number of interior knots (that is, not counting the lower and upper limits) plus the order of the polynomials:

nbasis = norder + length(breaks) - 2

The consistency of the values of nbasis, norder and breaks is checked, and an error message results if this equation is not satisfied.

B-splines are a set of special spline functions that can be used to construct any such piecewise polynomial by computing the appropriate linear combination. They derive their computational convenience from the fact that any B-spline basis function is nonzero over at most m adjacent intervals. The number of basis functions is given by the rule above for the number of degrees of freedom.

The number of intervals controls the flexibility of the spline; the more knots, the more flexible the resulting spline will be. But the position of the knots also plays a role. Where do we position the knots? There is room for judgment here, but two considerations must be kept in mind: (1) you usually want at least one argument value between two adjacent knots, and (2) there should be more knots where the curve needs to have sharp curvatures such as a sharp peak or valley or an abrupt change of level, but only a few knots are required where the curve is changing very slowly.

This function automatically includes norder replicates of the end points rangeval. By contrast, the analogous functions splineDesign and spline.des in the splines package do NOT automatically replicate the end points. To compare answers, the end knots must be replicated manually when using splineDesign or spline.des.

Value

a basis object of the type bspline

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

basisfd, create.constant.basis, create.exponential.basis, create.fourier.basis, create.monomial.basis, create.polygonal.basis, create.power.basis splineDesign spline.des

Examples

##
## The simplest basis currently available with this function:
##
bspl1.1 <- create.bspline.basis(norder=1)
oldpar <- par(no.readonly=TRUE)
plot(bspl1.1)
# 1 basis function, order 1 = degree 0 = step function:
# should be the same as above:
b1.1 <- create.bspline.basis(0:1, nbasis=1, norder=1, breaks=0:1)

all.equal(bspl1.1, b1.1)

bspl2.2 <- create.bspline.basis(norder=2)
plot(bspl2.2)
bspl3.3 <- create.bspline.basis(norder=3)
plot(bspl3.3)
bspl4.4 <- create.bspline.basis()
plot(bspl4.4)
bspl1.2 <- create.bspline.basis(norder=1, breaks=c(0,.5, 1))
plot(bspl1.2)
# 2 bases, order 1 = degree 0 = step functions:
# (1) constant 1 between 0 and 0.5 and 0 otherwise
# (2) constant 1 between 0.5 and 1 and 0 otherwise.
bspl2.3 <- create.bspline.basis(norder=2, breaks=c(0,.5, 1))
plot(bspl2.3)
# 3 bases:  order 2 = degree 1 = linear
# (1) line from (0,1) down to (0.5, 0), 0 after
# (2) line from (0,0) up to (0.5, 1), then down to (1,0)
# (3) 0 to (0.5, 0) then up to (1,1).
bspl3.4 <- create.bspline.basis(norder=3, breaks=c(0,.5, 1))
plot(bspl3.4)
# 4 bases:  order 3 = degree 2 = parabolas.
# (1) (x-.5)^2 from 0 to .5, 0 after
# (2) 2*(x-1)^2 from .5 to 1, and a parabola
#     from (0,0 to (.5, .5) to match
# (3 & 4) = complements to (2 & 1).
bSpl4. <- create.bspline.basis(c(-1,1))
plot(bSpl4.)
# Same as bSpl4.23 but over (-1,1) rather than (0,1).
# set up the b-spline basis for the lip data, using 23 basis functions,
#   order 4 (cubic), and equally spaced knots.
#  There will be 23 - 4 = 19 interior knots at 0.05, ..., 0.95
lipbasis <- create.bspline.basis(c(0,1), 23)
plot(lipbasis)
bSpl.growth <- create.bspline.basis(growth$age)
# cubic spline (order 4)
bSpl.growth6 <- create.bspline.basis(growth$age,norder=6)
# quintic spline (order 6)
##
## Nonnumeric rangeval
##
# Date
July4.1776 <- as.Date('1776-07-04')
Apr30.1789 <- as.Date('1789-04-30')
AmRev <- c(July4.1776, Apr30.1789)
BspRevolution <- create.bspline.basis(AmRev)
# POSIXct
July4.1776ct <- as.POSIXct1970('1776-07-04')
Apr30.1789ct <- as.POSIXct1970('1789-04-30')
AmRev.ct <- c(July4.1776ct, Apr30.1789ct)
BspRev.ct <- create.bspline.basis(AmRev.ct)
par(oldpar)

Create a Constant Basis

Description

Create a constant basis object, defining a single basis function whose value is everywhere 1.0.

Usage

create.constant.basis(rangeval=c(0, 1), names="const", axes=NULL)

Arguments

Value

a basis object with type component const.

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

basisfd, create.bspline.basis, create.exponential.basis, create.fourier.basis, create.monomial.basis, create.polygonal.basis, create.power.basis

Examples

basisobj <- create.constant.basis(c(-1,1))

Create an Exponential Basis

Description

Create an exponential basis object defining a set of exponential functions with rate constants in argument ratevec.

Usage

create.exponential.basis(rangeval=c(0,1), nbasis=NULL, ratevec=NULL,
                         dropind=NULL, quadvals=NULL, values=NULL,
                         basisvalues=NULL, names='exp', axes=NULL)

Arguments

Details

Exponential functions are of the type $exp(bx)$ where $b$ is the rate constant. If $b = 0$, the constant function is defined.

Value

a basis object with the type expon.

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

basisfd, create.bspline.basis, create.constant.basis, create.fourier.basis, create.monomial.basis, create.polygonal.basis, create.power.basis

Examples

#  Create an exponential basis over interval [0,5]
#  with basis functions 1, exp(-t) and exp(-5t)
basisobj <- create.exponential.basis(c(0,5),3,c(0,-1,-5))
#  plot the basis
oldpar <- par(no.readonly=TRUE)
plot(basisobj)
par(oldpar)

Create a Fourier Basis

Description

Create an Fourier basis object defining a set of Fourier functions with specified period.

Usage

create.fourier.basis(rangeval=c(0, 1), nbasis=3,
              period=diff(rangeval), dropind=NULL, quadvals=NULL,
              values=NULL, basisvalues=NULL, names=NULL,
              axes=NULL)

Arguments

Details

Functional data objects are constructed by specifying a set of basis functions and a set of coefficients defining a linear combination of these basis functions. The Fourier basis is a system that is usually used for periodic functions. It has the advantages of very fast computation and great flexibility. If the data are considered to be nonperiod, the Fourier basis is usually preferred. The first Fourier basis function is the constant function. The remainder are sine and cosine pairs with integer multiples of the base period. The number of basis functions generated is always odd.

Value

a basis object with the type fourier.

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

basisfd, create.bspline.basis, create.constant.basis, create.exponential.basis, create.monomial.basis, create.polygonal.basis, create.power.basis

Examples

# Create a minimal Fourier basis for annual data
#  using 3 basis functions
yearbasis3 <- create.fourier.basis(c(0,365),
                    axes=list("axesIntervals") )
#  plot the basis
oldpar <- par(no.readonly=TRUE)
plot(yearbasis3)

# Identify the months with letters
plot(yearbasis3, axes=list('axesIntervals', labels=monthLetters))

# The same labels as part of the basis object
yearbasis3. <- create.fourier.basis(c(0,365),
       axes=list("axesIntervals", labels=monthLetters) )
plot(yearbasis3.)

# set up the Fourier basis for the monthly temperature data,
#  using 9 basis functions with period 12 months.
monthbasis <- create.fourier.basis(c(0,12), 9, 12.0)

#  plot the basis
plot(monthbasis)

# Create a false Fourier basis using 1 basis function.
falseFourierBasis <- create.fourier.basis(nbasis=1)
#  plot the basis:  constant
plot(falseFourierBasis)
par(oldpar)

Create a Monomial Basis

Description

Creates a set of basis functions consisting of powers of the argument.

Usage

create.monomial.basis(rangeval=c(0, 1), nbasis=NULL,
         exponents=NULL, dropind=NULL, quadvals=NULL,
         values=NULL, basisvalues=NULL, names='monomial',
         axes=NULL)

Arguments

Value

a basis object with the type monom.

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

basisfd, link{create.basis} create.bspline.basis, create.constant.basis, create.fourier.basis, create.exponential.basis, create.polygonal.basis, create.power.basis

Examples

##
## simplest example: one constant 'basis function'
##
m0 <- create.monomial.basis(nbasis=1)
plot(m0)

##
## Create a monomial basis over the interval [-1,1]
##  consisting of the first three powers of t
##
basisobj <- create.monomial.basis(c(-1,1), 5)
#  plot the basis
oldpar <- par(no.readonly=TRUE)
plot(basisobj)

##
## rangeval of class Date or POSIXct
##
# Date
invasion1 <- as.Date('1775-09-04')
invasion2 <- as.Date('1812-07-12')
earlyUS.Canada <- c(invasion1, invasion2)
BspInvade1 <- create.monomial.basis(earlyUS.Canada)

# POSIXct
AmRev.ct <- as.POSIXct1970(c('1776-07-04', '1789-04-30'))
BspRev1.ct <- create.monomial.basis(AmRev.ct)
par(oldpar)

Create a Polygonal Basis

Description

A basis is set up for constructing polygonal lines, consisting of straight line segments that join together.

Usage

create.polygonal.basis(rangeval=NULL, argvals=NULL, dropind=NULL,
        quadvals=NULL, values=NULL, basisvalues=NULL, names='polygon',
        axes=NULL)

Arguments

Details

The actual basis functions consist of triangles, each with its apex over an argument value. Note that in effect the polygonal basis is identical to a B-spline basis of order 2 and a knot or break value at each argument value. The range of the polygonal basis is set to the interval defined by the smallest and largest argument values.

Value

a basis object with the type polyg.

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

basisfd, create.bspline.basis, create.constant.basis, create.exponential.basis, create.fourier.basis, create.monomial.basis, create.power.basis

Examples

#  Create a polygonal basis over the interval [0,1]
#  with break points at 0, 0.1, ..., 0.95, 1
(basisobj <- create.polygonal.basis(seq(0,1,0.1)))
#  plot the basis
oldpar <- par(no.readonly=TRUE)
plot(basisobj)
par(oldpar)

Create a Power Basis Object

Description

The basis system is a set of powers of argument $x$. That is, a basis function would be x^exponent, where exponent is a vector containing a set of powers or exponents. The power basis would normally only be used for positive values of x, since the power of a negative number is only defined for nonnegative integers, and the exponents here can be any real numbers.

Usage

create.power.basis(rangeval=c(0, 1), nbasis=NULL, exponents=NULL,
            dropind=NULL, quadvals=NULL, values=NULL,
            basisvalues=NULL, names='power', axes=NULL)

Arguments

Details

The power basis differs from the monomial basis in two ways. First, the powers may be nonintegers. Secondly, they may be negative. Consequently, a power basis is usually used with arguments that only take positive values, although a zero value can be tolerated if none of the powers are negative.

Value

a basis object of type power.

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

basisfd, create.bspline.basis, create.constant.basis, create.exponential.basis, create.fourier.basis, create.monomial.basis, create.polygonal.basis,

Examples

#  Create a power basis over the interval [1e-7,1]
#  with powers or exponents -1, -0.5, 0, 0.5 and 1
basisobj <- create.power.basis(c(1e-7,1), 5, seq(-1,1,0.5))
#  plot the basis
oldpar <- par(no.readonly=TRUE)
plot(basisobj)
par(oldpar)

Continuously Stirred Tank Reactor

Description

Functions for solving the Continuously Stirred Tank Reactor (CSTR) Ordinary Differential Equations (ODEs). A solution for observations where metrology error is assumed to be negligible can be obtained via lsoda(y, Time, CSTR2, parms); CSTR2 calls CSTR2in. When metrology error can not be ignored, use CSTRfn (which calls CSTRfitLS). To estimate parameters in the CSTR differential equation system (kref, EoverR, a, and / or b), pass CSTRres to nls. If nls fails to converge, first use optim or nlminb with CSTRsse, then pass the estimates to nls.

Usage

CSTR2in(Time, condition =
   c('all.cool.step', 'all.hot.step', 'all.hot.ramp', 'all.cool.ramp',
     'Tc.hot.exponential', 'Tc.cool.exponential', 'Tc.hot.ramp',
     'Tc.cool.ramp', 'Tc.hot.step', 'Tc.cool.step'),
   tau=1)
CSTR2(Time, y, parms)

CSTRfitLS(coef, datstruct, fitstruct, lambda, gradwrd=FALSE)
CSTRfn(parvec, datstruct, fitstruct, CSTRbasis, lambda, gradwrd=TRUE)
CSTRres(kref=NULL, EoverR=NULL, a=NULL, b=NULL,
        datstruct, fitstruct, CSTRbasis, lambda, gradwrd=FALSE)
CSTRsse(par, datstruct, fitstruct, CSTRbasis, lambda)

Arguments

Details

Ramsay et al. (2007) considers the following differential equation system for a continuously stirred tank reactor (CSTR):

dC/dt = (-betaCC(T, F.in)*C + F.in*C.in)

dT/dt = (-betaTT(Fcvec, F.in)*T + betaTC(T, F.in)*C + alpha(Fcvec)*T.co)

where

betaCC(T, F.in) = kref*exp(-1e4*EoverR*(1/T - 1/Tref)) + F.in

betaTT(Fcvec, F.in) = alpha(Fcvec) + F.in

betaTC(T, F.in) = (-delH/(rho*Cp))*betaCC(T, F.in)

alpha(Fcvec) = (a*Fcvec^(b+1) / (K1*(Fcvec + K2*Fcvec^b)))

K1 = V*rho*Cp

K2 = 1/(2*rhoc*Cpc)

The four functions CSTR2in, CSTR2, CSTRfitLS, and CSTRfn compute coefficients of basis vectors for two different solutions to this set of differential equations. Functions CSTR2in and CSTR2 work with 'lsoda' to provide a solution to this system of equations. Functions CSTSRitLS and CSTRfn are used to estimate parameters to fit this differential equation system to noisy data. These solutions are conditioned on specified values for kref, EoverR, a, and b. The other function, CSTRres, support estimation of these parameters using 'nls'.

CSTR2in translates a character string 'condition' into a data.frame containing system inputs for which the reaction of the system is desired. CSTR2 calls CSTR2in and then computes the corresponding predicted first derivatives of CSTR system outputs according to the right hand side of the system equations. CSTR2 can be called by 'lsoda' in the 'deSolve' package to actually solve the system of equations. To solve the CSTR equations for another set of inputs, the easiest modification might be to change CSTR2in to return the desired inputs. Another alternative would be to add an argument 'input.data.frame' that would be used in place of CSTR2in when present.

CSTRfitLS computes standardized residuals for systems outputs Conc, Temp or both as specified by fitstruct[["fit"]], a logical vector of length 2. The standardization is sqrt(datstruct[["Cwt"]]) and / or sqrt(datstruct[["Twt"]]) for Conc and Temp, respectively. CSTRfitLS also returns standardized deviations from the predicted first derivatives for Conc and Temp.

CSTRfn uses a Gauss-Newton optimization to estimates the coefficients of CSTRbasis to minimize the weighted sum of squares of residuals returned by CSTRfitLS.

CSTRres provides an interface between 'nls' and 'CSTRfn'. It gets the parameters to be estimated via the official function arguments, kref, EoverR, a, and / or b. The subset of these parameters to estimate must be specified both directly in the function call to 'nls' and indirectly via fitstruct[["estimate"]]. CSTRres gets the other CSTRfn arguments (datstruct, fitstruct, CSTRbasis, and lambda) via the 'data' argument of 'nls'.

CSTRsse computes sum of squares of residuals for use with optim or nlminb.

Value

References

Ramsay, J. O., Hooker, G., Cao, J. and Campbell, D. (2007) Parameter estimation for differential equations: A generalized smoothing approach (with discussion). Journal of the Royal Statistical Society, Series B, 69, 741-796.

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

lsoda nls

Examples

###
###
### 1.  lsoda(y, times, func=CSTR2, parms=...)
###
###
#  The system of two nonlinear equations has five forcing or
#  input functions.
#  These equations are taken from
#  Marlin, T. E. (2000) Process Control, 2nd Edition, McGraw Hill,
#  pages 899-902.
##
##  Set up the problem
##
fitstruct <- list(V    = 1.0,#  volume in cubic meters
                  Cp   = 1.0,#  concentration in cal/(g.K)
                  rho  = 1.0,#  density in grams per cubic meter
                  delH = -130.0,# cal/kmol
                  Cpc  = 1.0,#  concentration in cal/(g.K)
                  rhoc = 1.0,#  cal/kmol
                  Tref = 350)#  reference temperature
#  store true values of known parameters
EoverRtru = 0.83301#   E/R in units K/1e4
kreftru   = 0.4610 #   reference value
atru      = 1.678#     a in units (cal/min)/K/1e6
btru      = 0.5#       dimensionless exponent

#% enter these parameter values into fitstruct

fitstruct[["kref"]]   = kreftru#
fitstruct[["EoverR"]] = EoverRtru#  kref = 0.4610
fitstruct[["a"]]      = atru#       a in units (cal/min)/K/1e6
fitstruct[["b"]]      = btru#       dimensionless exponent

Tlim  = 64#    reaction observed over interval [0, Tlim]
delta = 1/12#  observe every five seconds
tspan = seq(0, Tlim, delta)#

coolStepInput <- CSTR2in(tspan, 'all.cool.step')

#  set constants for ODE solver

#  cool condition solution
#  initial conditions

Cinit.cool = 1.5965#  initial concentration in kmol per cubic meter
Tinit.cool = 341.3754# initial temperature in deg K
yinit = c(Conc = Cinit.cool, Temp=Tinit.cool)

#  load cool input into fitstruct

fitstruct[["Tcin"]] = coolStepInput[, "Tcin"];

#  solve  differential equation with true parameter values

if (require(deSolve)) {
coolStepSoln <- lsoda(y=yinit, times=tspan, func=CSTR2,
  parms=list(fitstruct=fitstruct, condition='all.cool.step', Tlim=Tlim) )
}
###
###
### 2.  CSTRfn
###
###

# See the script in '~R\library\fda\scripts\CSTR\CSTR_demo.R'
#  for more examples.

Compute a Cumulative Distribution Functional Data Object

Description

Function smooth.morph() maps a sorted set of variable values inside a closed interval into a set of equally-spaced probabilities in [0,1].

Usage

cumfd(xrnd, xrng, nbreaks=7, nfine=101)

Arguments

Details

Only the values of x within the interior of xrng are used in order to avoid distortion due to boundary inflation or deflation.

Value

A named list of length 2 containing:

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

smooth.morph, landmarkreg, register.fd

Examples

#  see the use of smooth.morph in landmarkreg.R
xrnd <- rbeta(50, 2, 5)
xrng <- c(0,1)
hist(xrnd)
range(xrnd)
cdfd <- cumfd(xrnd, xrng)

Plot Cycles for a Periodic Bivariate Functional Data Object

Description

A plotting function for data such as the knee-hip angles in the gait data or temperature-precipitation curves for the weather data.

Usage

cycleplot.fd(fdobj, matplt=TRUE, nx=201, ...)

Arguments

Value

None

Side Effects

A plot of the cycles

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

plot.fd, plotfit.fd,

Create smooth functions that fit scatterplot data.

Description

The function converts scatter plot data into one or a set of curves that do a satisfactory job of representing their corresponding abscissa and ordinate values by smooth curves. It returns a list object of class fdSmooth that contains curves defined by functional data objects of class fd as well as a variety of other results related to the data smoothing process.

The function tries to do as much for the user as possible before setting up a call to function smooth.basisPar. This is achieved by an initial call to function argvalsSwap that examines arguments argvals that contains abscissa values and y that contains ordinate values.

If the arguments are provided in an order that is not the one above, Data2fd attempts to swap argument positions to provide the correct order. A warning message is returned if this swap takes place. Any such automatic decision, though, has the possibility of being wrong, and the results should be carefully checked.

Preferably, the order of the arguments should be respected: argvals comes first and y comes second.

Be warned that the result that the fdSmooth object that is returned may not define a satisfactory smooth of the data, and consequently that it may be necessary to use the more sophisticated data smoothing function smooth.basis instead. Its help file provides a great deal more information than is provided here.

Usage

  Data2fd(argvals=NULL, y=NULL, basisobj=NULL, nderiv=NULL,
          lambda=3e-8/diff(as.numeric(range(argvals))),
          fdnames=NULL, covariates=NULL, method="chol")

Arguments

Details

This function tends to be used in rather simple applications where there is no need to control the roughness of the resulting curve with any great finesse. The roughness is essentially controlled by how many basis functions are used. In more sophisticated applications, it would be better to use the function smooth.basisPar.

It may happen that a value in argvals is outside the basis object interval due to rounding error and other causes of small violations. The code tests for this and pulls these near values into the interval if they are within 1e-7 times the interval width.

Value

an S3 list object of class fdSmooth defined in file smooth.basis1.R containing:

fd

the functional data object defining smooth B-spline curves that fit the data

df

the degrees of freedom in the fits

gcv

a generalized cross-validation value

beta

the semi-parametric coefficient array

SSE

Error sum of squares

penmat

the symmetric matrix defining roughness penalty

argvals

the argument values

y

the array of ordinate values in the scatter-plot

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

smooth.basisPar, smooth.basis, smooth.basis1, project.basis, smooth.fd, smooth.monotone, smooth.pos, day.5

Examples

##
## Simplest possible example:  constant function
##
# 1 basis, order 1 = degree 0 = constant function
b1.1 <- create.bspline.basis(nbasis=1, norder=1)
# data values: 1 and 2, with a mean of 1.5
y12 <- 1:2
# smooth data, giving a constant function with value 1.5
fd1.1 <- Data2fd(y12, basisobj=b1.1)
oldpar <- par(no.readonly=TRUE)
plot(fd1.1)
# now repeat the analysis with some smoothing, which moves the
# toward 0.
fd1.5 <- Data2fd(y12, basisobj=b1.1, lambda=0.5)
#  values of the smooth:
# fd1.1 = sum(y12)/(n+lambda*integral(over arg=0 to 1 of 1))
#         = 3 / (2+0.5) = 1.2
#. JR ... Data2fd returns an fdsmooth object and a member of the
#         is an fd smooth object.  Calls to functions expecting 
#         an fd object require attaching $fd to the fdsmooth object
#         this is required in lines 268, 311 and 337
eval.fd(seq(0, 1, .2), fd1.5)
##
## step function smoothing
##
# 2 step basis functions: order 1 = degree 0 = step functions
b1.2 <- create.bspline.basis(nbasis=2, norder=1)
#  fit the data without smoothing
fd1.2 <- Data2fd(1:2, basisobj=b1.2)
# plot the result:  A step function:  1 to 0.5, then 2
op <- par(mfrow=c(2,1))
plot(b1.2, main='bases')
plot(fd1.2, main='fit')
par(op)
##
## Simple oversmoothing
##
# 3 step basis functions: order 1 = degree 0 = step functions
b1.3 <- create.bspline.basis(nbasis=3, norder=1)
#  smooth the data with smoothing
fd1.3 <- Data2fd(y12, basisobj=b1.3, lambda=0.5)
#  plot the fit along with the points
plot(0:1, c(0, 2), type='n')
points(0:1, y12)
lines(fd1.3)
# Fit = penalized least squares with penalty =
#          = lambda * integral(0:1 of basis^2),
#            which shrinks the points towards 0.
# X1.3 = matrix(c(1,0, 0,0, 0,1), 2)
# XtX = crossprod(X1.3) = diag(c(1, 0, 1))
# penmat = diag(3)/3
#        = 3x3 matrix of integral(over arg=0:1 of basis[i]*basis[j])
# Xt.y = crossprod(X1.3, y12) = c(1, 0, 2)
# XtX + lambda*penmat = diag(c(7, 1, 7)/6
# so coef(fd1.3.5) = solve(XtX + lambda*penmat, Xt.y)
#                  = c(6/7, 0, 12/7)
##
## linear spline fit
##
# 3 bases, order 2 = degree 1
b2.3 <- create.bspline.basis(norder=2, breaks=c(0, .5, 1))
# interpolate the values 0, 2, 1
fd2.3 <- Data2fd(c(0,2,1), basisobj=b2.3, lambda=0)
#  display the coefficients
round(fd2.3$coefs, 4)
# plot the results
op <- par(mfrow=c(2,1))
plot(b2.3, main='bases')
plot(fd2.3, main='fit')
par(op)
# apply some smoothing
fd2.3. <- Data2fd(c(0,2,1), basisobj=b2.3, lambda=1)
op <- par(mfrow=c(2,1))
plot(b2.3, main='bases')
plot(fd2.3., main='fit', ylim=c(0,2))
par(op)
all.equal(
 unclass(fd2.3)[-1], 
 unclass(fd2.3.)[-1])
##** CONCLUSION:  
##** The only differences between fd2.3 and fd2.3.
##** are the coefficients, as we would expect.  

##
## quadratic spline fit
##
# 4 bases, order 3 = degree 2 = continuous, bounded, locally quadratic
b3.4 <- create.bspline.basis(norder=3, breaks=c(0, .5, 1))
# fit values c(0,4,2,3) without interpolation
fd3.4 <- Data2fd(c(0,4,2,3), basisobj=b3.4, lambda=0)
round(fd3.4$coefs, 4)
op <- par(mfrow=c(2,1))
plot(b3.4)
plot(fd3.4)
points(c(0,1/3,2/3,1), c(0,4,2,3))
par(op)
#  try smoothing
fd3.4. <- Data2fd(c(0,4,2,3), basisobj=b3.4, lambda=1)
round(fd3.4.$coef, 4)
op <- par(mfrow=c(2,1))
plot(b3.4)
plot(fd3.4., ylim=c(0,4))
points(seq(0,1,len=4), c(0,4,2,3))
par(op)
##
##  Two simple Fourier examples
##
gaitbasis3 <- create.fourier.basis(nbasis=5)
gaitfd3    <- Data2fd(seq(0,1,len=20), gait, basisobj=gaitbasis3)
# plotfit.fd(gait, seq(0,1,len=20), gaitfd3)
#    set up the fourier basis
daybasis <- create.fourier.basis(c(0, 365), nbasis=65)
#  Make temperature fd object
#  Temperature data are in 12 by 365 matrix tempav
#    See analyses of weather data.
tempfd <- Data2fd(CanadianWeather$dailyAv[,,"Temperature.C"],
                  day.5, daybasis)
#  plot the temperature curves
par(mfrow=c(1,1))
plot(tempfd)
##
## argvals of class Date and POSIXct
##
#  These classes of time can generate very large numbers when converted to 
#  numeric vectors.  For basis systems such as polynomials or splines,
#  severe rounding error issues can arise if the time interval for the 
#  data is very large.  To offset this, it is best to normalize the
#  numeric version of the data before analyzing them.
#  Date class time unit is one day, divide by 365.25.
invasion1 <- as.Date('1775-09-04')
invasion2 <- as.Date('1812-07-12')
earlyUS.Canada <- as.numeric(c(invasion1, invasion2))/365.25
BspInvasion <- create.bspline.basis(earlyUS.Canada)
earlyYears  <- seq(invasion1, invasion2, length.out=7)
earlyQuad   <- (as.numeric(earlyYears-invasion1)/365.25)^2
earlyYears  <- as.numeric(earlyYears)/365.25
fitQuad <- Data2fd(earlyYears, earlyQuad, BspInvasion)
# POSIXct: time unit is one second, divide by 365.25*24*60*60
rescale     <- 365.25*24*60*60
AmRev.ct    <- as.POSIXct1970(c('1776-07-04', '1789-04-30'))
BspRev.ct   <- create.bspline.basis(as.numeric(AmRev.ct)/rescale)
AmRevYrs.ct <- seq(AmRev.ct[1], AmRev.ct[2], length.out=14)
AmRevLin.ct <- as.numeric(AmRevYrs.ct-AmRev.ct[1])
AmRevYrs.ct <- as.numeric(AmRevYrs.ct)/rescale
AmRevLin.ct <- as.numeric(AmRevLin.ct)/rescale
fitLin.ct   <- Data2fd(AmRevYrs.ct, AmRevLin.ct, BspRev.ct)
par(oldpar)

Numeric and character vectors to facilitate working with dates

Description

Numeric and character vectors to simplify functional data computations and plotting involving dates.

Format

dayOfYear

a numeric vector = 1:365 with names 'jan01' to 'dec31'.

dayOfYearShifted

a numeric vector = c(182:365, 1:181) with names 'jul01' to 'jun30'.

day.5

a numeric vector = dayOfYear-0.5 = 0.5, 1.5, ..., 364.5

daysPerMonth

a numeric vector of the days in each month (ignoring leap years) with names = month.abb

monthEnd

a numeric vector of cumsum(daysPerMonth) with names = month.abb

monthEnd.5

a numeric vector of the middle of the last day of each month with names = month.abb = c(Jan=30.5, Feb=58.5, ..., Dec=364.5)

monthBegin.5

a numeric vector of the middle of the first day of each month with names - month.abb = c(Jan=0.5, Feb=31.5, ..., Dec=334.5)

monthMid

a numeric vector of the middle of the month = (monthBegin.5 + monthEnd.5)/2

monthLetters

A character vector of c("j", "F", "m", "A", "M", "J", "J", "A", "S", "O", "N", "D"), with 'month.abb' as the names.

weeks

a numeric vector of length 53 marking 52 periods of approximately 7 days each throughout the year = c(0, 365/52, ..., 365)

Details

Miscellaneous vectors often used in 'fda' scripts.

Source

Ramsay, James O., and Silverman, Bernard W. (2006), Functional Data Analysis, 2nd ed., Springer, New York, pp. 5, 47-53.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

axisIntervals month.abb

Examples

daybasis65 <- create.fourier.basis(c(0, 365), 65)
daytempfd <- with(CanadianWeather, smooth.basisPar(day.5,
    dailyAv[,,"Temperature.C"], daybasis65)$fd )
oldpar <- par(axes=FALSE)
plot(daytempfd)
axisIntervals(1)
# axisIntervals by default uses
# monthBegin.5, monthEnd.5, monthMid, and month.abb
axis(2)
par(oldpar)

Compute a Probability Density Function

Description

Like the regular S-PLUS function density, this function computes a probability density function for a sample of values of a random variable. However, in this case the density function is defined by a functional parameter object WfdParobj along with a normalizing constant C.

The density function $p(x)$ has the form p(x) = C exp[W(x)] where function $W(x)$ is defined by the functional data object WfdParobj.

Usage

## S3 method for class 'fd'
density(x, WfdParobj, conv=0.0001, iterlim=20,
           active=1:nbasis, dbglev=0, ...)

Arguments

Details

The goal of the function is provide a smooth density function estimate that approaches some target density by an amount that is controlled by the linear differential operator Lfdobj and the penalty parameter. For example, if the second derivative of $W(t)$ is penalized heavily, this will force the function to approach a straight line, which in turn will force the density function itself to be nearly normal or Gaussian. Similarly, to each textbook density function there corresponds a $W(t)$, and to each of these in turn their corresponds a linear differential operator that will, when apply to $W(t)$, produce zero as a result. To plot the density function or to evaluate it, evaluate Wfdobj, exponentiate the resulting vector, and then divide by the normalizing constant C.

Value

a named list of length 4 containing:

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

intensity.fd

Examples

#  set up range for density
rangeval <- c(-3,3)
#  set up some standard normal data
x <- rnorm(50)
#  make sure values within the range
x[x < -3] <- -2.99
x[x >  3] <-  2.99
#  set up basis for W(x)
basisobj <- create.bspline.basis(rangeval, 11)
#  set up initial value for Wfdobj
Wfd0 <- fd(matrix(0,11,1), basisobj)
WfdParobj <- fdPar(Wfd0)
#  estimate density
denslist <- density.fd(x, WfdParobj)
#  plot density
oldpar <- par(no.readonly=TRUE)
xval <- seq(-3,3,.2)
wval <- eval.fd(xval, denslist$Wfdobj)
pval <- exp(wval)/denslist$C
plot(xval, pval, type="l", ylim=c(0,0.4))
points(x,rep(0,50))
par(oldpar)

Compute a Derivative of a Functional Data Object

Description

A derivative of a functional data object, or the result of applying a linear differential operator to a functional data object, is then converted to a functional data object. This is intended for situations where a derivative is to be manipulated as a functional data object rather than simply evaluated.

Usage

## S3 method for class 'fd'
deriv(expr, Lfdobj=int2Lfd(1), ...)

Arguments

Details

Typically, a derivative has more high frequency variation or detail than the function itself. The basis defining the function is used, and therefore this must have enough basis functions to represent the variation in the derivative satisfactorily.

Value

a functional data object for the derivative

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

getbasismatrix, eval.basis deriv

Examples

#  Estimate the acceleration functions for growth curves
#  See the analyses of the growth data.
#  Set up the ages of height measurements for Berkeley data
age <- c( seq(1, 2, 0.25), seq(3, 8, 1), seq(8.5, 18, 0.5))
#  Range of observations
rng <- c(1,18)
#  Set up a B-spline basis of order 6 with knots at ages
knots  <- age
norder <- 6
nbasis <- length(knots) + norder - 2
hgtbasis <- create.bspline.basis(rng, nbasis, norder, knots)
#  Set up a functional parameter object for estimating
#  growth curves.  The 4th derivative is penalyzed to
#  ensure a smooth 2nd derivative or acceleration.
Lfdobj <- 4
lambda <- 10^(-0.5)   #  This value known in advance.
growfdPar <- fdPar(fd(matrix(0,nbasis,1),hgtbasis), Lfdobj, lambda)
#  Smooth the data.  The data for the boys and girls
#  are in matrices hgtm and hgtf, respectively.
hgtmfd <- smooth.basis(age, growth$hgtm, growfdPar)$fd
hgtffd <- smooth.basis(age, growth$hgtf, growfdPar)$fd
#  Compute the acceleration functions
accmfd <- deriv.fd(hgtmfd, 2)
accffd <- deriv.fd(hgtffd, 2)
#  Plot the two sets of curves
oldpar <- par(mfrow=c(2,1))
plot(accmfd)
plot(accffd)
par(oldpar)

Degrees of Freedom for Residuals from a Functional Regression

Description

Effective degrees of freedom for residuals, being the trace of the idempotent hat matrix transforming observations into residuals from the fit.

Usage

## S3 method for class 'fRegress'
df.residual(object, ...)

Arguments

Details

1. Determine N = number of observations

2. df.model <- object$df

3. df.residual <- (N - df.model)

4. Add attributes

Value

The numeric value of the residual degrees-of-freedom extracted from object with the following attributes:

Author(s)

Spencer Graves

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York. Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York. Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York. Hastie, Trevor, Tibshirani, Robert, and Friedman, Jerome (2001) The Elements of Statistical Learning: Data Mining, Inference, and Prediction, Springer, New York.

See Also

fRegress df.residual

Convert Degrees of Freedom to a Smoothing Parameter Value

Description

The degree of roughness of an estimated function is controlled by a smoothing parameter $lambda$ that directly multiplies the penalty. However, it can be difficult to interpret or choose this value, and it is often easier to determine the roughness by choosing a value that is equivalent of the degrees of freedom used by the smoothing procedure. This function converts a degrees of freedom value into a multiplier $lambda$.

Usage

df2lambda(argvals, basisobj, wtvec=rep(1, n), Lfdobj=0,
          df=nbasis)

Arguments

Details

The conversion requires a one-dimensional optimization and may be therefore computationally intensive.

Value

a positive smoothing parameter value $lambda$

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

lambda2df, lambda2gcv

Examples

#  Smooth growth curves using a specified value of
#  degrees of freedom.
#  Set up the ages of height measurements for Berkeley data
age <- c( seq(1, 2, 0.25), seq(3, 8, 1), seq(8.5, 18, 0.5))
#  Range of observations
rng <- c(1,18)
#  Set up a B-spline basis of order 6 with knots at ages
knots  <- age
norder <- 6
nbasis <- length(knots) + norder - 2
hgtbasis <- create.bspline.basis(rng, nbasis, norder, knots)
#  Find the smoothing parameter equivalent to 12
#  degrees of freedom
lambda <- df2lambda(age, hgtbasis, df=12)
#  Set up a functional parameter object for estimating
#  growth curves.  The 4th derivative is penalyzed to
#  ensure a smooth 2nd derivative or acceleration.
Lfdobj <- 4
growfdPar <- fdPar(fd(matrix(0,nbasis,1),hgtbasis), Lfdobj, lambda)
#  Smooth the data.  The data for the girls are in matrix
#  hgtf.
hgtffd <- smooth.basis(age, growth$hgtf, growfdPar)$fd
#  Plot the curves
oldpar <- par(no.readonly=TRUE)
plot(hgtffd)
par(oldpar)

Get subdirectories

Description

If you want only subfolders and no files, use dirs. With recursive = FALSE, dir returns both folders and files. With recursive = TRUE, it returns only files.

Usage

dirs(path='.', pattern=NULL, exclude=NULL, all.files=FALSE,
     full.names=FALSE, recursive=FALSE, ignore.case=FALSE) 

Arguments

Details

1. mainDir <- dir(...) without recurse

2. Use file.info to restrict mainDir to only directories.

3. If !recursive, return the restricted mainDir. Else, if length(mainDir) > 0, create dirList to hold the results of the recursion and call dirs for each component of mainDir. Then unlist and return the result.

Value

A character vector of the desired subdirectories.

Author(s)

Spencer Graves

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

dir, file.info

Eigenanalysis preserving dimnames

Description

Compute eigenvalues and vectors, assigning names to the eigenvalues and dimnames to the eigenvectors.

Usage

Eigen(x, symmetric, only.values = FALSE, valuenames)

Arguments

Details

1. Check 'symmetric'

2. ev <- eigen(x, symmetric, only.values = FALSE, EISPACK = FALSE); see eigen for more details.

3. rNames = rownames(x); if this is NULL, rNames = if(symmetric) paste('x', 1:nrow(x), sep=”) else paste('xcol', 1:nrow(x)).

4. Parse 'valuenames', assign to names(ev[['values']]).

5. dimnames(ev[['vectors']]) <- list(rNames, valuenames)

NOTE: This naming convention is fairly obvious if 'x' is symmetric. Otherwise, dimensional analysis suggests problems with almost any naming convention. To see this, consider the following simple example:

X <- matrix(1:4, 2, dimnames=list(LETTERS[1:2], letters[3:4]))

X.inv <- solve(X)

One way of interpreting this is to assume that colnames are really reciprocals of the units. Thus, in this example, X[1,1] is in units of 'A/c' and X.inv[1,1] is in units of 'c/A'. This would make any matrix with the same row and column names potentially dimensionless. Since eigenvalues are essentially the diagonal of a diagonal matrix, this would mean that eigenvalues are dimensionless, and their names are merely placeholders.

Value

a list with components values and (if only.values = FALSE) vectors, as described in eigen.

Author(s)

Spencer Graves

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

eigen, svd qr chol

Examples

X <- matrix(1:4, 2, dimnames=list(LETTERS[1:2], letters[3:4]))
eigen(X)
Eigen(X)
Eigen(X, valuenames='eigval')

Y <- matrix(1:4, 2, dimnames=list(letters[5:6], letters[5:6]))
Eigen(Y)

Eigen(Y, symmetric=TRUE)
# only the lower triangle is used;
# the upper triangle is ignored.  

Stability Analysis for Principle Differential Analysis

Description

Performs a stability analysis of the result of pda.fd, returning the real and imaginary parts of the eigenfunctions associated with the linear differential operator.

Usage

eigen.pda(pdaList,plotresult=TRUE,npts=501,...)

Arguments

Details

Conducts an eigen decomposition of the linear differential equation implied by the result of pda.fd. Imaginary eigenvalues indicate instantaneous oscillatory behavior. Positive real eigenvalues indicate exponential increase, negative real eigenvalues correspond to exponential decay. If the principle differential analysis also included the estimation of a forcing function, the limiting stable points are also tracked.

Value

Returns a list with elements

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

pda.fd plot.pda.fd pda.overlay

Examples

#  A pda analysis of the handwriting data

# reduce the size to reduce the compute time for the example
ni <- 281
indx <- seq(1, 1401, length=ni)
fdaarray = handwrit[indx,,]
fdatime  <- seq(0, 2.3, len=ni)

#  basis for coordinates

fdarange <- c(0, 2.3)
breaks = seq(0,2.3,length.out=116)
norder = 6
fdabasis = create.bspline.basis(fdarange,norder=norder,breaks=breaks)
nbasis <- fdabasis$nbasis

#  parameter object for coordinates

fdaPar = fdPar(fd(matrix(0,nbasis,1),fdabasis),int2Lfd(4),1e-8)

#  coordinate functions and a list tontaining them

Xfd = smooth.basis(fdatime, fdaarray[,,1], fdaPar)$fd
Yfd = smooth.basis(fdatime, fdaarray[,,2], fdaPar)$fd

xfdlist = list(Xfd, Yfd)

#  basis and parameter object for weight functions

fdabasis2 = create.bspline.basis(fdarange,norder=norder,nbasis=31)
fdafd2    = fd(matrix(0,31,2),fdabasis2)
pdaPar    = fdPar(fdafd2,1,1e-8)

pdaParlist = list(pdaPar, pdaPar)

bwtlist = list( list(pdaParlist,pdaParlist), list(pdaParlist,pdaParlist) )

Predicting electricity demand in Adelaide from temperature

Description

The data sets used in this demonstration analysis consist of half-hourly electricity demands from Sunday to Saturday in Adelaide, Australia, between July 6, 1976 and March 31, 2007. Also provided in the same format and times are the temperatures in degrees Celsius at the Adelaide airport. The shapes of the demand curves for each day resemble those for temperature, and the goal is to see how well a concurrent functional regress model fit by function fRegress can fit the demand curves.

Format

There is a data object for each day of the week, and in each object, the member y, such as mondaydemand$y, is a 48 by 508 matrix of electricity demand values, one for each of 48 half-hourly points, and for each of 508 weeks.

Details

The demonstration is designed to show how to use the main functional regression function fRegress in order to fit a functional dependent variable (electricity demand) from an obviously important input or covariate functional variable (temperature). In the texts cited, this model is referred to as a concurrent functional regression because the model assumes that changes in the input functional variable directly and immediately change the dependent functional variable.

Also illustrated is the estimation of pointwise confidence intervals for the regression coefficient functions and for the fitting functions that approximate electricity demand. These steps involve functons fRegress.stderr and predict.fRegress.

The data are supplied through the CRAN system by package fds, and further information as well as many other datasets are to be found there.

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

Values of Basis Functions or their Derivatives

Description

A set of basis functions are evaluated at a vector of argument values. If a linear differential object is provided, the values are the result of applying the the operator to each basis function.

Usage

eval.basis(evalarg, basisobj, Lfdobj=0, returnMatrix=FALSE)
## S3 method for class 'basisfd'
predict(object, newdata=NULL, Lfdobj=0,
                          returnMatrix=FALSE, ...)

Arguments

Details

If a linear differential operator object is supplied, the basis must be such that the highest order derivative can be computed. If a B-spline basis is used, for example, its order must be one larger than the highest order of derivative required.

Value

a matrix of basis function values with rows corresponding to argument values and columns to basis functions.

predict.basisfd is a convenience wrapper for eval.basis.

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

getbasismatrix, eval.fd, plot.basisfd

Examples

##
## 1.  B-splines
##
# The simplest basis currently available:
# a single step function
bspl1.1 <- create.bspline.basis(norder=1, breaks=0:1)
eval.bspl1.1 <- eval.basis(seq(0, 1, .2), bspl1.1)

# check
eval.bspl1.1. <- matrix(rep(1, 6), 6,
       dimnames=list(NULL, 'bspl') )

all.equal(eval.bspl1.1, eval.bspl1.1.)


# The second simplest basis:
# 2 step functions, [0, .5], [.5, 1]
bspl1.2 <- create.bspline.basis(norder=1, breaks=c(0,.5, 1))
eval.bspl1.2 <- eval.basis(seq(0, 1, .2), bspl1.2)

# Second order B-splines (degree 1:  linear splines)
bspl2.3 <- create.bspline.basis(norder=2, breaks=c(0,.5, 1))
eval.bspl2.3 <- eval.basis(seq(0, 1, .1), bspl2.3)
# 3 bases:  order 2 = degree 1 = linear
# (1) line from (0,1) down to (0.5, 0), 0 after
# (2) line from (0,0) up to (0.5, 1), then down to (1,0)
# (3) 0 to (0.5, 0) then up to (1,1).

##
## 2.  Fourier
##
# The false Fourier series with 1 basis function
falseFourierBasis <- create.fourier.basis(nbasis=1)
eval.fFB <- eval.basis(seq(0, 1, .2), falseFourierBasis)

# Simplest real Fourier basis with 3 basis functions
fourier3 <- create.fourier.basis()
eval.fourier3 <- eval.basis(seq(0, 1, .2), fourier3)

# 3 basis functions on [0, 365]
fourier3.365 <- create.fourier.basis(c(0, 365))
eval.F3.365 <- eval.basis(day.5, fourier3.365)
oldpar <- par(no.readonly=TRUE)
matplot(eval.F3.365, type="l")

# The next simplest Fourier basis (5  basis functions)
fourier5 <- create.fourier.basis(nbasis=5)
eval.F5 <- eval.basis(seq(0, 1, .1), fourier5)
matplot(eval.F5, type="l")

# A more complicated example
dayrng <- c(0, 365)

nbasis <- 51
norder <- 6

weatherBasis <- create.fourier.basis(dayrng, nbasis)
basisMat <- eval.basis(day.5, weatherBasis)

matplot(basisMat[, 1:5], type="l")

##
## 3.  predict.basisfd
##
basisMat. <- predict(weatherBasis, day.5)

all.equal(basisMat, basisMat.)


##
## 4.  Date and POSIXct
##
# Date
July4.1776 <- as.Date('1776-07-04')
Apr30.1789 <- as.Date('1789-04-30')
AmRev <- c(July4.1776, Apr30.1789)
BspRevolution <- create.bspline.basis(AmRev)
AmRevYears <- seq(July4.1776, Apr30.1789, length.out=14)
AmRevBases <- predict(BspRevolution, AmRevYears)
matplot(AmRevYears, AmRevBases, type='b')
# Image is correct, but
# matplot does not recogize the Date class of x

# POSIXct
AmRev.ct <- as.POSIXct1970(c('1776-07-04', '1789-04-30'))
BspRev.ct <- create.bspline.basis(AmRev.ct)
AmRevYrs.ct <- seq(AmRev.ct[1], AmRev.ct[2], length.out=14)
AmRevBas.ct <- predict(BspRev.ct, AmRevYrs.ct)
matplot(AmRevYrs.ct, AmRevBas.ct, type='b')
# Image is correct, but
# matplot does not recognize the POSIXct class of x
par(oldpar)

Values a Two-argument Functional Data Object

Description

A vector of argument values for the first argument s of the functional data object to be evaluated.

Usage

eval.bifd(sevalarg, tevalarg, bifd, sLfdobj=0, tLfdobj=0)

Arguments

Value

an array of 2, 3, or 4 dimensions containing the function values. The first dimension corresponds to the argument values in sevalarg, the second to argument values in tevalarg, the third if present to replications, and the fourth if present to functions.

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

Examples

# every-other-day basis to save test time
daybasis   <- create.fourier.basis(c(0,365), 183)
harmLcoef  <- c(0,(2*pi/365)^2,0)
harmLfd    <- vec2Lfd(harmLcoef, c(0,365))
templambda <- 1.0
tempfdPar  <- fdPar(fd(matrix(0,daybasis$nbasis,1), daybasis), harmLfd, 
                    lambda=1)
tempfd     <- smooth.basis(day.5, CanadianWeather$dailyAv[,,"Temperature.C"], 
                           tempfdPar)$fd
#    define the variance-covariance bivariate fd object
tempvarbifd <- var.fd(tempfd)
#    evaluate the variance-covariance surface and plot
weektime    <- seq(0,365,len=53)
tempvarmat  <- eval.bifd(weektime,weektime,tempvarbifd)
#    make a perspective plot of the variance function
oldpar <- par(no.readonly=TRUE)
persp(tempvarmat)
par(oldpar)

Values of a Functional Data Object

Description

Evaluate a functional data object at specified argument values, or evaluate a derivative or the result of applying a linear differential operator to the functional object.

Usage

eval.fd(evalarg, fdobj, Lfdobj=0, returnMatrix=FALSE)
## S3 method for class 'fd'
predict(object, newdata=NULL, Lfdobj=0, returnMatrix=FALSE,
                     ...)
## S3 method for class 'fdPar'
predict(object, newdata=NULL, Lfdobj=0,
                     returnMatrix=FALSE, ...)
## S3 method for class 'fdSmooth'
predict(object, newdata=NULL, Lfdobj=0,
                     returnMatrix=FALSE, ...)
## S3 method for class 'fdSmooth'
fitted(object, returnMatrix=FALSE, ...)
## S3 method for class 'fdSmooth'
residuals(object, returnMatrix=FALSE, ...)

Arguments

Details

eval.fd evaluates Lfdobj of fdobj at evalarg.

predict.fd is a convenience wrapper for eval.fd. If newdata is NULL and fdobj[['basis']][['type']] is bspline, newdata = unique(knots(fdojb,interior=FALSE)); otherwise, newdata = fdobj[['basis']][['rangeval']].

predict.fdSmooth, fitted.fdSmooth and residuals.fdSmooth are other wrappers for eval.fd.

Value

an array of 2 or 3 dimensions containing the function values. The first dimension corresponds to the argument values in evalarg, the second to replications, and the third if present to functions.

Author(s)

Soren Hosgaard wrote an initial version of predict.fdSmooth, fitted.fdSmooth, and residuals.fdSmooth.

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

getbasismatrix, eval.bifd, eval.penalty, eval.monfd, eval.posfd

Examples

oldpar <- par(no.readonly=TRUE)
##
## eval.fd
##
#    set up the fourier basis
daybasis <- create.fourier.basis(c(0, 365), nbasis=65)
#  Make temperature fd object
#  Temperature data are in 12 by 365 matrix tempav
#  See analyses of weather data.
#  Set up sampling points at mid days
#  Convert the data to a functional data object
tempfd <- smooth.basis(day.5,  CanadianWeather$dailyAv[,,"Temperature.C"],
                       daybasis)$fd
#   set up the harmonic acceleration operator
Lbasis  <- create.constant.basis(c(0, 365))
Lcoef   <- matrix(c(0,(2*pi/365)^2,0),1,3)
bfdobj  <- fd(Lcoef,Lbasis)
bwtlist <- fd2list(bfdobj)
harmaccelLfd <- Lfd(3, bwtlist)
#   evaluate the value of the harmonic acceleration
#   operator at the sampling points
Ltempmat <- eval.fd(day.5, tempfd, harmaccelLfd)

#  Confirm that it still works with
#  evalarg = a matrix with only one column
#  when fdobj[['coefs']] is a matrix with multiple columns

Ltempmat. <- eval.fd(matrix(day.5, ncol=1), tempfd, harmaccelLfd)
#  confirm that the two answers are the same


all.equal(Ltempmat, Ltempmat.)


#  Plot the values of this operator
matplot(day.5, Ltempmat, type="l")

##
## predict.fd
##
predict(tempfd) # end points only at 35 locations
str(predict(tempfd, day.5)) # 365 x 35 matrix
str(predict(tempfd, day.5, harmaccelLfd))

# cublic splie with knots at 0, .5, 1
bspl3 <- create.bspline.basis(c(0, .5, 1))
plot(bspl3) # 5 bases
fd.bspl3 <- fd(c(0, 0, 1, 0, 0), bspl3)
pred3 <- predict(fd.bspl3)

pred3. <- matrix(c(0, .5, 0), 3)
dimnames(pred3.) <- list(NULL, 'reps 1')

all.equal(pred3, pred3.)


pred.2 <- predict(fd.bspl3, c(.2, .8))

pred.2. <- matrix(.176, 2, 1)
dimnames(pred.2.) <- list(NULL, 'reps 1')

all.equal(pred.2, pred.2.)


##
## predict.fdSmooth
##
lipSm9 <- smooth.basisPar(liptime, lip, lambda=1e-9)$fd
plot(lipSm9)

##
## with evalarg of class Date and POSIXct
##
# Date
July4.1776 <- as.Date('1776-07-04')
Apr30.1789 <- as.Date('1789-04-30')
AmRev <- c(July4.1776, Apr30.1789)
BspRevolution <- create.bspline.basis(AmRev)

AmRevYears <- seq(July4.1776, Apr30.1789, length.out=14)
(AmRevLinear <- as.numeric(AmRevYears-July4.1776))
fitLin <- smooth.basis(AmRevYears, AmRevLinear, BspRevolution)
AmPred <- predict(fitLin, AmRevYears)

# POSIXct
AmRev.ct <- as.POSIXct1970(c('1776-07-04', '1789-04-30'))
BspRev.ct <- create.bspline.basis(AmRev.ct)
AmRevYrs.ct <- seq(AmRev.ct[1], AmRev.ct[2], length.out=14)
(AmRevLin.ct <- as.numeric(AmRevYrs.ct-AmRev.ct[2]))
fitLin.ct <- smooth.basis(AmRevYrs.ct, AmRevLin.ct, BspRev.ct)
AmPred.ct <- predict(fitLin.ct, AmRevYrs.ct)
par(oldpar)

Values of a Monotone Functional Data Object

Description

Evaluate a monotone functional data object at specified argument values, or evaluate a derivative of the functional object.

Usage

eval.monfd(evalarg, Wfdobj, Lfdobj=int2Lfd(0), returnMatrix=FALSE)
## S3 method for class 'monfd'
predict(object, newdata=NULL, Lfdobj=0, returnMatrix=FALSE, ...)
## S3 method for class 'monfd'
fitted(object, ...)
## S3 method for class 'monfd'
residuals(object, ...)

Arguments

Details

A monotone function data object $h(t)$ is defined by $h(t) = [D^{-1} exp Wfdobj](t)$. In this equation, the operator $D^{-1}$ means taking the indefinite integral of the function to which it applies. Note that this equation implies that the monotone function has a value of zero at the lower limit of the arguments. To actually fit monotone data, it will usually be necessary to estimate an intercept and a regression coefficient to be applied to $h(t)$, usually with the least squares regression function lsfit. The function Wfdobj that defines the monotone function is usually estimated by monotone smoothing function smooth.monotone.

eval.monfd only computes the standardized monotone form. predict.monfd computes the scaled version using with(object, beta[1] + beta[2]*eval.monfd(...)) if Lfdobj = 0 or beta[2]*eval.monfd(...) if Lfdobj > 0.

Value

a matrix containing the monotone function values. The first dimension corresponds to the argument values in evalarg and the second to replications.

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

eval.fd, smooth.monotone eval.posfd

Examples

oldpar <- par(no.readonly=TRUE)
#  Estimate the acceleration functions for growth curves
#  See the analyses of the growth data.
#  Set up the ages of height measurements for Berkeley data
age <- c( seq(1, 2, 0.25), seq(3, 8, 1), seq(8.5, 18, 0.5))
#  Range of observations
rng <- c(1,18)
#  First set up a basis for monotone smooth
#  We use b-spline basis functions of order 6
#  Knots are positioned at the ages of observation.
norder <- 6
nage   <- length(age)
nbasis <- nage + norder - 2
wbasis <- create.bspline.basis(rng, nbasis, norder, age)
#  starting values for coefficient
cvec0 <- matrix(0,nbasis,1)
Wfd0  <- fd(cvec0, wbasis)
#  set up functional parameter object
Lfdobj    <- 3          #  penalize curvature of acceleration
lambda    <- 10^(-0.5)  #  smoothing parameter
growfdPar <- fdPar(Wfd0, Lfdobj, lambda)
#  Smooth the data for the first girl
hgt1 <- growth$hgtf[,1]
#   set conv = 0.1 and iterlim=1 to reduce the compute time
#   required for this test on CRAN;
#   We would not do this normally.
result <- smooth.monotone(age, hgt1, growfdPar, conv=0.1,
                          iterlim=1)
#  Extract the functional data object and regression
#  coefficients
Wfd  <- result$Wfdobj
beta <- result$beta
#  Evaluate the fitted height curve over a fine mesh
agefine <- seq(1,18,len=60)
hgtfine <- beta[1] + beta[2]*eval.monfd(agefine, Wfd)
#  Plot the data and the curve
plot(age, hgt1, type="p")
lines(agefine, hgtfine)
#  Evaluate the acceleration curve
accfine <- beta[2]*eval.monfd(agefine, Wfd, 2)
#  Plot the acceleration curve
plot(agefine, accfine, type="l")
lines(c(1,18),c(0,0),lty=4)

##
## using predict.monfd
##
hgtfit <- with(result, beta[1]+beta[2]*eval.monfd(argvals, Wfdobj))
hgtfit. <- fitted(result)

all.equal(hgtfit, hgtfit.)


accfine. <- predict(result, agefine, Lfdobj=2)

all.equal(accfine, accfine.)


growthResid <- resid(result)

all.equal(growthResid, with(result, y-hgtfit.))

par(oldpar)

Evaluate a Basis Penalty Matrix

Description

A basis roughness penalty matrix is the matrix containing the possible inner products of pairs of basis functions. These inner products are typically defined in terms of the value of a derivative or of a linear differential operator applied to the basis function. The basis penalty matrix plays an important role in the computation of functions whose roughness is controlled by a roughness penalty.

Usage

eval.penalty(basisobj, Lfdobj=int2Lfd(0), rng=rangeval)

Arguments

Details

The inner product can be computed exactly for many types of bases if $m$ is an integer. These include B-spline, fourier, exponential, monomial, polynomial and power bases. In other cases, and for noninteger operators, the inner products are computed by an iterative numerical integration method called Richard extrapolation using the trapezoidal rule.

If the penalty matrix must be evaluated repeatedly, computation can be greatly speeded up by avoiding the use of this function, and instead using quadrature points and weights defined by Simpson's rule.

Value

a square symmetric matrix whose order is equal to the number of basis functions defined by the basis function object basisobj . If Lfdobj is $m$ or a linear differential operator of order $m$, the rank of the matrix should be at least approximately equal to its order minus $m$.

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

getbasispenalty, eval.basis,

Evaluate a Positive Functional Data Object

Description

Evaluate a positive functional data object at specified argument values, or evaluate a derivative of the functional object.

Usage

eval.posfd(evalarg, Wfdobj, Lfdobj=int2Lfd(0))
## S3 method for class 'posfd'
predict(object, newdata=NULL, Lfdobj=0, ...)
## S3 method for class 'posfd'
fitted(object, ...)
## S3 method for class 'posfd'
residuals(object, ...)

Arguments

Details

A positive function data object $h(t)$ is defined by $h(t) =[exp Wfd](t)$. The function Wfdobj that defines the positive function is usually estimated by positive smoothing function smooth.pos

Value

a matrix containing the positive function values. The first dimension corresponds to the argument values in evalarg and the second to replications.

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

eval.fd, eval.monfd

Examples

harmaccelLfd <- vec2Lfd(c(0,(2*pi/365)^2,0), c(0, 365))
smallbasis   <- create.fourier.basis(c(0, 365), 65)
index        <- (1:35)[CanadianWeather$place == "Vancouver"]
VanPrec      <- CanadianWeather$dailyAv[,index, "Precipitation.mm"]
lambda       <- 1e4
dayfdPar     <- fdPar(fd(matrix(0,smallbasis$nbasis,1), smallbasis), 
                      harmaccelLfd, lambda)
VanPrecPos   <- smooth.pos(day.5, VanPrec, dayfdPar)
#  compute fitted values using eval.posfd()
VanPrecPosFit1 <- eval.posfd(day.5, VanPrecPos$Wfdobj)
#  compute fitted values using predict()
VanPrecPosFit2 <- predict(VanPrecPos, day.5)

all.equal(VanPrecPosFit1, VanPrecPosFit2)

#  compute fitted values using fitted()
VanPrecPosFit3 <- fitted(VanPrecPos)
#  compute residuals
VanPrecRes <- resid(VanPrecPos)

all.equal(VanPrecRes, VanPrecPos$y-VanPrecPosFit3)

Values of a Functional Data Object Defining Surprisal Curves.

Description

A surprisal vector of length M is minus the log to a positive integer base M of a set of M multinomial probabilities. Surprisal curves are functions of a one-dimensional index set, such that at any value of the index set the values of the curves are a surprisal vector. See Details below for further explanations.

Usage

eval.surp(evalarg, Wfdobj, nderiv = 0)

Arguments

Details

A surprisal M-vector is information measured in M-bits. Since a multinomial probability vector must sum to one, it follows that the surprisal vector S must satisfy the constraint log_M(sum(M^(-S)) = 0. That is, surprisal vectors lie within a curved M-1-dimensional manifold.

Surprisal curves are defined by a set of unconstrained M-1 B-spline functional data objects defined over an index set that are transformed into surprisal curves defined over the index set.

Let C be a K by M-1 coefficient matrix defining the B-spline curves, where K is the number of B-spline basis functions.

Let a M by M-1 matrix Z have orthonormal columns. Matrices satisfying these constraints are generated by function zerobasis().

Let N by K matrix be a matrix of B-spline basis values evaluated at N evaluation points using function eval.basis().

Let N by M matrix X = B * C * t(Z).

Then the N by M matrix S of surprisal values is S = -X + outer(log(rowSums(M^X))/log(M),rep(1,M)).

Value

A N by M matrix S of surprisal values at points evalarg, or their first or second derivatives.

Author(s)

Juan Li and James Ramsay

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

smooth.surp

Examples

#  see example in man/smooth.surp.Rd

Evaluate the Diagonal of a Bivariate Functional Data Object

Description

Bivariate function data objects are functions of two arguments, $f(s,t)$. It can be useful to evaluate the function for argument values satisfying $s=t$, such as evaluating the univariate variance function given the bivariate function that defines the variance-covariance function or surface. A linear differential operator can be applied to function $f(s,t)$ considered as a univariate function of either object holding the other object fixed.

Usage

evaldiag.bifd(evalarg, bifdobj, sLfd=int2Lfd(0), tLfd=int2Lfd(0))

Arguments

Value

a vector or matrix of diagonal function values.

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

var.fd, eval.bifd

Exponential Basis Function Values

Description

Evaluates a set of exponential basis functions, or a derivative of these functions, at a set of arguments.

Usage

expon(x, ratevec=1, nderiv=0)

Arguments

Details

There are no restrictions on the rate constants.

Value

a matrix of basis function values with rows corresponding to argument values and columns to basis functions.

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

exponpen

Powers of a functional data ('fd') object

Description

Exponentiate a functional data object where feasible.

Usage

## S3 method for class 'fd'
e1 ^ e2
exponentiate.fd(e1, e2, tolint=.Machine$double.eps^0.75,
  basisobj=e1$basis,
  tolfd=sqrt(.Machine$double.eps)*
          sqrt(sum(e1$coefs^2)+.Machine$double.eps)^abs(e2),
  maxbasis=NULL, npoints=NULL)

Arguments

Details

If e1 has a B-spline basis, this uses the B-spline algorithm.

Otherwise it throws an error unless it finds one of the following special cases:

e2 = 0

Return an fd object with a constant basis that is everywhere 1

e2 is a positive integer to within tolint

Multiply e1 by itself e2 times

e2 is positive and e1 has a Fourier basis

e120 <- e1^floor(e2)

outBasis <- e120$basis

rng <- outBasis$rangeval

Time <- seq(rng[1], rng[2], npoints)

e1.2 <- predict(e1, Time)^e2

fd1.2 <- smooth.basis(Time, e1.2, outBasis)$

d1.2 <- (e1.2 - predict(fd1.2, Time))

if(all(abs(d1.2)<tolfd))return(fd1.2)

Else if(outBasis$nbasis<maxbasis) increase the size of outBasis and try again.

Else write a warning with the max(abs(d1.2)) and return fd1.2.

Value

A function data object approximating the desired power.

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

arithmetic.fd basisfd, basisfd.product

Examples

##
## sin^2
##

basis3 <- create.fourier.basis(nbasis=3)
oldpar <- par(no.readonly=TRUE)
plot(basis3)
# max = sqrt(2), so
# integral of the square of each basis function (from 0 to 1) is 1
integrate(function(x)sin(2*pi*x)^2, 0, 1) # = 0.5

# sin(theta)
fdsin <- fd(c(0,sqrt(0.5),0), basis3)
plot(fdsin)

fdsin2 <- fdsin^2

# check
fdsinsin <- fdsin*fdsin
# sin^2(pi*time) = 0.5*(1-cos(2*pi*theta) basic trig identity
plot(fdsinsin) # good


all.equal(fdsin2, fdsinsin)

par(oldpar)

Exponential Penalty Matrix

Description

Computes the matrix defining the roughness penalty for functions expressed in terms of an exponential basis.

Usage

exponpen(basisobj, Lfdobj=int2Lfd(2))

Arguments

Details

A roughness penalty for a function $x(t)$ is defined by integrating the square of either the derivative of $ x(t) $ or, more generally, the result of applying a linear differential operator $L$ to it. The most common roughness penalty is the integral of the square of the second derivative, and this is the default. To apply this roughness penalty, the matrix of inner products of the basis functions (possibly after applying the linear differential operator to them) defining this function is necessary. This function just calls the roughness penalty evaluation function specific to the basis involved.

Value

a symmetric matrix of order equal to the number of basis functions defined by the exponential basis object. Each element is the inner product of two exponential basis functions after applying the derivative or linear differential operator defined by Lfdobj.

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

expon, eval.penalty, getbasispenalty

Examples

#  set up an exponential basis with 3 basis functions
ratevec  <- c(0, -1, -5)
basisobj <- create.exponential.basis(c(0,1),3,ratevec)
#  compute the 3 by 3 matrix of inner products of
#  second derivatives
penmat <- exponpen(basisobj)

Functional Boxplots

Description

Produces functional boxplots or enhanced functional boxplots of the given functional data. It can also be used to carry out functional data ordering based on band depth.

Usage

fbplot(fit, x = NULL, method = "MBD", depth = NULL, plot = TRUE,
	 prob = 0.5, color = 6, outliercol = 2, barcol = 4,
	 fullout=FALSE, factor=1.5,xlim=c(1,nrow(fit)),
	 ylim=c(min(fit)-.5*diff(range(fit)),max(fit)+.5*diff(range(fit))),...)
## S3 method for class 'fd'
boxplot(x, z=NULL, ...)
## S3 method for class 'fdPar'
boxplot(x, z=NULL, ...)
## S3 method for class 'fdSmooth'
boxplot(x, z=NULL, ...)

Arguments

Details

For functional data, the band depth (BD) or modified band depth (MBD) allows for ordering a sample of curves from the center outwards and, thus, introduces a measure to define functional quantiles and the centrality or outlyingness of an observation. A smaller rank is associated with a more central position with respect to the sample curves. BD usually provides many ties (curves have the same depth values), but MBD does not. "BD2" uses two curves to determine a band. The method "Both" uses "BD2" first and then uses "MBD" to break ties. The method "Both" uses BD2 first and then uses MBD to break ties. The computation is carried out by the fast algorithm proposed by Sun et. al. (2012).

Value

Author(s)

Ying Sun sunwards@stat.osu.edu

Marc G. Genton marc.genton@kaust.edu.sa

References

Sun, Y., Genton, M. G. and Nychka, D. (2012), "Exact fast computation of band depth for large functional datasets: How quickly can one million curves be ranked?" Stat, 1, 68-74.

Sun, Y. and Genton, M. G. (2011), "Functional Boxplots," Journal of Computational and Graphical Statistics, 20, 316-334.

Lopez-Pintado, S. and Romo, J. (2009), "On the concept of depth for functional data," Journal of the American Statistical Association, 104, 718-734.

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

Examples

##
## 1.  generate 50 random curves with some covariance structure
##     model 1 without outliers
##
cov.fun=function(d,k,c,mu){
        k*exp(-c*d^mu)
}
n=50
p=30
t=seq(0,1,len=p)
d=dist(t,upper=TRUE,diag=TRUE)
d.matrix=as.matrix(d)
#covariance function in time
t.cov=cov.fun(d.matrix,1,1,1)
# Cholesky Decomposition
L=chol(t.cov)
mu=4*t
e=matrix(rnorm(n*p),p,n)
ydata  = mu+t(L)%*%e

#functional boxplot
oldpar <- par(no.readonly=TRUE)
fbplot(ydata,method='MBD',ylim=c(-11,15))

# The same using boxplot.fd
boxplot.fd(ydata, method='MBD', ylim=c(-11, 15))

# same with default ylim
boxplot.fd(ydata)

##
## 2.  as an fd object
##
T      = dim(ydata)[1]
time   = seq(0,T,len=T)
ybasis = create.bspline.basis(c(0,T), 23)
Yfd    = smooth.basis(time, ydata, ybasis)$fd
boxplot(Yfd)

##
## 3.  as an fdPar object
##
Ypar <- fdPar(Yfd)
boxplot(Ypar)

##
## 4.  Smoothed version
##
Ysmooth <- smooth.fdPar(Yfd)
boxplot(Ysmooth)

##
## 5.  model 2 with outliers
##
#magnitude
k=6
#randomly introduce outliers
C=rbinom(n,1,0.1)
s=2*rbinom(n,1,0.5)-1
cs.m=matrix(C*s,p,n,byrow=TRUE)

e=matrix(rnorm(n*p),p,n)
y=mu+t(L)%*%e+k*cs.m

#functional boxplot
fbplot(y,method='MBD',ylim=c(-11,15))
par(oldpar)

Define a Functional Data Object

Description

This is the constructor function for objects of the fd class. Each function that sets up an object of this class must call this function. This includes functions smooth.basis, density.fd, and so forth that estimate functional data objects that smooth or otherwise represent data. Ordinarily, users of the functional data analysis software will not need to call this function directly, but these notes are valuable to understanding the components of a list of class fd.

Usage

fd(coef=NULL, basisobj=NULL, fdnames=NULL)

Arguments

Details

To check that an object is of this class, use function is.fd.

Normally only developers of new functional data analysis functions will actually need to use this function.

Value

A functional data object (i.e., having class fd), which is a list with components named coefs, basis, and fdnames.

Source

Ramsay, James O., and Silverman, Bernard W. (2006), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York

See Also

smooth.basis smooth.fdPar smooth.basisPar density.fd create.bspline.basis arithmetic.fd

Examples

##
## default
##
fd()
oldpar <- par(no.readonly=TRUE)
##
## The simplest b-spline basis:  order 1, degree 0, zero interior knots:
##       a single step function
##
bspl1.1    <- create.bspline.basis(norder=1, breaks=0:1)
fd.bspl1.1 <- fd(0, basisobj=bspl1.1)

fd.bspl1.1a <- fd(basisobj=bspl1.1)

all.equal(fd.bspl1.1, fd.bspl1.1a)

# TRUE

# the following three lines shown an error in a non-cran check:
# if(!CRAN()) {
#   fd.bspl1.1b <- fd(0)
# }

##
## Cubic spline:  4  basis functions
##
bspl4 <- create.bspline.basis(nbasis=4)
plot(bspl4)
parab4.5 <- fd(c(3, -1, -1, 3)/3, bspl4)
# = 4*(x-.5)^2
plot(parab4.5)

##
## Fourier basis
##
f3 <- fd(c(0,0,1), create.fourier.basis())
plot(f3)
# range over +/-sqrt(2), because
# integral from 0 to 1 of cos^2 = 1/2
# so multiply by sqrt(2) to get
# its square to integrate to 1.

##
## subset of an fd object
##

gaitbasis3 <- create.fourier.basis(nbasis=5)
gaittime = (1:20)/21
gaitfd3    <- smooth.basis(gaittime, gait, gaitbasis3)$fd
gaitfd3[1]
par(oldpar)

Convert a univariate functional data object to a list

Description

Convert a univariate functional data object to a list for input to Lfd.

Usage

fd2list(fdobj)

Arguments

Value

a list as required for the second argument of Lfd.

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

Lfd

Examples

Lbasis  = create.constant.basis(c(0,365));  #  create a constant basis
Lcoef   = matrix(c(0,(2*pi/365)^2,0),1,3)   #  set up three coefficients
wfdobj  = fd(Lcoef,Lbasis)      # define an FD object for weight functions
wfdlist = fd2list(wfdobj)       # convert the FD object to a cell object
harmaccelLfd = Lfd(3, wfdlist)  #  define the operator object

Functions for statistical analyses of functions

Description

The data analyses that we use for data matrices or dataFrames, such as means, covariances, linear regressions, principal and canonical component analyses, can also be applied to samples of functions or curves. This package provides functional versions of these analyses along with plotting and other assessment tools.

The functions themselves are often a consequence of smoothing discrete data values over domains like time, space, and other continuous variables. Data smoothing tools are also provided.

But there are transformations of functions that have no meaning for data matrices for which rows may be re-ordered. Derivatives and integrals are often used to set up dynamic models, and can also play a constructive role in the data smoothing process. Methods for these and other functional operations are also available in this package.

There are now many texts and papers on functional data analysis. The two resources provided by the buildeers of the fda package are:

James Ramsay and Bernard Silverman (2005) Functional Data Analysis. New York, Springer.

James Ramsay, Giles Hooker and Spencer Graves (2009) Functional Data Analysis with R and Matlab. New York: Springer.

Another relevant package is the package Data2LD that offers ann introduction and methods for constructing linear and nonlinear differential equation models, along with the text:

James Ramsay and Giles Hooker, Dynamic Data Analysis, New York: Springer.

Extract plot labels and names for replicates and variables

Description

Extract plot labels and, if available, names for each replicate and variable

Usage

fdlabels(fdnames, nrep, nvar) 

Arguments

Details

xlabel <- if(length(fdnames[[1]])>1) names(fdnames)[1] else fdnames[[1]]

ylabel <- if(length(fdnames[[3]])>1) names(fdnames)[3] else fdnames[[3]]

casenames <- if(length(fdnames[[2]])== nrep)fdnames[[2]] else NULL

varnames <- if(length(fdnames[[3]])==nvar)fdnames[[3]] else NULL

Value

A list of xlabel, ylabel, casenames, and varnames

Author(s)

Jim Ramsay

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

plot.fd

Define a Functional Parameter Object

Description

Functional parameter objects are used as arguments to functions that estimate functional parameters, such as smoothing functions like smooth.basis. A functional parameter object is a functional data object with additional slots specifying a roughness penalty, a smoothing parameter and whether or not the functional parameter is to be estimated or held fixed. Functional parameter objects are used as arguments to functions that estimate functional parameters.

Usage

  fdPar(fdobj=NULL, Lfdobj=NULL, lambda=0, estimate=TRUE, penmat=NULL)

Arguments

Details

Functional parameters are often needed to specify initial values for iteratively refined estimates, as is the case in functions register.fd and smooth.monotone.

Often a list of functional parameters must be supplied to a function as an argument, and it may be that some of these parameters are considered known and must remain fixed during the analysis. This is the case for functions fRegress and pda.fd, for example.

Value

a functional parameter object (i.e., an object of class fdPar), which is a list with the following components:

Source

Ramsay, James O., and Silverman, Bernard W. (2006), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York

See Also

cca.fd, density.fd, fRegress, intensity.fd, pca.fd, smooth.fdPar, smooth.basis, smooth.monotone, int2Lfd

Examples

oldpar <- par(no.readonly=TRUE)
##
## Simple example
##
#  set up range for density
rangeval <- c(-3,3)
#  set up some standard normal data
x <- rnorm(50)
#  make sure values within the range
x[x < -3] <- -2.99
x[x >  3] <-  2.99
#  set up basis for W(x)
basisobj <- create.bspline.basis(rangeval, 11)
#  set up initial value for Wfdobj
Wfd0 <- fd(matrix(0,11,1), basisobj)
WfdParobj <- fdPar(Wfd0)

WfdP3 <- fdPar(seq(-3, 3, length=11))

##
##  smooth the Canadian daily temperature data
##
#    set up the fourier basis
nbasis   <- 365
dayrange <- c(0,365)
daybasis <- create.fourier.basis(dayrange, nbasis)
dayperiod <- 365
harmaccelLfd <- vec2Lfd(c(0,(2*pi/365)^2,0), dayrange)
#  Make temperature fd object
#  Temperature data are in 12 by 365 matrix tempav
#    See analyses of weather data.
#  Set up sampling points at mid days
daytime  <- (1:365)-0.5
#  Convert the data to a functional data object
daybasis65 <- create.fourier.basis(dayrange, nbasis, dayperiod)
templambda <- 1e1
tempfdPar  <- fdPar(fdobj=daybasis65, Lfdobj=harmaccelLfd,
                    lambda=templambda)

#FIXME
#tempfd <- smooth.basis(CanadianWeather$tempav, daytime, tempfdPar)$fd
#  Set up the harmonic acceleration operator
Lbasis  <- create.constant.basis(dayrange);
Lcoef   <- matrix(c(0,(2*pi/365)^2,0),1,3)
bfdobj  <- fd(Lcoef,Lbasis)
bwtlist <- fd2list(bfdobj)
harmaccelLfd <- Lfd(3, bwtlist)
#  Define the functional parameter object for
#  smoothing the temperature data
lambda   <- 0.01  #  minimum GCV estimate
#tempPar <- fdPar(daybasis365, harmaccelLfd, lambda)
#  smooth the data
#tempfd <- smooth.basis(daytime, CanadialWeather$tempav, tempPar)$fd
#  plot the temperature curves
#plot(tempfd)

##
## with rangeval of class Date and POSIXct
##
par(oldpar)

Convert fd or basisfd Objects to fdPar Objects

Description

If the input is an fd object, default fdPar parameters are added to convert the object to an fdPar object. If a basisfd object, it is first converted to an fd object with a nbasis by ncurve coefficient matrix of zeros, and then converted to a fdPar object.

Usage

fdParcheck(fdParobj, ncurve=NULL)

Arguments

Details

The previous version of this function did not include the argument ncurve, and set up an fd object with a square matrix of zeros. This could cause an error in many of the functions that called it.

Value

Returns a fdPar object.

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

Fourier Basis Function Values

Description

Evaluates a set of Fourier basis functions, or a derivative of these functions, at a set of arguments.

Usage

fourier(x, nbasis=n, period=span, nderiv=0)

Arguments

Value

a matrix of function values. The number of rows equals the number of arguments, and the number of columns equals the number of basis functions.

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

fourierpen

Examples

#  set up a set of 11 argument values
x <- seq(0,1,0.1)
names(x) <- paste("x", 0:10, sep="")
#  compute values for five Fourier basis functions
#  with the default period (1) and derivative (0)
(basismat <- fourier(x, 5))

# Create a false Fourier basis, i.e., nbasis = 1
# = a constant function
fourier(x, 1)

Fourier Penalty Matrix

Description

Computes the matrix defining the roughness penalty for functions expressed in terms of a Fourier basis.

Usage

fourierpen(basisobj, Lfdobj=int2Lfd(2))

Arguments

Details

A roughness penalty for a function $x(t)$ is defined by integrating the square of either the derivative of $ x(t) $ or, more generally, the result of applying a linear differential operator $L$ to it. The most common roughness penalty is the integral of the square of the second derivative, and this is the default. To apply this roughness penalty, the matrix of inner products of the basis functions (possibly after applying the linear differential operator to them) defining this function is necessary. This function just calls the roughness penalty evaluation function specific to the basis involved.

Value

a symmetric matrix of order equal to the number of basis functions defined by the Fourier basis object. Each element is the inner product of two Fourier basis functions after applying the derivative or linear differential operator defined by Lfdobj.

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

fourier, eval.penalty, getbasispenalty

Examples

#  set up a Fourier basis with 13 basis functions
#  and and period 1.0.
basisobj <- create.fourier.basis(c(0,1),13)
#  compute the 13 by 13 matrix of inner products
#  of second derivatives
penmat <- fourierpen(basisobj)

Permutation F-test for functional linear regression.

Description

Fperm.fd creates a null distribution for a test of no effect in functional linear regression. It makes generic use of fRegress and permutes the yfdPar input.

Usage

Fperm.fd(yfdPar, xfdlist, betalist, wt=NULL, nperm=200,
         argvals=NULL, q=0.05, plotres=TRUE, ...)

Arguments

Details

An F-statistic is calculated as the ratio of residual variance to predicted variance. The observed F-statistic is returned along with the permutation distribution.

If yfdPar is a fd object, the maximal value of the pointwise F-statistic is calculated. The pointwise F-statistics are also returned.

The default of setting q = 0.95 is, by now, fairly standard. The default nperm = 200 may be small, depending on the amount of computing time available.

If argvals is not specified and yfdPar is a fd object, it defaults to 101 equally-spaced points on the range of yfdPar.

Value

A list with the following components:

Side Effects

a plot of the functional observations

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

fRegress, Fstat.fd

Examples

oldpar <- par(no.readonly=TRUE)
##
## 1.  yfdPar = vector
##
annualprec <- log10(apply(
    CanadianWeather$dailyAv[,,"Precipitation.mm"], 2,sum))

#  set up a smaller basis using only 40 Fourier basis functions
#  to save some computation time

smallnbasis <- 40
smallbasis  <- create.fourier.basis(c(0, 365), smallnbasis)
tempfd      <- smooth.basis(day.5, CanadianWeather$dailyAv[,,"Temperature.C"],
                            smallbasis)$fd
constantfd <- fd(matrix(1,1,35), create.constant.basis(c(0, 365)))

xfdlist <- vector("list",2)
xfdlist[[1]] <- constantfd
xfdlist[[2]] <- tempfd[1:35]

betalist   <- vector("list",2)
#  set up the first regression function as a constant
betabasis1 <- create.constant.basis(c(0, 365))
betafd1    <- fd(0, betabasis1)
betafdPar1 <- fdPar(betafd1)
betalist[[1]] <- betafdPar1

nbetabasis  <- 35
betabasis2  <- create.fourier.basis(c(0, 365), nbetabasis)
betafd2     <- fd(matrix(0,nbetabasis,1), betabasis2)

lambda          <- 10^12.5
harmaccelLfd365 <- vec2Lfd(c(0,(2*pi/365)^2,0), c(0, 365))
betafdPar2      <- fdPar(betafd2, harmaccelLfd365, lambda)
betalist[[2]]   <- betafdPar2

# Should use the default nperm = 200
# but use 10 to save test time for illustration

F.res2 = Fperm.fd(annualprec, xfdlist, betalist, nperm=100)

##
## 2.  yfdpar = Functional data object (class fd)
##
# The very simplest example is the equivalent of the permutation
# t-test on the growth data.

# First set up a basis system to hold the smooths

if (!CRAN()) {

Knots  <- growth$age
norder <- 6
nbasis <- length(Knots) + norder - 2
hgtbasis <- create.bspline.basis(range(Knots), nbasis, norder, Knots)

# Now smooth with a fourth-derivative penalty and a very small smoothing
# parameter

Lfdobj    <- 4
lambda    <- 1e-2
growfd    <- fd(matrix(0,nbasis,1),hgtbasis)
growfdPar <- fdPar(growfd, Lfdobj, lambda)

hgtfd <- smooth.basis(growth$age,
                      cbind(growth$hgtm,growth$hgtf),growfdPar)$fd

# Now set up factors for fRegress:

cbasis = create.constant.basis(range(Knots))

maleind = c(rep(1,ncol(growth$hgtm)),rep(0,ncol(growth$hgtf)))

constfd = fd( matrix(1,1,length(maleind)),cbasis)
maleindfd = fd( matrix(maleind,1,length(maleind)),cbasis)

xfdlist = list(constfd,maleindfd)

# The fdPar object for the coefficients and call Fperm.fd

betalist = list(fdPar(hgtfd,2,1e-6),fdPar(hgtfd,2,1e-6))

# Should use nperm = 200 or so,
# but use 10 to save test time

Fres = Fperm.fd(hgtfd,xfdlist,betalist,nperm=100)

par(oldpar)
}

Functional Regression Analysis

Description

This function carries out a functional regression analysis, where either the dependent variable or one or more independent variables are functional. Non-functional variables may be used on either side of the equation. In a simple problem where there is a single scalar independent covariate with values z_i, i=1,\ldots,N and a single functional covariate with values x_i(t), the two versions of the model fit by fRegress are the scalar dependent variable model

y_i = \beta_1 z_i + \int x_i(t) \beta_2(t) \, dt + e_i

and the concurrent functional dependent variable model

y_i(t) = \beta_1(t) z_i + \beta_2(t) x_i(t) + e_i(t).

In these models, the final term e_i or e_i(t) is a residual, lack of fit or error term.

In the concurrent functional linear model for a functional dependent variable, all functional variables are all evaluated at a common time or argument value t. That is, the fit is defined in terms of the behavior of all variables at a fixed time, or in terms of "now" behavior.

All regression coefficient functions \beta_j(t) are considered to be functional. In the case of a scalar dependent variable, the regression coefficient for a scalar covariate is converted to a functional variable with a constant basis. All regression coefficient functions can be forced to be smooth through the use of roughness penalties, and consequently are specified in the argument list as functional parameter objects.

Usage

fRegress(y, ...)
## S3 method for class 'fd'
fRegress(y, xfdlist, betalist, wt=NULL,
                     y2cMap=NULL, SigmaE=NULL, returnMatrix=FALSE, 
                        method=c('fRegress', 'model'), sep='.', ...)
## S3 method for class 'double'
fRegress(y, xfdlist, betalist, wt=NULL,
                     y2cMap=NULL, SigmaE=NULL, returnMatrix=FALSE, ...)
## S3 method for class 'formula'
fRegress(y, data=NULL, betalist=NULL, wt=NULL,
                 y2cMap=NULL, SigmaE=NULL,
                 method='fRegress', sep='.', ...)
## S3 method for class 'character'
fRegress(y, data=NULL, betalist=NULL, wt=NULL,
                 y2cMap=NULL, SigmaE=NULL,
                 method='fRegress', sep='.', ...)

Arguments

Details

Alternative forms of functional regression can be categorized with traditional least squares using the following 2 x 2 table:

For fRegress.numeric, the numeric response is assumed to be the sum of integrals of xfd * beta for all functional xfd terms.

fRegress.fd or .fdPar produces a concurrent regression with each beta being also a (univariate) function.

linmod predicts a functional response from a convolution integral, estimating a bivariate regression function.

In the computation of regression function estimates in fRegress, all independent variables are treated as if they are functional. If argument xfdlist contains one or more vectors, these are converted to functional data objects having the constant basis with coefficients equal to the elements of the vector.

Needless to say, if all the variables in the model are scalar, do NOT use this function. Instead, use either lm or lsfit.

These functions provide a partial implementation of Ramsay and Silverman (2005, chapters 12-20).

Value

These functions return either a standard fRegress fit object or or a model specification:

Author(s)

J. O. Ramsay, Giles Hooker, and Spencer Graves

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

fRegress.stderr, fRegress.CV, Fperm.fd, Fstat.fd, linmod

Examples

oldpar <- par(no.readonly=TRUE)
###
###
###   vector response with functional explanatory variable  
###
###

#  data are in Canadian Weather object
#  print the names of the data
print(names(CanadianWeather))
#  set up log10 of annual precipitation for 35 weather stations
annualprec <- 
    log10(apply(CanadianWeather$dailyAv[,,"Precipitation.mm"], 2,sum))
# The simplest 'fRegress' call is singular with more bases
# than observations, so we use only 25 basis functions, for this example
smallbasis  <- create.fourier.basis(c(0, 365), 25)
# The covariate is the temperature curve for each station.
tempfd <- 
 smooth.basis(day.5, CanadianWeather$dailyAv[,,"Temperature.C"], smallbasis)$fd
##
## formula interface:  specify the model by a formula, the method
## fRegress.formula automatically sets up the regression coefficient functions,
## a constant function for the intercept, 
## and a higher dimensional function
## for the inner product with temperature
##

precip.Temp1 <- fRegress(annualprec ~ tempfd, method="fRegress")

#  the output is a list with class name fRegress, display names
names(precip.Temp1)
#[c1] "yvec"           "xfdlist"        "betalist"       "betaestlist"    "yhatfdobj"     
# [6] "Cmat"           "Dmat"           "Cmatinv"        "wt"             "df"            
#[11] "GCV"            "OCV"            "y2cMap"         "SigmaE"         "betastderrlist"
#[16] "bvar"           "c2bMap"       

#  the vector of fits to the data is object  precip.Temp1$yfdPar,
#  but since the dependent variable is a vector, so is the fit
annualprec.fit1 <- precip.Temp1$yhatfdobj
#  plot the data and the fit
plot(annualprec.fit1, annualprec, type="p", pch="o")
lines(annualprec.fit1, annualprec.fit1, lty=2)
#  print root mean squared error
RMSE <- round(sqrt(mean((annualprec-annualprec.fit1)^2)),3)
print(paste("RMSE =",RMSE))
#  plot the estimated regression function
plot(precip.Temp1$betaestlist[[2]])
#  This isn't helpful either, the coefficient function is too
#  complicated to interpret.
#  display the number of basis functions used:
print(precip.Temp1$betaestlist[[2]]$fd$basis$nbasis)
#  25 basis functions to fit 35 values, no wonder we over-fit the data

##
## Get the default setup and modify it
## the "model" value of the method argument causes the analysis
## to produce a list vector of arguments for calling the
## fRegress function
##

precip.Temp.mdl1 <- fRegress(annualprec ~ tempfd, method="model")
# First confirm we get the same answer as above by calling
# function fRegress() with these arguments:
precip.Temp.m <- do.call('fRegress', precip.Temp.mdl1)

all.equal(precip.Temp.m, precip.Temp1)


#  set up a smaller basis for beta2 than for temperature so that we
#  get a more parsimonious fit to the data

nbetabasis2 <- 21  #  not much less, but we add some roughness penalization
betabasis2  <- create.fourier.basis(c(0, 365), nbetabasis2)
betafd2     <- fd(rep(0, nbetabasis2), betabasis2)
# add smoothing
betafdPar2  <- fdPar(betafd2, lambda=10)

# replace the regress coefficient function with this fdPar object

precip.Temp.mdl2 <- precip.Temp.mdl1
precip.Temp.mdl2[['betalist']][['tempfd']] <- betafdPar2

# Now do re-fit the data

precip.Temp2 <- do.call('fRegress', precip.Temp.mdl2)

# Compare the two fits:
#  degrees of freedom
precip.Temp1[['df']] # 26
precip.Temp2[['df']] # 22
#  root-mean-squared errors:
RMSE1 <- round(sqrt(mean(with(precip.Temp1, (yhatfdobj-yvec)^2))),3)
RMSE2 <- round(sqrt(mean(with(precip.Temp2, (yhatfdobj-yvec)^2))),3)
print(c(RMSE1, RMSE2))
#  display further results for the more parsimonious model
annualprec.fit2 <- precip.Temp2$yhatfdobj
plot(annualprec.fit2, annualprec, type="p", pch="o")
lines(annualprec.fit2, annualprec.fit2, lty=2)
#  plot the estimated regression function
plot(precip.Temp2$betaestlist[[2]])
#  now we see that it is primarily the temperatures in the
#  early winter that provide the fit to log precipitation by temperature

##
## Manual construction of xfdlist and betalist
##

xfdlist <- list(const=rep(1, 35), tempfd=tempfd)

# The intercept must be constant for a scalar response
betabasis1 <- create.constant.basis(c(0, 365))
betafd1    <- fd(0, betabasis1)
betafdPar1 <- fdPar(betafd1)

betafd2     <- fd(matrix(0,7,1), create.bspline.basis(c(0, 365),7))
# convert to an fdPar object
betafdPar2  <- fdPar(betafd2)

betalist <- list(const=betafdPar1, tempfd=betafdPar2)

precip.Temp3   <- fRegress(annualprec, xfdlist, betalist)
annualprec.fit3 <- precip.Temp3$yhatfdobj
#  plot the data and the fit
plot(annualprec.fit3, annualprec, type="p", pch="o")
lines(annualprec.fit3, annualprec.fit3)
plot(precip.Temp3$betaestlist[[2]])

###
###
###  functional response with vector explanatory variables  
###
###

##
## simplest:  formula interface
##

daybasis65 <- create.fourier.basis(rangeval=c(0, 365), nbasis=65,
                  axes=list('axesIntervals'))
Temp.fd <- with(CanadianWeather, smooth.basisPar(day.5,
                dailyAv[,,'Temperature.C'], daybasis65)$fd)
TempRgn.f <- fRegress(Temp.fd ~ region, CanadianWeather)

##
## Get the default setup and possibly modify it
##

TempRgn.mdl <- fRegress(Temp.fd ~ region, CanadianWeather, method='model')

# make desired modifications here
# then run

TempRgn.m <- do.call('fRegress', TempRgn.mdl)

# no change, so match the first run

all.equal(TempRgn.m, TempRgn.f)


##
## More detailed set up
##

region.contrasts <- model.matrix(~factor(CanadianWeather$region))
rgnContr3 <- region.contrasts
dim(rgnContr3) <- c(1, 35, 4)
dimnames(rgnContr3) <- list('', CanadianWeather$place, c('const',
   paste('region', c('Atlantic', 'Continental', 'Pacific'), sep='.')) )

const365 <- create.constant.basis(c(0, 365))
region.fd.Atlantic <- fd(matrix(rgnContr3[,,2], 1), const365)
# str(region.fd.Atlantic)
region.fd.Continental <- fd(matrix(rgnContr3[,,3], 1), const365)
region.fd.Pacific <- fd(matrix(rgnContr3[,,4], 1), const365)
region.fdlist <- list(const=rep(1, 35),
     region.Atlantic=region.fd.Atlantic,
     region.Continental=region.fd.Continental,
     region.Pacific=region.fd.Pacific)
# str(TempRgn.mdl$betalist)

###
###
###  functional response with functional explanatory variable  
###
###

##
##  predict knee angle from hip angle;  
##     from demo('gait', package='fda')

##
## formula interface
##
gaittime   <- as.matrix((1:20)/21)
gaitrange  <- c(0,20)
gaitbasis  <- create.fourier.basis(gaitrange, nbasis=21)
gaitnbasis <- gaitbasis$nbasis
gaitcoef   <- matrix(0,gaitnbasis,dim(gait)[2])
harmaccelLfd <- vec2Lfd(c(0, (2*pi/20)^2, 0), rangeval=gaitrange)
gaitfd     <- smooth.basisPar(gaittime, gait, gaitbasis, 
                          Lfdobj=harmaccelLfd, lambda=1e-2)$fd
hipfd  <- gaitfd[,1]
kneefd <- gaitfd[,2]

knee.hip.f <- fRegress(kneefd ~ hipfd)

##
## manual set-up
##

#  set up the list of covariate objects
const  <- rep(1, dim(kneefd$coef)[2])
xfdlist  <- list(const=const, hipfd=hipfd)

beta0 <- with(kneefd, fd(gaitcoef, gaitbasis, fdnames))
beta1 <- with(hipfd,  fd(gaitcoef, gaitbasis, fdnames))

betalist  <- list(const=fdPar(beta0), hipfd=fdPar(beta1))

fRegressout <- fRegress(kneefd, xfdlist, betalist)
par(oldpar)

Computes Cross-validated Error Sum of Integrated Squared Errors for a Functional Regression Model

Description

For a functional regression model, a cross-validated error sum of squares is computed. For a functional dependent variable this is the sum of integrated squared errors. For a scalar response, this function has been superseded by the OCV and gcv elements returned by fRegress. This function aids the choice of smoothing parameters in this model using the cross-validated error sum of squares criterion.

Usage

#fRegress.CV(y, xfdlist, betalist, wt=NULL, CVobs=1:N,
#            returnMatrix=FALSE, ...)

#NOTE:  The following is required by CRAN rules that
# function names like "as.numeric" must follow the documentation
# standards for S3 generics, even when they are not.
# Please ignore the following line:
## S3 method for class 'CV'
fRegress(y, xfdlist, betalist, wt=NULL, CVobs=1:N,
            returnMatrix=FALSE, ...)

Arguments

Value

A list containing

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

fRegress, fRegress.stderr

Examples

  #. See the analyses of the Canadian daily weather data.

Compute Standard errors of Coefficient Functions Estimated by Functional Regression Analysis

Description

Function fRegress carries out a functional regression analysis of the concurrent kind, and estimates a regression coefficient function corresponding to each independent variable, whether it is scalar or functional. This function uses the list that is output by fRegress to provide standard error functions for each regression function. These standard error functions are pointwise, meaning that sampling standard deviation functions only are computed, and not sampling covariances.

Usage

## S3 method for class 'stderr'
fRegress(y, y2cMap, SigmaE, returnMatrix=FALSE, ...)

Arguments

Value

a named list of length 3 containing:

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

fRegress, fRegress.CV

Examples

#See the weather data analyses in the file daily.ssc for
#examples of the use of function fRegress.stderr.

F-statistic for functional linear regression.

Description

Fstat.fd calculates a pointwise F-statistic for functional linear regression.

Usage

 Fstat.fd(y,yhat,argvals=NULL)

Arguments

Details

An F-statistic is calculated as the ratio of residual variance to predicted variance.

If argvals is not specified and yfdPar is a fd object, it defaults to 101 equally-spaced points on the range of yfdPar.

Value

A list with components

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

fRegress Fstat.fd

Hip and knee angle while walking

Description

Hip and knee angle in degrees through a 20 point movement cycle for 39 boys

Format

An array of dim c(20, 39, 2) giving the "Hip Angle" and "Knee Angle" for 39 repetitions of a 20 point gait cycle.

Details

The components of dimnames(gait) are as follows:

[[1]] standardized gait time = seq(from=0.025, to=0.975, by=0.05)

[[2]] subject ID = "boy1", "boy2", ..., "boy39"

[[3]] gait variable = "Hip Angle" or "Knee Angle"

Examples

oldpar <- par(no.readonly=TRUE)
plot(gait[,1, 1], gait[, 1, 2], type="b")
par(oldpar)

Generalized eigenanalysis

Description

Find matrices L and M to maximize

tr(L'AM) / sqrt(tr(L'BL) tr(M'CM'))

where A = a p x q matrix, B = p x p symmetric, positive definite matrix, B = q x q symmetric positive definite matrix, L = p x s matrix, and M = q x s matrix, where s = the number of non-zero generalized eigenvalues of A.

Usage

geigen(Amat, Bmat, Cmat)

Arguments

Value

list(values, Lmat, Mmat)

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

eigen

Examples

A <- matrix(1:6, 2)
B <- matrix(c(2, 1, 1, 2), 2)
C <- diag(1:3)
ABC <- geigen(A, B, C)

Values of Basis Functions or their Derivatives

Description

Evaluate a set of basis functions or their derivatives at a set of argument values.

Usage

getbasismatrix(evalarg, basisobj, nderiv=0, returnMatrix=FALSE)

Arguments

Value

a matrix of basis function or derivative values. Rows correspond to argument values and columns to basis functions.

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

eval.fd

Examples

##
## Minimal example:  a B-spline of order 1, i.e., a step function
## with 0 interior knots:
##
bspl1.1 <- create.bspline.basis(norder=1, breaks=0:1)
m <- getbasismatrix(seq(0, 1, .2), bspl1.1)

# check
m. <- matrix(rep(1, 6), 6,
    dimnames=list(NULL, 'bspl') )

all.equal(m, m.)


##
## Date and POSIXct
##
# Date
July4.1776    <- as.Date('1776-07-04')
Apr30.1789    <- as.Date('1789-04-30')
AmRev         <- c(July4.1776, Apr30.1789)
BspRevolution <- create.bspline.basis(AmRev)
AmRevYears    <- as.numeric(seq(July4.1776, Apr30.1789, length.out=14))
AmRevMatrix   <- getbasismatrix(AmRevYears, BspRevolution)
matplot(AmRevYears, AmRevMatrix, type='b')

# POSIXct
AmRev.ct    <- as.POSIXct1970(c('1776-07-04', '1789-04-30'))
BspRev.ct   <- create.bspline.basis(AmRev.ct)
AmRevYrs.ct <- as.numeric(seq(AmRev.ct[1], AmRev.ct[2], length.out=14))
AmRevMat.ct <- getbasismatrix(AmRevYrs.ct, BspRev.ct)
matplot(AmRevYrs.ct, AmRevMat.ct, type='b')

Evaluate a Roughness Penalty Matrix

Description

A basis roughness penalty matrix is the matrix containing the possible inner products of pairs of basis functions. These inner products are typically defined in terms of the value of a derivative or of a linear differential operator applied to the basis function. The basis penalty matrix plays an important role in the computation of functions whose roughness is controlled by a roughness penalty.

Usage

getbasispenalty(basisobj, Lfdobj=NULL)

Arguments

Details

A roughness penalty for a function $x(t)$ is defined by integrating the square of either the derivative of $ x(t) $ or, more generally, the result of applying a linear differential operator $L$ to it. The most common roughness penalty is the integral of the square of the second derivative, and this is the default. To apply this roughness penalty, the matrix of inner products of the basis functions defining this function is necessary. This function just calls the roughness penalty evaluation function specific to the basis involved.

Value

a symmetric matrix of order equal to the number of basis functions defined by the B-spline basis object. Each element is the inner product of two B-spline basis functions after taking the derivative.

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

eval.penalty

Examples

#  set up a B-spline basis of order 4 with 13 basis functions
#  and knots at 0.0, 0.1,..., 0.9, 1.0.
basisobj <- create.bspline.basis(c(0,1),13)
#  compute the 13 by 13 matrix of inner products of second derivatives
penmat <- getbasispenalty(basisobj)
#  set up a Fourier basis with 13 basis functions
#  and and period 1.0.
basisobj <- create.fourier.basis(c(0,1),13)
#  compute the 13 by 13 matrix of inner products of second derivatives
penmat <- getbasispenalty(basisobj)

Extract the range from a basis object

Description

Extracts the 'range' component from basis object 'basisobj'.

Usage

getbasisrange(basisobj) 

Arguments

Value

a numeric vector of length 2

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

Berkeley Growth Study data

Description

A list containing the heights of 39 boys and 54 girls from age 1 to 18 and the ages at which they were collected.

Format

This list contains the following components:

hgtm

a 31 by 39 numeric matrix giving the heights in centimeters of 39 boys at 31 ages.

hgtf

a 31 by 54 numeric matrix giving the heights in centimeters of 54 girls at 31 ages.

age

a numeric vector of length 31 giving the ages at which the heights were measured.

Details

The ages are not equally spaced.

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

Tuddenham, R. D., and Snyder, M. M. (1954) "Physical growth of California boys and girls from birth to age 18", University of California Publications in Child Development, 1, 183-364.

Examples

with(growth, matplot(age, hgtf[, 1:10], type="b"))

Cursive handwriting samples

Description

20 cursive samples of 1401 (x, y,) coordinates for writing "fda"

Usage

handwrit
handwritTime

Format

handwrit

An array of dimensions (1401, 20, 2) giving 1401 pairs of (x, y) coordinates for each of 20 replicates of cursively writing "fda"

handwritTime

seq(0, 2300, length=1401) = sampling times

Details

These data are the X-Y coordinates of 20 replications of writing the script "fda". The subject was Jim Ramsay. Each replication is represented by 1401 coordinate values. The scripts have been extensively pre-processed. They have been adjusted to a common length that corresponds to 2.3 seconds or 2300 milliseconds, and they have already been registered so that important features in each script are aligned.

This analysis is designed to illustrate techniques for working with functional data having rather high frequency variation and represented by thousands of data points per record. Comments along the way explain the choices of analysis that were made.

The final result of the analysis is a third order linear differential equation for each coordinate forced by a constant and by time. The equations are able to reconstruct the scripts to a fairly high level of accuracy, and are also able to accommodate a substantial amount of the variation in the observed scripts across replications. by contrast, a second order equation was found to be completely inadequate.

An interesting surprise in the results is the role placed by a 120 millisecond cycle such that sharp features such as cusps correspond closely to this period. This 110-120 msec cycle seems is usually seen in human movement data involving rapid movements, such as speech, juggling and so on.

These 20 records have already been normalized to a common time interval of 2300 milliseconds and have beeen also registered so that prominent features occur at the same times across replications. Time will be measured in (approximate) milliseconds and space in meters. The data will require a small amount of smoothing, since an error of 0.5 mm is characteristic of the OPTOTRAK 3D measurement system used to collect the data.

Milliseconds were chosen as a time scale in order to make the ratio of the time unit to the inter-knot interval not too far from one. Otherwise, smoothing parameter values may be extremely small or extremely large.

The basis functions will be B-splines, with a spline placed at each knot. One may question whether so many basis functions are required, but this decision is found to be essential for stable derivative estimation up to the third order at and near the boundaries.

Order 7 was used to get a smooth third derivative, which requires penalizing the size of the 5th derivative, which in turn requires an order of at least 7. This implies norder + no. of interior knots = 1399 + 7 = 1406 basis functions.

The smoothing parameter value 1e8 was chosen to obtain a fitting error of about 0.5 mm, the known error level in the OPTOTRACK equipment.

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

Examples

oldpar <- par(no.readonly=TRUE)
plot(handwrit[, 1, 1], handwrit[, 1, 2], type="l")
par(oldpar)

Tibia Length for One Baby

Description

Measurement of the length of the tibia for the first 40 days of life for one infant.

Usage

data(infantGrowth)

Format

A matrix with three columns:

day

age in days

tibiaLength

The average of five measurements of tibia length in millimeters

sd.length

The standard deviation of five measurements of tibia length in millimeters

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

Hermanussen, M., Thiel, C., von Bueren, E., de los Angeles Rol. de Lama, M., Perez Romero, A., Arizonavarreta Ruiz, C., Burmeister, J., Tresguerras, J. A. F. (1998) Micro and macro perspectives in auxology: Findings and considerations upon the variability of short term and individual growth and the stability of population derived parameters, Annals of Human Biology, 25: 359-395.

Examples

data(infantGrowth)
oldpar <- par(no.readonly=TRUE)
plot(tibiaLength~day, infantGrowth)
par(oldpar)

Inner products of Functional Data Objects.

Description

Computes a matrix of inner products for each pairing of a replicate for the first argument with a replicate for the second argument. This is perhaps the most important function in the functional data library. Hardly any analysis fails to use inner products in some way, and many employ multiple inner products. While in certain cases these may be computed exactly, this is a more general function that approximates the inner product approximately when required. The inner product is defined by two derivatives or linear differential operators that are applied to the first two arguments. The range used to compute the inner product may be contained within the range over which the functions are defined. A weight functional data object may also be used to define weights for the inner product.

Usage

inprod(fdobj1, fdobj2,
       Lfdobj1=int2Lfd(0), Lfdobj2=int2Lfd(0), rng = range1, wtfd = 0)

Arguments

Details

The approximation method is Richardson extrapolation using numerical integration by the trapezoidal rule. At each iteration, the number of values at which the functions are evaluated is doubled, and a polynomial extrapolation method is used to estimate the converged integral values as well as an error tolerance. Convergence is declared when the relative error falls below EPS for all products. The extrapolation method generally saves at least one and often two iterations relative to un-extrapolated trapezoidal integration. Functional data analyses will seldom need to use inprod directly, but code developers should be aware of its pivotal role. Future work may require more sophisticated and specialized numerical integration methods. inprod computes the definite integral, but some functions such as smooth.monotone and register.fd also need to compute indefinite integrals. These use the same approximation scheme, but usually require more accuracy, and hence more iterations. When one or both arguments are basis objects, they are converted to functional data objects using identity matrices as the coefficient matrices. inprod is only called when there is no faster or exact method available. In cases where there is, it has been found that the approximation is good to about four to five significant digits, which is sufficient for most applications. Perhaps surprisingly, in the case of B-splines, the exact method is not appreciably faster, but of course is more accurate. inprod calls function eval.fd perhaps thousands of times, so high efficiency for this function and the functions that it calls is important.

Value

a matrix of inner products. The number of rows is the number of functions or basis functions in argument fd1, and the number of columns is the same thing for argument fd2.

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

Press, et, al, $Numerical Recipes$.

See Also

eval.penalty,

Compute Inner Products B-spline Expansions.

Description

Computes the matrix of inner products when both functions are represented by B-spline expansions and when both derivatives are integers. This function is called by function inprod, and is not normally used directly.

Usage

inprod.bspline(fdobj1, fdobj2=fdobj1, nderiv1=0, nderiv2=0)

Arguments

Value

a matrix of inner products with number of rows equal to the number of replications of the first argument and number of columns equal to the number of replications of the second object.

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

Convert Integer to Linear Differential Operator

Description

This function turns an integer specifying an order of a derivative into the equivalent linear differential operator object. It is also useful for checking that an object is of the "Lfd" class.

Usage

int2Lfd(m=0)

Arguments

Details

Smoothing is achieved by penalizing the integral of the square of the derivative of order m over rangeval:

m = 0 penalizes the squared difference from 0 of the function

1 = penalize the square of the slope or velocity

2 = penalize the squared acceleration

3 = penalize the squared rate of change of acceleration

4 = penalize the squared curvature of acceleration?

Value

a linear differential operator object of the "Lfd" class that is equivalent to the integer argument.

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

Examples

# Lfd to penalize the squared acceleration
# typical for smoothing a cubic spline (order 4)
int2Lfd(2)

# Lfd to penalize the curvature of acceleration
# used with splines of order 6
# when it is desired to study velocity and acceleration
int2Lfd(4)

Intensity Function for Point Process

Description

The intensity $mu$ of a series of event times that obey a homogeneous Poisson process is the mean number of events per unit time. When this event rate varies over time, the process is said to be nonhomogeneous, and $mu(t)$, and is estimated by this function intensity.fd.

Usage

intensity.fd(x, WfdParobj, conv=0.0001, iterlim=20,
             dbglev=1, returnMatrix=FALSE)

Arguments

Details

The intensity function $I(t)$ is almost the same thing as a probability density function $p(t)$ estimated by function densify.fd. The only difference is the absence of the normalizing constant $C$ that a density function requires in order to have a unit integral. The goal of the function is provide a smooth intensity function estimate that approaches some target intensity by an amount that is controlled by the linear differential operator Lfdobj and the penalty parameter in argument WfdPar. For example, if the first derivative of $W(t)$ is penalized heavily, this will force the function to approach a constant, which in turn will force the estimated Poisson process itself to be nearly homogeneous. To plot the intensity function or to evaluate it, evaluate Wfdobj, exponentiate the resulting vector.

Value

A named list of length 4 containing:

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

density.fd

Examples

oldpar <- par(no.readonly=TRUE)
#  Generate 101 Poisson-distributed event times with
#  intensity or rate two events per unit time
N  <- 101
mu <- 2
#  generate 101 uniform deviates
uvec <- runif(rep(0,N))
#  convert to 101 exponential waiting times
wvec <- -log(1-uvec)/mu
#  accumulate to get event times
tvec <- cumsum(wvec)
tmax <- max(tvec)
#  set up an order 4 B-spline basis over [0,tmax] with
#  21 equally spaced knots
tbasis <- create.bspline.basis(c(0,tmax), 23)
#  set up a functional parameter object for W(t),
#  the log intensity function.  The first derivative
#  is penalized in order to smooth toward a constant
lambda <- 10
Wfd0 <- fd(matrix(0,23,1),tbasis)
WfdParobj <- fdPar(Wfd0, 1, lambda)
#  estimate the intensity function
Wfdobj <- intensity.fd(tvec, WfdParobj)$Wfdobj
#  get intensity function values at 0 and event times
events <- c(0,tvec)
intenvec <- exp(eval.fd(events,Wfdobj))
#  plot intensity function
plot(events, intenvec, type="b")
lines(c(0,tmax),c(mu,mu),lty=4)
par(oldpar)

Confirm Object is Class "Basisfd"

Description

Check that an argument is a basis object.

Usage

is.basis(basisobj)

Arguments

Value

a logical value: TRUE if the class is correct, FALSE otherwise.

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

is.fd, is.fdPar, is.Lfd

Confirm that two objects of class "Basisfd" are identical

Description

Check all the slots of two basis objects to see that they are identical.

Usage

is.eqbasis(basisobj1, basisobj2)

Arguments

Value

a logical value: TRUE if the two basis objects are identical, FALSE otherwise.

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

basisfd,

Confirm Object has Class "fd"

Description

Check that an argument is a functional data object.

Usage

is.fd(fdobj)

Arguments

Value

a logical value: TRUE if the class is correct, FALSE otherwise.

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

is.basis, is.fdPar, is.Lfd

Confirm Object has Class "fdPar"

Description

Check that an argument is a functional parameter object.

Usage

is.fdPar(fdParobj)

Arguments

Value

a logical value: TRUE if the class is correct, FALSE otherwise.

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

is.basis, is.fd, is.Lfd

Confirm Object has Class "fdSmooth"

Description

Check that an argument is a functional parameter object.

Usage

is.fdSmooth(fdSmoothobj)

Arguments

Value

a logical value: TRUE if the class is correct, FALSE otherwise.

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

is.basis, is.fd, is.Lfd, is.fdPar

Confirm Object has Class "Lfd"

Description

Check that an argument is a linear differential operator object.

Usage

is.Lfd(Lfdobj)

Arguments

Value

a logical value: TRUE if the class is correct, FALSE otherwise.

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

is.basis, is.fd, is.fdPar

Extract the knots from a function basis or data object

Description

Extract either all or only the interior knots from an object of class basisfd, fd, or fdSmooth.

Usage

## S3 method for class 'fd'
knots(Fn, interior=TRUE, ...)
## S3 method for class 'fdSmooth'
knots(Fn, interior=TRUE, ...)
## S3 method for class 'basisfd'
knots(Fn, interior=TRUE, ...)

Arguments

Details

The interior knots of a bspline basis are stored in the params component. The remaining knots are in the rangeval component, with multiplicity norder(Fn).

Value

Numeric vector. If 'interior' is TRUE, this is the params component of the bspline basis. Otherwise, params is bracketed by rep(rangeval, norder(basisfd).

Author(s)

Spencer Graves

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

fd, create.bspline.basis,

Convert Smoothing Parameter to Degrees of Freedom

Description

The degree of roughness of an estimated function is controlled by a smoothing parameter $lambda$ that directly multiplies the penalty. However, it can be difficult to interpret or choose this value, and it is often easier to determine the roughness by choosing a value that is equivalent of the degrees of freedom used by the smoothing procedure. This function converts a multipler $lambda$ into a degrees of freedom value.

Usage

lambda2df(argvals, basisobj, wtvec=rep(1, n), Lfdobj=NULL, lambda=0)

Arguments

Value

the equivalent degrees of freedom value.

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

df2lambda

Compute GCV Criterion

Description

The generalized cross-validation or GCV criterion is often used to select an appropriate smoothing parameter value, by finding the smoothing parameter that minimizes GCV. This function locates that value.

Usage

  lambda2gcv(log10lambda, argvals, y, fdParobj, wtvec=rep(1,length(argvals)))

Arguments

Details

Currently, lambda2gcv

Value

1. fdParobj[['lambda']] <- 10^log10lambda

2. smoothlist <- smooth.basks(argvals, y, fdParobj, wtvec)

3. return(smoothlist[['gcv']])

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

smooth.basis fdPar

Landmark Registration of Functional Observations with Differing Ranges

Description

It is common to see that among a set of functions certain prominent features such peaks and valleys, called $landmarks$, do not occur at the same times, or other argument values. This is called $phase variation$, and it can be essential to align these features before proceeding with further functional data analyses.

This function uses the timings of these features to align or register the curves. The registration involves estimating a nonlinear transformation of the argument continuum for each functional observation. This transformation is called a warping function. It must be strictly increasing and smooth.

Warning: As of March 2022, landmark registration cannot be done using function smooth.basis instead of function smooth.morph. The warping function must be strictly monotonic, and we have found that using smooth.basis too often violates this contraint. Function smooth.morph ensures monotonicity.

Usage

landmarkreg(unregfd, ximarks, x0marks, x0lim, 
             WfdPar=NULL, WfdPar0=NULL, ylambda=1e-10)

Arguments

Details

A value of an arbitrary strictly monotone function at a point x can be defined as the indefinite integral from a fixed lower boundary to x of the exponential of an unconstrained function value W(x).

We use B-spline basis functions to define function W, and optimize the coefficients of its B-spline expansion with respect to the data we are fitting. The optimized core function W is output along with the registered functions, the warping function qnd the inverse warping function.

It is essential that the location of every landmark be clearly defined in each of the curves as well as the template function. If this is not the case, consider using the continuous registration function register.fd.

Not much curvature is usually required in the warping functions, so a rather lower power basis, usually B-splines, is suitable for defining the functional parameter argument WfdPar. A registration with a few prominent landmarks is often a good preliminary to using the more sophisticated but more lengthy process in register.fd.

Value

a named list of length 4 with components:

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

register.fd, landmarkreg, smooth.morph

Define a Linear Differential Operator Object

Description

A linear differential operator of order $m$ is defined, usually to specify a roughness penalty.

Usage

Lfd(nderiv=0, bwtlist=vector("list", 0))

Arguments

Details

To check that an object is of this class, use functions is.Lfd or int2Lfd.

Linear differential operator objects are often used to define roughness penalties for smoothing towards a "hypersmooth" function that is annihilated by the operator. For example, the harmonic acceleration operator used in the analysis of the Canadian daily weather data annihilates linear combinations of $1, sin(2 pi t/365)$ and $cos(2 pi t/365)$, and the larger the smoothing parameter, the closer the smooth function will be to a function of this shape.

Function pda.fd estimates a linear differential operator object that comes as close as possible to annihilating a functional data object.

A linear differential operator of order $m$ is a linear combination of the derivatives of a functional data object up to order $m$. The derivatives of orders 0, 1, ..., $m-1$ can each be multiplied by a weight function $b(t)$ that may or may not vary with argument $t$.

If the notation $D^j$ is taken to mean "take the derivative of order $j$", then a linear differental operator $L$ applied to function $x$ has the expression

$Lx(t) = b_0(t) x(t) + b_1(t)Dx(t) + ... + b_{m-1}(t) D^{m-1} x(t) + D^mx(t)$

There are print, summary, and plot methods for objects of class Lfd.

Value

a linear differential operator object

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

int2Lfd, vec2Lfd, fdPar, pda.fd plot.Lfd

Examples

#  Set up the harmonic acceleration operator
dayrange  <- c(0,365)
Lbasis  <- create.constant.basis(dayrange,
                  axes=list("axesIntervals"))
Lcoef   <- matrix(c(0,(2*pi/365)^2,0),1,3)
bfdobj  <- fd(Lcoef,Lbasis)
bwtlist <- fd2list(bfdobj)
harmaccelLfd <- Lfd(3, bwtlist)

Add Lines from Functional Data to a Plot

Description

Lines defined by functional observations are added to an existing plot.

Usage

## S3 method for class 'fd'
lines(x, Lfdobj=int2Lfd(0), nx=201, ...)
## S3 method for class 'fdSmooth'
lines(x, Lfdobj=int2Lfd(0), nx=201, ...)

Arguments

Side Effects

Lines added to an existing plot.

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

plot.fd, plotfit.fd

Examples

oldpar <- par(no.readonly=TRUE)
##
## plot a fit with 3 levels of smoothing
##
x <- seq(-1,1,0.02)
y <- x + 3*exp(-6*x^2) + sin(1:101)/2
# sin not rnorm to make it easier to compare
# results across platforms 

result4.0  <- smooth.basisPar(argvals=x, y=y, lambda=1)
result4.m4 <- smooth.basisPar(argvals=x, y=y, lambda=1e-4)

plot(x, y)
lines(result4.0$fd)
lines(result4.m4$fd,      col='blue')
lines.fdSmooth(result4.0, col='red') 

plot(x, y, xlim=c(0.5, 1))
lines.fdSmooth(result4.0)
lines.fdSmooth(result4.m4, col='blue')
# no visible difference from the default?  
par(oldpar)

Fit Fully Functional Linear Model

Description

A functional dependent variable y_i(t) is approximated by a single functional covariate x_i(s) plus an intercept function \alpha(t), and the covariate can affect the dependent variable for all values of its argument. The equation for the model is

y_i(t) = \beta_0(t) + \int \beta_1(s,t) x_i(s) ds + e_i(t)

for i = 1,...,N. The regression function \beta_1(s,t) is a bivariate function. The final term e_i(t) is a residual, lack of fit or error term. There is no need for values s and t to be on the same continuum.

Usage

linmod(xfdobj, yfdobj, betaList, wtvec=NULL)

Arguments

Value

a named list of length 3 with the following entries:

Source

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009) Functional Data Analysis in R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

bifdPar, fRegress

Examples

#See the prediction of precipitation using temperature as
#the independent variable in the analysis of the daily weather
#data, and the analysis of the Swedish mortality data.

Lip motion

Description

51 measurements of the position of the lower lip every 7 milliseconds for 20 repitions of the syllable 'bob'.

Usage

lip
lipmarks
liptime

Format

lip

a matrix of dimension c(51, 20) giving the position of the lower lip every 7 milliseconds for 350 milliseconds.

lipmarks

a matrix of dimension c(20, 2) giving the positions of the 'leftElbow' and 'rightElbow' in each of the 20 repetitions of the syllable 'bob'.

liptime

time in seconds from the start = seq(0, 0.35, 51) = every 7 milliseconds.

Details

These are rather simple data, involving the movement of the lower lip while saying "bob". There are 20 replications and 51 sampling points. The data are used to illustrate two techniques: landmark registration and principal differental analysis. Principal differential analysis estimates a linear differential equation that can be used to describe not only the observed curves, but also a certain number of their derivatives. For a rather more elaborate example of principal differential analysis, see the handwriting data.

See the lip demo.

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

Examples

#  See the lip demo.

Search along a line for a minimum within an optimisation algorithm.

Description

This is a version of the function in Numerical Recipes. It is initialized with a function value and gradient, and it does a series of quadratic searches until a convergence criterion is reached. This version includes code for display the progress of iteration for debugging purposes.

Usage

lnsrch(xold, fold, g, p, func, dataList, stpmax, itermax=20, TOLX=1e-10, dbglev=0)

Arguments

Value

A named list containing:

Author(s)

Juan Li and James Ramsay

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

Press, W. H., Taukolsky, S. A., Vetterline, W. T. and Flannery, B. P. (2020) Numerical Recipes, Third Edition, Cambridge.

See Also

smooth.surp

Plot Columns of Matrices

Description

Plot the columns of one matrix against the columns of another.

Usage

matplot(x, ...)
## Default S3 method:
matplot(x, y, type = "p", lty = 1:5, lwd = 1,
    lend = par("lend"), pch = NULL, col = 1:6, cex = NULL, bg = NA,
    xlab = NULL, ylab = NULL, xlim = NULL, ylim = NULL, ..., add = FALSE,
    verbose = getOption("verbose"))
## S3 method for class 'Date'
matplot(x, y, type = "p", lty = 1:5, lwd = 1,
    lend = par("lend"), pch = NULL, col = 1:6, cex = NULL, bg = NA,
    xlab = NULL, ylab = NULL, xlim = NULL, ylim = NULL, ..., add = FALSE,
    verbose = getOption("verbose"))
## S3 method for class 'POSIXct'
matplot(x, y, type = "p", lty = 1:5, lwd = 1,
    lend = par("lend"), pch = NULL, col = 1:6, cex = NULL, bg = NA,
    xlab = NULL, ylab = NULL, xlim = NULL, ylim = NULL, ..., add = FALSE,
    verbose = getOption("verbose"))

Arguments

Details

Note that for multivariate data, a suitable array must first be defined using the par function.

matplot.default calls matplot. The other methods are needed, because the default methods ignore the Date or POSIXct character of x, labeling the horizontal axis as numbers, thereby placing it on the user to translate the numbers of days or seconds since the start of the epoch into dates (and possibly times for POSIXct x).

Side Effects

a plot of the functional observations

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

matplot, plot, points, lines, matrix, par

Examples

oldpar <- par(no.readonly=TRUE)
##
## matplot.Date, matplot.POSIXct
##
# Date
invasion1 <- as.Date('1775-09-04')
invasion2 <- as.Date('1812-07-12')
earlyUS.Canada <- c(invasion1, invasion2)
Y <- matrix(1:4, 2, 2)
matplot(earlyUS.Canada, Y)

# POSIXct
AmRev.ct <- as.POSIXct1970(c('1776-07-04', '1789-04-30'))
matplot(AmRev.ct, Y)

##
## matplot.default (copied from matplot{graphics})
##
matplot((-4:5)^2, main = "Quadratic") # almost identical to plot(*)
sines <- outer(1:20, 1:4, function(x, y) sin(x / 20 * pi * y))
matplot(sines, pch = 1:4, type = "o", col = rainbow(ncol(sines)))
matplot(sines, type = "b", pch = 21:23, col = 2:5, bg = 2:5,
             main = "matplot(...., pch = 21:23, bg = 2:5)")

x <- 0:50/50
matplot(x, outer(x, 1:8, function(x, k) sin(k*pi * x)),
             ylim = c(-2,2), type = "plobcsSh",
             main= "matplot(,type = \"plobcsSh\" )")
## pch & type =  vector of 1-chars :
matplot(x, outer(x, 1:4, function(x, k) sin(k*pi * x)),
             pch = letters[1:4], type = c("b","p","o"))

lends <- c("round","butt","square")
matplot(matrix(1:12, 4), type="c", lty=1, lwd=10, lend=lends)
text(cbind(2.5, 2*c(1,3,5)-.4), lends, col= 1:3, cex = 1.5)

table(iris$Species) # is data.frame with 'Species' factor
iS <- iris$Species == "setosa"
iV <- iris$Species == "versicolor"
op <- par(bg = "bisque")
matplot(c(1, 8), c(0, 4.5), type= "n", xlab = "Length", ylab = "Width",
             main = "Petal and Sepal Dimensions in Iris Blossoms")
matpoints(iris[iS,c(1,3)], iris[iS,c(2,4)], pch = "sS", col = c(2,4))
matpoints(iris[iV,c(1,3)], iris[iV,c(2,4)], pch = "vV", col = c(2,4))
legend(1, 4, c("    Setosa Petals", "    Setosa Sepals",
                    "Versicolor Petals", "Versicolor Sepals"),
            pch = "sSvV", col = rep(c(2,4), 2))

nam.var <- colnames(iris)[-5]
nam.spec <- as.character(iris[1+50*0:2, "Species"])
iris.S <- array(NA, dim = c(50,4,3),
                     dimnames = list(NULL, nam.var, nam.spec))
for(i in 1:3) iris.S[,,i] <- data.matrix(iris[1:50+50*(i-1), -5])

matplot(iris.S[,"Petal.Length",], iris.S[,"Petal.Width",], pch="SCV",
             col = rainbow(3, start = .8, end = .1),
             sub = paste(c("S", "C", "V"), dimnames(iris.S)[[3]],
                         sep = "=", collapse= ",  "),
             main = "Fisher's Iris Data")
par(op)
par(oldpar)

Mean of Functional Data

Description

Evaluate the mean of a set of functions in a functional data object.

Usage

## S3 method for class 'fd'
mean(x, ...)

Arguments

Value

a functional data object with a single replication that contains the mean of the functions in the object fd.

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

stddev.fd, var.fd, sum.fd, center.fd mean

Examples

oldpar <- par(no.readonly=TRUE)
##
## 1.  univeriate:  lip motion
##
liptime  <- seq(0,1,.02)
liprange <- c(0,1)

#  -------------  create the fd object -----------------
#       use 31 order 6 splines so we can look at acceleration

nbasis <- 51
norder <- 6
lipbasis <- create.bspline.basis(liprange, nbasis, norder)

#  ------------  apply some light smoothing to this object  -------

lipLfdobj <- int2Lfd(4)
lipLambda <- 1e-12
lipfdPar  <- fdPar(fd(matrix(0,nbasis,1),lipbasis), lipLfdobj, lipLambda)

lipfd <- smooth.basis(liptime, lip, lipfdPar)$fd
names(lipfd$fdnames) = c("Normalized time", "Replications", "mm")

lipmeanfd <- mean.fd(lipfd)
plot(lipmeanfd)

##
## 2.  Trivariate:  CanadianWeather
##
dayrng <- c(0, 365)

nbasis <- 51
norder <- 6

weatherBasis <- create.fourier.basis(dayrng, nbasis)

weather.fd <- smooth.basis(day.5, CanadianWeather$dailyAv,
                           weatherBasis)$fd

str(weather.fd.mean <- mean.fd(weather.fd))
par(oldpar)

melanoma 1936-1972

Description

These data from the Connecticut Tumor Registry present age-adjusted numbers of melanoma skin-cancer incidences per 100,000 people in Connectict for the years from 1936 to 1972.

Format

A data frame with 37 observations on the following 2 variables.

year

Years 1936 to 1972.

incidence

Rate of melanoma cancer per 100,000 population.

Details

This is a copy of the 'melanoma' dataset in the 'lattice' package. It is unrelated to the object of the same name in the 'boot' package.

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

melanoma melanoma

Examples

oldpar <- par(no.readonly=TRUE)
plot(melanoma[, -1], type="b")
par(oldpar)

Evaluate the a monotone function

Description

Evaluate a monotone function defined as the indefinite integral of $exp(W(t))$ where $W$ is a function defined by a basis expansion. Function $W$ is the logarithm of the derivative of the monotone function.

Usage

  monfn(argvals, Wfdobj, basislist=vector("list", JMAX), returnMatrix=FALSE)

Arguments

Details

This function evaluates a strictly monotone function of the form

h(x) = [D^{-1} exp(Wfdobj)](x),

where D^{-1} means taking the indefinite integral. The interval over which the integration takes places is defined in the basis object in Wfdobj.

Value

A numerical vector or matrix containing the values the warping function h.

Author(s)

J. O. Ramsay

References

Ramsay, James O., Hooker, G. and Graves, S. (2009), Functional Data Analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

mongrad, landmarkreg, smooth.morph

Examples

oldpar <- par(no.readonly=TRUE)
## basically this example resembles part of landmarkreg.R that uses monfn.R to
## estimate the warping function.

## Specify the curve subject to be registered
n=21
tbreaks = seq(0, 2*pi, len=n)
xval <- sin(tbreaks)
rangeval <- range(tbreaks)

## Establish a B-spline basis for the curve
wbasis <- create.bspline.basis(rangeval=rangeval, breaks=tbreaks)
Wfd0   <- fd(matrix(0,wbasis$nbasis,1),wbasis)
WfdPar <- fdPar(Wfd0, 1, 1e-4)
fdObj  <- smooth.basis(tbreaks, xval, WfdPar)$fd

## Set the mean landmark times. Note that the objective of the warping
## function is to transform the curve such that the landmarks of the curve
## occur at the designated mean landmark times.

## Specify the mean landmark times: tbreak[8]=2.2 and tbreaks[13]=3.76
meanmarks <- c(rangeval[1], tbreaks[8], tbreaks[13], rangeval[2])
## Specify landmark locations of the curve: tbreaks[6] and tbreaks[16]
cmarks <- c(rangeval[1], tbreaks[6], tbreaks[16], rangeval[2])

## Establish a B-basis object for the warping function
Wfd <- smooth.morph(x=meanmarks, y=cmarks, ylim=rangeval, 
                    WfdPar=WfdPar)$Wfdobj

## Estimate the warping function
h = monfn(tbreaks, Wfd)

## scale using a linear equation h such that h(0)=0 and h(END)=END
b <- (rangeval[2]-rangeval[1])/ (h[n]-h[1])
a <- rangeval[1] - b*h[1]
h <- a + b*h
plot(tbreaks, h, xlab="Time", ylab="Transformed time", type="b")
par(oldpar)

Evaluate the gradient of a monotone function

Description

Evaluates the gradient of a monotone function with respect to the coefficients defining the log-first derivative $W(t)$ at each of a set of argument values.

Usage

mongrad(x, Wfdobj, basislist=vector("list",JMAX), 
                    returnMatrix=FALSE)

Arguments

Value

A matrix with as many rows as argument values and as many columns as basis functions defining $W$.

Author(s)

J. O. Ramsay

References

Ramsay, James O., Hooker, G. and Graves, S. (2009), Functional Data Analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

monfn, monhess, landmarkreg, smooth.morph

Evaluate the Hessian matrix of a monotone function

Description

Evaluates the hessian or second derivative matrix of a monotone function with respect to the coefficients definingthe log-first derivative $W(t)$ at each of a set of argument values.

Usage

monhess(x, Wfd, basislist)

Arguments

Value

A three dimensional array with first dimension corresponding to argument values and second and third dimensions to number of basis functions defining $W$.

Author(s)

J. O. Ramsay

References

Ramsay, James O., Hooker, G. and Graves, S. (2009), Functional Data Analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

monfn, mongrad, landmarkreg, smooth.morph

Evaluate Monomial Basis

Description

Computes the values of the powers of argument t.

Usage

monomial(evalarg, exponents=1, nderiv=0, argtrans=c(0,1))

Arguments

Value

a matrix of values of basis functions. Rows correspond to argument values and columns to basis functions.

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

power, expon, fourier, polyg, bsplineS

Examples

# set up a monomial basis for the first five powers
nbasis   <- 5
basisobj <- create.monomial.basis(c(-1,1),nbasis)
#  evaluate the basis
tval <- seq(-1,1,0.1)
basismat <- monomial(tval, 1:basisobj$nbasis)

Evaluate Monomial Roughness Penalty Matrix

Description

The roughness penalty matrix is the set of inner products of all pairs of a derivative of integer powers of the argument.

Usage

monomialpen(basisobj, Lfdobj=int2Lfd(2),
            rng=basisobj$rangeval)

Arguments

Value

a symmetric matrix of order equal to the number of monomial basis functions.

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

exponpen, fourierpen, bsplinepen, polygpen

Examples

##
## set up a monomial basis for the first five powers
##
nbasis   <- 5
basisobj <- create.monomial.basis(c(-1,1),nbasis)
#  evaluate the rougness penalty matrix for the
#  second derivative.
penmat <- monomialpen(basisobj, 2)

##
## with rng of class Date and POSIXct
##
# Date
invasion1 <- as.Date('1775-09-04')
invasion2 <- as.Date('1812-07-12')
earlyUS.Canada <- c(invasion1, invasion2)
BspInvade1 <- create.monomial.basis(earlyUS.Canada)
invadmat <- monomialpen(BspInvade1)

# POSIXct
AmRev.ct <- as.POSIXct1970(c('1776-07-04', '1789-04-30'))
BspRev1.ct <- create.monomial.basis(AmRev.ct)
revmat <- monomialpen(BspRev1.ct)

Montreal Daily Temperature

Description

Temperature in degrees Celsius in Montreal each day from 1961 through 1994

Usage

data(MontrealTemp)

Format

A numeric array with dimnames = list(1961:1994, names(dayOfYear)).

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

CanadianWeather monthAccessories

Examples

data(MontrealTemp)

JanuaryThaw <- t(MontrealTemp[, 16:47])

Nondurable goods index

Description

US nondurable goods index time series, January 1919 to January 2000.

Format

An object of class 'ts'.

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

Examples

oldpar <- par(no.readonly=TRUE)
plot(nondurables, log="y")
par(oldpar)

Order of a B-spline

Description

norder = number of basis functions minus the number of interior knots.

Usage

norder(x, ...)
## S3 method for class 'fd'
norder(x, ...)
## S3 method for class 'basisfd'
norder(x, ...)
## Default S3 method:
norder(x, ...)

#norder.bspline(x, ...)

#NOTE:  The following is required by CRAN rules that
# function names like "as.numeric" must follow the documentation
# standards for S3 generics, even when they are not.
# Please ignore the following line:
## S3 method for class 'bspline'
norder(x, ...)

Arguments

Details

norder throws an error of basisfd[['type']] != 'bspline'.

Value

An integer giving the order of the B-spline.

Author(s)

Spencer Graves

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

create.bspline.basis

Examples

  bspl1.1 <- create.bspline.basis(norder=1, breaks=0:1)

  stopifnot(norder(bspl1.1)==1)

  stopifnot(norder(fd(0, basisobj=bspl1.1))==1)

  stopifnot(norder(fd(rep(0,4)))==4)

  stopifnot(norder(list(fd(rep(0,4))))==4)
  ## Not run: 
    norder(list(list(fd(rep(0,4)))))
    Error in norder.default(list(list(fd(rep(0, 4))))) :
    input is not a 'basisfd' object and does not have a 'basisfd'
    component.
  
## End(Not run)

  stopifnot(norder(create.bspline.basis(norder=1, breaks=c(0,.5, 1))) == 1)

  stopifnot(norder(create.bspline.basis(norder=2, breaks=c(0,.5, 1))) == 2)

  # Default B-spline basis:  Cubic spline:  degree 3, order 4,
  # 21 breaks, 19 interior knots.
  stopifnot(norder(create.bspline.basis()) == 4)

  # these five lines throw an error of for a nocran check 
  # if (!CRAN()) {
  #   norder(create.fourier.basis(c(0,12) ))
    # Error in norder.bspline(x) :
    # object x is of type = fourier;  'norder' is only defined for type = 'bspline'
  # }

Add names to an object

Description

Add names to an object from 'preferred' if available and 'default' if not.

Usage

objAndNames(object, preferred, default)

Arguments

Details

1. If length(dim(object))<1, names(object) are taken from 'preferred' if they are not NULL and have the correct length, else try 'default'.

2. Else for(lvl in 1:length(dim(object))) take dimnames[[lvl]] from 'preferred[[i]]' if they are not NULL and have the correct length, else try 'default[[lvl]].

Value

An object of the same class and structure as 'object' but with either names or dimnames added or changed.

Author(s)

Spencer Graves

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

bifd

Examples

# The following should NOT check 'anything' here
tst1 <- objAndNames(1:2, list(letters[1:2], LETTERS[1:2]), anything)
all.equal(tst1, c(a=1, b=2))

# The following should return 'object unchanged
tst2 <- objAndNames(1:2, NULL, list(letters))
all.equal(tst2, 1:2)


tst3 <- objAndNames(1:2, list("a", 2), list(letters[1:2]))
all.equal(tst3, c(a=1, b=2) )

# The following checks a matrix / array
tst4 <- array(1:6, dim=c(2,3))
tst4a <- tst4
dimnames(tst4a) <- list(letters[1:2], LETTERS[1:3])
tst4b <- objAndNames(tst4, 
       list(letters[1:2], LETTERS[1:3]), anything)
all.equal(tst4b, tst4a)

tst4c <- objAndNames(tst4, NULL,        
       list(letters[1:2], LETTERS[1:3]) )
all.equal(tst4c, tst4a)

Numerical Solution mth Order Differential Equation System

Description

The system of differential equations is linear, with possibly time-varying coefficient functions. The numerical solution is computed with the Runge-Kutta method.

Usage

odesolv(bwtlist, ystart=diag(rep(1,norder)),
        h0=width/100, hmin=width*1e-10, hmax=width*0.5,
        EPS=1e-4, MAXSTP=1000)

Arguments

Details

This function is required to compute a set of solutions of an estimated linear differential equation in order compute a fit to the data that solves the equation. Such a fit will be a linear combinations of m independent solutions.

Value

a named list of length 2 containing

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

pda.fd. For new applications, users are encouraged to consider deSolve. The deSolve package provides general solvers for ordinary and partial differential equations, as well as differential algebraic equations and delay differential equations.

Examples

#See the analyses of the lip data.

Functional Principal Components Analysis

Description

Functional Principal components analysis aims to display types of variation across a sample of functions. Principal components analysis is an exploratory data analysis that tends to be an early part of many projects. These modes of variation are called $principal components$ or $harmonics.$ This function computes these harmonics, the eigenvalues that indicate how important each mode of variation, and harmonic scores for individual functions. If the functions are multivariate, these harmonics are combined into a composite function that summarizes joint variation among the several functions that make up a multivariate functional observation.

Usage

pca.fd(fdobj, nharm = 2, harmfdPar=fdPar(fdobj),
       centerfns = TRUE)

Arguments

Value

an object of class "pca.fd" with these named entries:

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

cca.fd, pda.fd

Examples

oldpar <- par(no.readonly=TRUE)
#  carry out a PCA of temperature
#  penalize harmonic acceleration, use varimax rotation

daybasis65 <- create.fourier.basis(c(0, 365), nbasis=65, period=365)
nbasis <-65
harmaccelLfd <- vec2Lfd(c(0,(2*pi/365)^2,0), c(0, 365))
harmfdPar    <- fdPar(fd(matrix(0,nbasis,1), daybasis65), harmaccelLfd, 
                      lambda=1e5)
daytempfd <- smooth.basis(day.5, CanadianWeather$dailyAv[,,"Temperature.C"],
                   daybasis65, fdnames=list("Day", "Station", "Deg C"))$fd

daytemppcaobj <- pca.fd(daytempfd, nharm=4, harmfdPar)
daytemppcaVarmx <- varmx.pca.fd(daytemppcaobj)
#  plot harmonics
op <- par(mfrow=c(2,2))
plot.pca.fd(daytemppcaobj, cex.main=0.9)

plot.pca.fd(daytemppcaVarmx, cex.main=0.9)
par(op)

plot(daytemppcaobj$harmonics)
plot(daytemppcaVarmx$harmonics)
par(oldpar)

Estimate the functional principal components

Description

Carries out a functional PCA with regularization from the estimate of the covariance surface

Usage

  pcaPACE(covestimate, nharm, harmfdPar, cross)

Arguments

Value

an object of class "pca.fd" with these named entries:

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

Yao, F., Mueller, H.G., Wang, J.L. (2005), Functional data analysis for sparse longitudinal data, J. American Statistical Association, 100, 577-590.