library(fda)
## Loading required package: splines
## Loading required package: Matrix
##
## Attaching package: 'fda'
## The following object is masked from 'package:graphics':
##
## matplot
##
## simplest call
##
girlGrowthSm <- with(growth, smooth.basisPar(argvals=age, y=hgtf, lambda=0.1))
plot(girlGrowthSm$fd, xlab="age", ylab="height (cm)",
main="Girls in Berkeley Growth Study" )
## [1] "done"
plot(deriv(girlGrowthSm$fd), xlab="age", ylab="growth rate (cm / year)",
main="Girls in Berkeley Growth Study" )
## [1] "done"
plot(deriv(girlGrowthSm$fd, 2), xlab="age",
ylab="growth acceleration (cm / year^2)",
main="Girls in Berkeley Growth Study" )
## [1] "done"
# Undersmoothed with lambda = 0
##
## Another simple call
##
lipSm <- smooth.basisPar(liptime, lip, lambda=1e-9)$fd
plot(lipSm)
## [1] "done"
##
## A third example
##
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
# set up a saturated B-spline basis
basisobj101 <- create.bspline.basis(x)
fdParobj101 <- fdPar(basisobj101, 2, lambda=1)
result101 <- smooth.basis(x, y, fdParobj101)
resultP <- smooth.basisPar(argvals=x, y=y, fdobj=basisobj101, lambda=1)
all.equal(result101, resultP)
## [1] TRUE
# TRUE
result4 <- smooth.basisPar(argvals=x, y=y, fdobj=4, lambda=1)
all.equal(resultP, result4)
## [1] TRUE
# TRUE
result4. <- smooth.basisPar(argvals=x, y=y, lambda=1)
all.equal(resultP, result4.)
## [1] TRUE
# TRUE
with(result4, c(df, gcv)) # display df and gcv measures
## rep1
## 2.900354 0.484053
result4.4 <- smooth.basisPar(argvals=x, y=y, lambda=1e-4)
with(result4.4, c(df, gcv)) # display df and gcv measures
## rep1
## 19.9640501 0.1621682
# less smoothing, more degrees of freedom, better fit
plot(result4.4)
## [1] "done"
lines(result4, col='green')
lines(result4$fd, col='green') # same as lines(result4, ...)
##
## fdnames?
##
girlGrow12 <- with(growth, smooth.basisPar(argvals=age, y=hgtf[, 1:2],
fdnames=c('age', 'girl', 'height'), lambda=0.1) )
girlGrow12. <- with(growth, smooth.basisPar(argvals=age, y=hgtf[, 1:2],
fdnames=list(age=age, girl=c('Carol', 'Sally'), value='height'),
lambda = 0.1) )
##
## Fourier basis with harmonic acceleration operator
##
daybasis65 <- create.fourier.basis(rangeval=c(0, 365), nbasis=65)
daytemp.fdSmooth <- with(CanadianWeather, smooth.basisPar(day.5,
dailyAv[,,"Temperature.C"],
daybasis65, fdnames=list("Day", "Station", "Deg C")) )
library(fda)
##
## 1. B-spline example
##
# Shows the effects of two levels of smoothing
# where the size of the third derivative is penalized.
# The null space contains quadratic functions.
x <- seq(-1,1,0.02)
y <- x + 3*exp(-6*x^2) + rnorm(rep(1,101))*0.2
# set up a saturated B-spline basis
basisobj <- create.bspline.basis(c(-1,1),81)
# convert to a functional data object that interpolates the data.
result <- smooth.basis(x, y, basisobj)
yfd <- result$fd
# set up a functional parameter object with smoothing
# parameter 1e-6 and a penalty on the 2nd derivative.
yfdPar <- fdPar(yfd, 2, 1e-6)
yfd1 <- smooth.fd(yfd, yfdPar)
yfd1. <- smooth.fdPar(yfd, 2, 1e-6)
all.equal(yfd1, yfd1.)
## [1] TRUE
# TRUE
# set up a functional parameter object with smoothing
# parameter 1 and a penalty on the 2nd derivative.
yfd2 <- smooth.fdPar(yfd, 2, 1)
# plot the data and smooth
plot(x,y) # plot the data
lines(yfd1, lty=1) # add moderately penalized smooth
lines(yfd2, lty=3) # add heavily penalized smooth
legend(-1,3,c("0.000001","1"),lty=c(1,3))
# plot the data and smoothing using function plotfit.fd
plotfit.fd(y, x, yfd1) # plot data and smooth
##
## 2. Fourier basis with harmonic acceleration operator
##
daybasis65 <- create.fourier.basis(rangeval=c(0, 365), nbasis=65)
daySmooth <- with(CanadianWeather, smooth.basis(day.5,
dailyAv[,,"Temperature.C"],
daybasis65, fdnames=list("Day", "Station", "Deg C")) )
daySmooth2 <- smooth.fdPar(daySmooth$fd, lambda=1e5)
op <- par(mfrow=c(2,1))
plot(daySmooth)
## [1] "done"
plot(daySmooth2)
## [1] "done"
par(op)
library(fda)
data(StatSciChinese)
i <- 3
StatSci1 <- StatSciChinese[, i, ]
# Where does the pen leave the paper?
plot(StatSci1[, 3], type='l')
thresh <- quantile(StatSci1[, 3], .8)
abline(h=thresh)
sel1 <- (StatSci1[, 3] < thresh)
StatSci1[!sel1, 1:2] <- NA
plot(StatSci1[, 1:2], type='l')
mark <- seq(1, 601, 12)
points(StatSci1[mark, 1], StatSci1[mark, 2])
library(fda)
# This tests the difference between boys and girls heights in the
# Berkeley growth data.
# First set up a basis system to hold the smooths
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
growfdPar <- fdPar(hgtbasis, Lfdobj, lambda)
hgtmfd <- smooth.basis(growth$age, growth$hgtm, growfdPar)$fd
hgtffd <- smooth.basis(growth$age, growth$hgtf, growfdPar)$fd
# Call tperm.fd
tres <- tperm.fd(hgtmfd,hgtffd)
library(fda)
##
## Example with 2 different bases
##
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 )
# Compare with structure described above under 'value':
str(tempPrecVar3.5 <- var.fd(tempfd3, precfd5))
## List of 4
## $ coefs : num [1:3, 1:5] 468.9 27.8 229.8 113.3 12.6 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:3] "const" "sin" "cos"
## .. ..$ : chr [1:5] "const" "sin1" "cos1" "sin2" ...
## $ sbasis :List of 10
## ..$ call : language basisfd(type = type, rangeval = rangeval, nbasis = nbasis, params = params, dropind = dropind, quadvals = qu| __truncated__
## ..$ type : chr "fourier"
## ..$ rangeval : num [1:2] 0 365
## ..$ nbasis : num 3
## ..$ params : num 365
## ..$ dropind : num(0)
## ..$ quadvals : NULL
## ..$ values : list()
## ..$ basisvalues: list()
## ..$ names : chr [1:3] "const" "sin" "cos"
## ..- attr(*, "class")= chr "basisfd"
## $ tbasis :List of 10
## ..$ call : language basisfd(type = type, rangeval = rangeval, nbasis = nbasis, params = params, dropind = dropind, quadvals = qu| __truncated__
## ..$ type : chr "fourier"
## ..$ rangeval : num [1:2] 0 365
## ..$ nbasis : num 5
## ..$ params : num 365
## ..$ dropind : num(0)
## ..$ quadvals : NULL
## ..$ values : list()
## ..$ basisvalues: list()
## ..$ names : chr [1:5] "const" "sin1" "cos1" "sin2" ...
## ..- attr(*, "class")= chr "basisfd"
## $ bifdnames:List of 2
## ..$ : chr [1:3] "const" "sin" "cos"
## ..$ : chr [1:5] "const" "sin1" "cos1" "sin2" ...
## - attr(*, "class")= chr "bifd"
##
## Example with 2 variables, same bases
##
gaitbasis3 <- create.fourier.basis(nbasis=3)
str(gaitfd3 <- Data2fd(gait, basisobj=gaitbasis3))
## 'y' is missing, using 'argvals'
## 'argvals' is missing; using seq( 0 , 1 , length= 20 )
## List of 3
## $ coefs : num [1:3, 1:39, 1:2] 23.86 -3.92 11.87 25.99 -3.57 ...
## ..- attr(*, "dimnames")=List of 3
## .. ..$ : NULL
## .. ..$ : chr [1:39] "boy1" "boy2" "boy3" "boy4" ...
## .. ..$ : chr [1:2] "Hip Angle" "Knee Angle"
## $ basis :List of 10
## ..$ call : language basisfd(type = type, rangeval = rangeval, nbasis = nbasis, params = params, dropind = dropind, quadvals = qu| __truncated__
## ..$ type : chr "fourier"
## ..$ rangeval : num [1:2] 0 1
## ..$ nbasis : num 3
## ..$ params : num 1
## ..$ dropind : num(0)
## ..$ quadvals : NULL
## ..$ values : list()
## ..$ basisvalues: list()
## ..$ names : chr [1:3] "const" "sin" "cos"
## ..- attr(*, "class")= chr "basisfd"
## $ fdnames:List of 3
## ..$ time : chr [1:20] "0.025" "0.075" "0.125" "0.175" ...
## ..$ reps : chr [1:39] "boy1" "boy2" "boy3" "boy4" ...
## ..$ values: chr [1:2] "Hip Angle" "Knee Angle"
## - attr(*, "class")= chr "fd"
str(gaitVar.fd3 <- var.fd(gaitfd3))
## List of 4
## $ coefs : num [1:3, 1:3, 1, 1:3] 30.5 1.05 5.13 1.05 5.2 ...
## $ sbasis :List of 10
## ..$ call : language basisfd(type = type, rangeval = rangeval, nbasis = nbasis, params = params, dropind = dropind, quadvals = qu| __truncated__
## ..$ type : chr "fourier"
## ..$ rangeval : num [1:2] 0 1
## ..$ nbasis : num 3
## ..$ params : num 1
## ..$ dropind : num(0)
## ..$ quadvals : NULL
## ..$ values : list()
## ..$ basisvalues: list()
## ..$ names : chr [1:3] "const" "sin" "cos"
## ..- attr(*, "class")= chr "basisfd"
## $ tbasis :List of 10
## ..$ call : language basisfd(type = type, rangeval = rangeval, nbasis = nbasis, params = params, dropind = dropind, quadvals = qu| __truncated__
## ..$ type : chr "fourier"
## ..$ rangeval : num [1:2] 0 1
## ..$ nbasis : num 3
## ..$ params : num 1
## ..$ dropind : num(0)
## ..$ quadvals : NULL
## ..$ values : list()
## ..$ basisvalues: list()
## ..$ names : chr [1:3] "const" "sin" "cos"
## ..- attr(*, "class")= chr "basisfd"
## $ bifdnames:List of 3
## ..$ time : chr [1:20] "0.025" "0.075" "0.125" "0.175" ...
## ..$ reps : chr [1:39] "boy1" "boy2" "boy3" "boy4" ...
## ..$ values:List of 3
## .. ..$ : chr "Hip Angle vs Hip Angle"
## .. ..$ : chr "Knee Angle vs Hip Angle"
## .. ..$ : chr "Knee Angle vs Knee Angle"
## - attr(*, "class")= chr "bifd"
# Check the answers with manual computations
all.equal(var(t(gaitfd3$coefs[,,1])), gaitVar.fd3$coefs[,,,1])
## [1] TRUE
# TRUE
all.equal(var(t(gaitfd3$coefs[,,2])), gaitVar.fd3$coefs[,,,3])
## [1] TRUE
# TRUE
all.equal(var(t(gaitfd3$coefs[,,2]), t(gaitfd3$coefs[,,1])),
gaitVar.fd3$coefs[,,,2])
## [1] TRUE
# TRUE
# NOTE:
dimnames(gaitVar.fd3$coefs)[[4]]
## NULL
# [1] Hip-Hip
# [2] Knee-Hip
# [3] Knee-Knee
# If [2] were "Hip-Knee", then
# gaitVar.fd3$coefs[,,,2] would match
#var(t(gaitfd3$coefs[,,1]), t(gaitfd3$coefs[,,2]))
# *** It does NOT. Instead, it matches:
#var(t(gaitfd3$coefs[,,2]), t(gaitfd3$coefs[,,1])),
##
## The following produces contour and perspective plots
##
# Evaluate at a 53 by 53 grid for plotting
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 )
str(tempvarbifd <- var.fd(daytempfd))
## List of 4
## $ coefs : num [1:65, 1:65] 13474.9 1951.3 4700.7 43.4 -455.9 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:65] "const" "sin1" "cos1" "sin2" ...
## .. ..$ : chr [1:65] "const" "sin1" "cos1" "sin2" ...
## $ sbasis :List of 10
## ..$ call : language basisfd(type = type, rangeval = rangeval, nbasis = nbasis, params = params, dropind = dropind, quadvals = qu| __truncated__
## ..$ type : chr "fourier"
## ..$ rangeval : num [1:2] 0 365
## ..$ nbasis : num 65
## ..$ params : num 365
## ..$ dropind : num(0)
## ..$ quadvals : NULL
## ..$ values : list()
## ..$ basisvalues: list()
## ..$ names : chr [1:65] "const" "sin1" "cos1" "sin2" ...
## ..- attr(*, "class")= chr "basisfd"
## $ tbasis :List of 10
## ..$ call : language basisfd(type = type, rangeval = rangeval, nbasis = nbasis, params = params, dropind = dropind, quadvals = qu| __truncated__
## ..$ type : chr "fourier"
## ..$ rangeval : num [1:2] 0 365
## ..$ nbasis : num 65
## ..$ params : num 365
## ..$ dropind : num(0)
## ..$ quadvals : NULL
## ..$ values : list()
## ..$ basisvalues: list()
## ..$ names : chr [1:65] "const" "sin1" "cos1" "sin2" ...
## ..- attr(*, "class")= chr "basisfd"
## $ bifdnames:List of 2
## ..$ : chr [1:65] "const" "sin1" "cos1" "sin2" ...
## ..$ : chr [1:65] "const" "sin1" "cos1" "sin2" ...
## - attr(*, "class")= chr "bifd"
str(tempvarmat <- eval.bifd(weeks,weeks,tempvarbifd))
## num [1:53, 1:53] 79.7 86.2 86.5 82.3 90.1 ...
# dim(tempvarmat)= c(53, 53)
op <- par(mfrow=c(1,2), pty="s")
#contour(tempvarmat, xlab="Days", ylab="Days")
contour(weeks, weeks, tempvarmat,
xlab="Daily Average Temperature",
ylab="Daily Average Temperature",
main=paste("Variance 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, tempvarmat,
xlab="Days", ylab="Days", zlab="Covariance")
mtext("Temperature Covariance", line=-4, outer=TRUE)
par(op)
##
## The simplest basis currently available with this function:
##
bspl1.1 <- create.bspline.basis(norder=1)
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)
## [1] TRUE
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)
##
## irregular
##
Time <- c(1:20, 41:60)
Birreg <- create.bspline.irregular(Time)
plot(Birreg)
# check
bks <- quantile(Time, seq(0, 1, length=4))
Bspi <- create.bspline.basis(c(1, 60), nbasis=round(sqrt(40)),
breaks=bks)
all.equal(Birreg, Bspi)
## [1] TRUE
##
## 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)
daybasis65 <- create.fourier.basis(c(0, 365), 65)
daytempfd <- with(CanadianWeather, smooth.basisPar(day.5,
dailyAv[,,"Temperature.C"], daybasis65)$fd )
plot(daytempfd, axes=FALSE)
## [1] "done"
axisIntervals(1)
# axisIntervals by default uses
# monthBegin.5, monthEnd.5, monthMid, and month.abb
axis(2)
# 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(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
par(mfrow=c(2,1))
plot(accmfd)
## [1] "done"
plot(accffd)
## [1] "done"
# 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(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
plot(hgtffd)
## [1] "done"
# 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)
# parameter object for coordinates
fdaPar = fdPar(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) )
# do the second order pda
pdaList = pda.fd(xfdlist, bwtlist)
# plot the results
eigres = eigen.pda(pdaList)
# 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(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
persp(tempvarmat)
##
## 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)
y=mu+t(L)%*%e
#functional boxplot
fbplot(y,method='MBD',ylim=c(-11,15))
## $depth
## [1] 0.3708844 0.3536327 0.3788844 0.3129252 0.3505850 0.3281633 0.4835374
## [8] 0.4718912 0.1672925 0.3802993 0.3998367 0.1841088 0.3810612 0.2517007
## [15] 0.3997279 0.2059320 0.2524626 0.4866395 0.3638095 0.3436190 0.5140680
## [22] 0.1787211 0.4948027 0.3160272 0.3374150 0.5056871 0.4408707 0.4456599
## [29] 0.2877279 0.3626122 0.4875646 0.4189932 0.4644898 0.3814422 0.5047075
## [36] 0.4143673 0.4075102 0.4898503 0.3780680 0.3903129 0.4333605 0.2010340
## [43] 0.2080000 0.3242449 0.2923537 0.3305034 0.4368435 0.1927619 0.3497143
## [50] 0.1432925
##
## $outpoint
## integer(0)
##
## $medcurve
## [1] 21
# The same using boxplot.fd
boxplot.fd(y, method='MBD', ylim=c(-11, 15))
## $depth
## [1] 0.3708844 0.3536327 0.3788844 0.3129252 0.3505850 0.3281633 0.4835374
## [8] 0.4718912 0.1672925 0.3802993 0.3998367 0.1841088 0.3810612 0.2517007
## [15] 0.3997279 0.2059320 0.2524626 0.4866395 0.3638095 0.3436190 0.5140680
## [22] 0.1787211 0.4948027 0.3160272 0.3374150 0.5056871 0.4408707 0.4456599
## [29] 0.2877279 0.3626122 0.4875646 0.4189932 0.4644898 0.3814422 0.5047075
## [36] 0.4143673 0.4075102 0.4898503 0.3780680 0.3903129 0.4333605 0.2010340
## [43] 0.2080000 0.3242449 0.2923537 0.3305034 0.4368435 0.1927619 0.3497143
## [50] 0.1432925
##
## $outpoint
## integer(0)
##
## $medcurve
## [1] 21
# same with default ylim
boxplot.fd(y)
## $depth
## [1] 0.3708844 0.3536327 0.3788844 0.3129252 0.3505850 0.3281633 0.4835374
## [8] 0.4718912 0.1672925 0.3802993 0.3998367 0.1841088 0.3810612 0.2517007
## [15] 0.3997279 0.2059320 0.2524626 0.4866395 0.3638095 0.3436190 0.5140680
## [22] 0.1787211 0.4948027 0.3160272 0.3374150 0.5056871 0.4408707 0.4456599
## [29] 0.2877279 0.3626122 0.4875646 0.4189932 0.4644898 0.3814422 0.5047075
## [36] 0.4143673 0.4075102 0.4898503 0.3780680 0.3903129 0.4333605 0.2010340
## [43] 0.2080000 0.3242449 0.2923537 0.3305034 0.4368435 0.1927619 0.3497143
## [50] 0.1432925
##
## $outpoint
## integer(0)
##
## $medcurve
## [1] 21
##
## 2. as an fd object
##
Y <- Data2fd(y)
## 'y' is missing, using 'argvals'
## 'argvals' is missing; using seq( 0 , 1 , length= 30 )
boxplot(Y)
## $depth
## [1] 0.3806102 0.3584643 0.3710244 0.3092908 0.3604203 0.3217701 0.4924874
## [8] 0.4914367 0.1287452 0.3695211 0.4299131 0.1797777 0.3917317 0.2499010
## [15] 0.4012528 0.1632249 0.2379551 0.4982421 0.3563306 0.3579794 0.5236371
## [22] 0.1621903 0.4879289 0.3204930 0.3493312 0.5231198 0.4493756 0.4663973
## [29] 0.2861427 0.3720752 0.4983067 0.4109517 0.4702445 0.3875611 0.5257870
## [36] 0.4336310 0.4216367 0.4975470 0.3806749 0.4120186 0.4415033 0.1766417
## [43] 0.1742332 0.2955668 0.2868539 0.3297393 0.4358780 0.1698848 0.3536149
## [50] 0.1069549
##
## $outpoint
## integer(0)
##
## $medcurve
## [1] 35
##
## 3. as an fdPar object
##
Ypar <- fdPar(Y)
boxplot(Ypar)
## $depth
## [1] 0.3806102 0.3584643 0.3710244 0.3092908 0.3604203 0.3217701 0.4924874
## [8] 0.4914367 0.1287452 0.3695211 0.4299131 0.1797777 0.3917317 0.2499010
## [15] 0.4012528 0.1632249 0.2379551 0.4982421 0.3563306 0.3579794 0.5236371
## [22] 0.1621903 0.4879289 0.3204930 0.3493312 0.5231198 0.4493756 0.4663973
## [29] 0.2861427 0.3720752 0.4983067 0.4109517 0.4702445 0.3875611 0.5257870
## [36] 0.4336310 0.4216367 0.4975470 0.3806749 0.4120186 0.4415033 0.1766417
## [43] 0.1742332 0.2955668 0.2868539 0.3297393 0.4358780 0.1698848 0.3536149
## [50] 0.1069549
##
## $outpoint
## integer(0)
##
## $medcurve
## [1] 35
##
## 4. Smoothed version
##
Ysmooth <- smooth.fdPar(Y)
boxplot(Ysmooth)
## $depth
## rep1 rep2 rep3 rep4 rep5 rep6 rep7
## 0.3806102 0.3584643 0.3710244 0.3092908 0.3604203 0.3217701 0.4924874
## rep8 rep9 rep10 rep11 rep12 rep13 rep14
## 0.4914367 0.1287452 0.3695211 0.4299131 0.1797777 0.3917317 0.2499010
## rep15 rep16 rep17 rep18 rep19 rep20 rep21
## 0.4012528 0.1632249 0.2379551 0.4982421 0.3563306 0.3579794 0.5236371
## rep22 rep23 rep24 rep25 rep26 rep27 rep28
## 0.1621903 0.4879289 0.3204930 0.3493312 0.5231198 0.4493756 0.4663973
## rep29 rep30 rep31 rep32 rep33 rep34 rep35
## 0.2861427 0.3720752 0.4983067 0.4109517 0.4702445 0.3875611 0.5257870
## rep36 rep37 rep38 rep39 rep40 rep41 rep42
## 0.4336310 0.4216367 0.4975470 0.3806749 0.4120186 0.4415033 0.1766417
## rep43 rep44 rep45 rep46 rep47 rep48 rep49
## 0.1742332 0.2955668 0.2868539 0.3297393 0.4358780 0.1698848 0.3536149
## rep50
## 0.1069549
##
## $outpoint
## integer(0)
##
## $medcurve
## rep35
## 35
##
## 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))
## $depth
## [1] 0.05175510 0.42236735 0.40310204 0.38394558 0.11412245 0.44440816
## [7] 0.49768707 0.36914286 0.21855782 0.47363265 0.33763265 0.07591837
## [13] 0.41931973 0.48772789 0.48500680 0.27689796 0.31521088 0.06742857
## [19] 0.39929252 0.48114286 0.47613605 0.45365986 0.35238095 0.47107483
## [25] 0.29970068 0.21121088 0.14416327 0.47493878 0.41687075 0.48348299
## [31] 0.18868027 0.33229932 0.43085714 0.40761905 0.47951020 0.46666667
## [37] 0.49779592 0.40734694 0.36288435 0.42655782 0.29736054 0.46960544
## [43] 0.29648980 0.34274830 0.05436735 0.50421769 0.47428571 0.39444898
## [49] 0.42095238 0.23738776
##
## $outpoint
## [1] 1 5 12 18 27 45
##
## $medcurve
## [1] 46
##
## 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)
## Warning in create.fourier.basis(c(0, 365), smallnbasis): nbasis must be an
## odd integer; is 40; will be increased by 1
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=10)
##
## 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
# cut this example to reduce test time on CRAN
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
growfdPar <- fdPar(hgtbasis, 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(hgtbasis,2,1e-6),fdPar(hgtbasis,2,1e-6))
# Should use nperm = 200 or so,
# but use 10 to save test time
Fres = Fperm.fd(hgtfd,xfdlist,betalist,nperm=10)
}
##
## 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(lipbasis, lipLfdobj, lipLambda)
lipfd <- smooth.basis(liptime, lip, lipfdPar)$fd
names(lipfd$fdnames) = c("Normalized time", "Replications", "mm")
lipmeanfd <- mean.fd(lipfd)
plot(lipmeanfd)
## [1] "done"
##
## 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))
## List of 3
## $ coefs : num [1:51, 1, 1:3] 35.87 -70.07 -194.03 -3.97 -8.4 ...
## ..- attr(*, "dimnames")=List of 3
## .. ..$ : NULL
## .. ..$ : chr "mean"
## .. ..$ : chr [1:3] "mean Temperature.C" "mean Precipitation.mm" "mean log10precip"
## $ basis :List of 10
## ..$ call : language basisfd(type = type, rangeval = rangeval, nbasis = nbasis, params = params, dropind = dropind, quadvals = qu| __truncated__
## ..$ type : chr "fourier"
## ..$ rangeval : num [1:2] 0 365
## ..$ nbasis : num 51
## ..$ params : num 365
## ..$ dropind : num(0)
## ..$ quadvals : NULL
## ..$ values : list()
## ..$ basisvalues: list()
## ..$ names : chr [1:51] "const" "sin1" "cos1" "sin2" ...
## ..- attr(*, "class")= chr "basisfd"
## $ fdnames:List of 3
## ..$ time : chr [1:365] "jan01" "jan02" "jan03" "jan04" ...
## ..$ reps : chr "mean"
## ..$ values: chr [1:3] "mean Temperature.C" "mean Precipitation.mm" "mean log10precip"
## - attr(*, "class")= chr "fd"