Package 'fda'

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.

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

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)

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)) 
  # }

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)

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))

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)

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))

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)

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.

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)

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)

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)

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)

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)

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)

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)

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)

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]))

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.

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")
  #   }
  # }

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

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)

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))

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)

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)

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)

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)

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)

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.

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)

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,

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)