Revision fd92910f232786c2531c1c74e02ece6b49a2f76f authored by J. O. Ramsay on 03 November 2009, 00:00:00 UTC, committed by Gabor Csardi on 03 November 2009, 00:00:00 UTC
1 parent 8ce7a35
fda-ch02.R
###
###
### Ramsey & Silverman (2006) Functional Data Analysis, 2nd ed. (Springer)
###
### DISCLAIMER:
###
### These script files are incomplete and where complete
### do not always reproduce exactly the images in the books.
### Reason for this include the following:
###
### (1) The authors were given permission to publish the results
### of thier analyses but not always the data.
### (2) These script files were produced starting in 2007
### not by the authors to the books but by volunteer(s)
### (initially Spencer Graves). The published analyses
### were conducted over more than a decade, and the details
### of exactly how the published analyses were performed were
### not in some cases still available.
### (3) The software and capabilities of computers available
### have improved over time. In some cases, this produces
### a better fit today than at the times the published analyses
### were prepared.
###
### It is hoped that these script files may still be useful,
### even with their gaps and even when they do not reproduce exactly
### the published analyses.
###
### Spencer Graves
###
###
### ch. 2. Tools for exploring functional data
###
###
library(fda)
##
## Section 2.1. Introduction
##
# No figures in this section.
##
## Section 2.2. Some notation
##
# No figures in this section.
##
## Section 2.3. Summary statistics for functional data
##
# pp. 22-23, Figure 2.1.
# Mean and standard deviation of pinch force
# after alignment / registration.
matplot(pinchtime, pinch, type="l")
# experiment with lambda to get a plot that looks
# reasonably smooth but not too smooth
pinch.fd <- smooth.basisPar(pinchtime, pinch)
pinch.fd <- smooth.basisPar(pinchtime, pinch, lambda=1e-6)
plot(pinch.fd$fd)
# The default lambda = 1/diff(pinchtime) produces a straight line
# lambda = 1e-8 undersmooths
# lambda = 1e-6 comes close to matching Figure 2.1
# ... but with discrepancies
# suggesting that the 'pinch' data used for Figure 2.1
# may have differed by more than just an offset of 2
# from the 'pinch' data in the 'fda' package
plot(mean(pinch.fd$fd), xlab="Time (sec.)", ylab="Force mean (N)")
abline(v=.1, lty="dotted")
plot(sd(pinch.fd$fd), ylim=c(0, 1), xlab="Time (sec.)",
ylab="Force Std. Dev. (N)")
abline(v=.1, lty="dotted")
# Reconsider after ch. 7?
# pp. 23, Figure 2.2. Perspective and contour plots of pinch force correlations
# complete after ch. 7
# pp. 24-25, Figure 2.3.
# Contour plots of correlation and cross correlations of gait data
#tempPrecCor4.6 <- cor.fd(weeks, tempSmooth4$fd, fdobj2=precSmooth6$fd)
# In \demo\gaitDemo.R, lambda=1e-11 was determined as optimal smoothing
# to minimize the generalized cross validation (GCV) criterion
# 'gaittime' is retained as the first of the gait dimnames
(gaittime <- as.numeric(dimnames(gait)[[1]]))
# Establish a full basis;
# 11 or 15 basis function might do as well or better,
# but I haven't tried that.
gaitbasis21 <- create.fourier.basis(nbasis=21)
# A finite Fourier series must satisfy a harmonic
# differential equation, as described in the book.
# Give that a name.
harmaccelLfd01 <- vec2Lfd(c(0, 0, (2*pi)^2, 0))
# lambda = 1e-11 is the roughness penalty
# found to minimize the generalized cross validation criterion;
# see '~R\library\demo\gaitDemo.R'.
gaitfd <- smooth.basisPar(gaittime, gait,
gaitbasis21, Lfdobj=harmaccelLfd01, lambda=1e-11)$fd
gaitfd$fdnames
# Change some of the names
names(gaitfd$fdnames) = c("Normalized time", "Child", "Angle")
gaitfd$fdnames[[3]] = c("Hip", "Knee")
# Establish a finer grid so the plot appears smoother
(gaitTimeSm <- seq(0, 1, length=41))
gaitCor <- cor.fd(gaitTimeSm, gaitfd)
gaitCorLbls <- (-5:5)/5
contour(gaitTimeSm, gaitTimeSm, gaitCor[,,1,1], cex=1.2,
levels=gaitCorLbls, xlab="Time in Knee Angle Cycle",
ylab="Time in Knee Angle Cycle")
x01. <- seq(0, 1, .05)
text(x01., x01., 1)
title("Knee - Knee")
# Smoothed more than Figure 2.3,
# which seems too rough
contour(gaitTimeSm, gaitTimeSm, gaitCor[,,1,2], cex=1.2,
levels=gaitCorLbls, ylab="Time in Knee Angle Cycle",
xlab="Time in Hip Angle Cycle")
title("Hip - Knee")
contour(gaitTimeSm, gaitTimeSm, gaitCor[,,1,2], cex=1.2,
levels=gaitCorLbls, xlab="Time in Knee Angle Cycle",
ylab="Time in Hip Angle Cycle")
title("Knee - Hip")
contour(gaitTimeSm, gaitTimeSm, gaitCor[,,1,3], cex=1.2,
levels=gaitCorLbls, xlab="Time in Hip Angle Cycle",
ylab="Time in Hip Angle Cycle")
text(x01., x01., 1)
title("Hip - Hip")
##
## pp. 25-26, Figure 2.4. Correlations for Canadian temperature & precipitation
##
# from demo\CanadianWeatherDemo.R
#
# ----------------------------------------------------------------------
# Descriptive Statistics Functions
# ----------------------------------------------------------------------
# ------------- set up fourier basis ---------------------------
# It was decided that 65 basis functions captured enough of
# the detail in the temperature data: about one basis function
# per week for present purposes.
daybasis65 <- create.fourier.basis(c(0, 365), nbasis=65)
# ----------- set up the harmonic acceleration operator ----------
harmaccelLfd365 <- vec2Lfd(c(0,(2*pi/365)^2,0), c(0, 365))
# For other purposes, a saturated basis (365 basis functions) is used,
# where smoothing is defined by the GCV (Generalized Cross Validation)
# criterion.
# --------- create fd objects for temp. and prec. ---------------
qqnorm(CanadianWeather$dailyAv[,,"Temperature.C"], datax=TRUE)
# Short tailed distribution
# because it's a strong, deterministic annual cycle
# + relatively modest normal error
daytempfd <- with(CanadianWeather, data2fd(dailyAv[,,"Temperature.C"],
day.5, daybasis65, argnames=list("Day", "Station", "Deg C")) )
qqnorm(CanadianWeather$dailyAv[,,"Precipitation.mm"], datax=TRUE)
# Apparent lognormal distribution
qqnorm(CanadianWeather$dailyAv[,,"log10precip"], datax=TRUE)
# Looks more like the 'tempav' distribution.
# Too discrete in the lower tail,
# because the numbers were rounded off too grossly
# before computing the logarithms.
sum(CanadianWeather$dailyAv[,,"Precipitation.mm"]==0)
# 27 of 365 numbers are 0
quantile(CanadianWeather$dailyAv[,,"log10precip"])
# Use log10precip rather than Precipitation.mm directly.
# Smooth temparature
# CanadianWeatherDemo.R compares alternative smoothing levels
# with a full basis rather than daybasis65.
# The generalized cross validation criterion is minimized for
# lambda = 0.01. However, that ignores serial correlation and
# seems visually to undersmooth; lambda = 1e6 looks more realistic
# and will be used here, even though the reduction in the number
# of basis vectors from 365 to 65 provides some smoothing without
# a roughness penalty.
TempSm6 <- smooth.basisPar(argvals=day.5, y=CanadianWeather$dailyAv[,,1],
fdobj=daybasis65, Lfdobj=harmaccelLfd365, lambda=1e6)$fd
# The minimal GCV in CanadianWeatherDemo.R for log10precip
# was also lambda=1e6 and will be used here as it gives a reasonable
# physical image in the plots.
PrecipSm6 <- smooth.basisPar(argvals=day.5, y=CanadianWeather$dailyAv[,,3],
fdobj=daybasis65, Lfdobj=harmaccelLfd365, lambda=1e6)$fd
lvl.1 <- (-10:10)/10
contour(weeks, weeks, cor.fd(weeks, TempSm6), levels=lvl.1,
xlab="Average daily Temperature (C)",
ylab="Average daily Temperature (C)",
axes=FALSE, main="Temperature correlations")
axisIntervals(1)
axisIntervals(2)
contour(weeks, weeks, cor.fd(weeks, PrecipSm6), levels=lvl.1,
xlab="Average daily Temperature (C)",
ylab="Average daily Temperature (C)",
axes=FALSE, main="Correlations in log10(Precipitation(mm))")
axisIntervals(1)
axisIntervals(2)
contour(weeks, weeks, cor.fd(weeks, TempSm6, weeks, PrecipSm6),
levels=lvl.1, xlab="Average daily Temperature (C)",
ylab="Average daily Temperature (C)", axes=FALSE,
main="Correlations in log10(Precipitation(mm))")
axisIntervals(1)
axisIntervals(2)
# Places where January is warmer are quite likely to have more rain (r=.8).
# Places where June is warmer are only moderately to have more rain (r=.4).
# This is roughly similar to Figure 2.4.
# To get the exact image,
# we would need the right combination of the following:
# (1) The exact same basis set,
# which probably was create.fourier.basis,
# but might have been something else.
# (2) Exactly the same order 'nbasis'.
# (3) Exactly the same lambda and Lfdobj for smoothing.
# (4) Evalated at exactly the same points
# I don't know how to find that exact combination.
# However, I believe this is essentially equivalent
# to how that figure was created.
##
## Section 2.4. The anatomy of a function
##
# No figures in this section
##
## Section 2.5. Phase-plane plots of periodic effects
##
# Figure 2.5. Nondurable goods production in the US
plot(nondurables, log="y")
str(nondurables)
plot(log10(nondurables))
(log.nondur.lm <- lm(log(nondurables)~time(nondurables)))
(log10.nondur.lm <- lm(log10(nondurables)~time(nondurables)))
abline(nondur.lm, lty="dashed")
length(nondurables)
diff(range(time(nondurables)))
# 81
(nondurGrowth <- (nondurables[973]/nondurables[1]))
# 15.221 = (1+growthRate)^(81 years)
#log(1+growthRate) = log(15.221)/81 = 0.0336
# 1+growthRate = nondurGrowth^(1/81) = 1.0342
#1.016^81 = 3.62
#1.0342^81 = 15.24
log(nondurables[973]/nondurables[1])/81
# Figure 2.6. log nondurable goods 1964-1967
nonDur1964.67 <- window(nondurables, 1964, 1967)
plot(log10(nondur1964.1967), type="p", axes=FALSE, xlab="Year",
ylab=expression(paste(log[10], " nondurable goods index")) )
axis(2)
axis(1, 1964:1967)
axis(1, seq(1964, 1967, by=0.5), labels=FALSE)
# For theory on smoothing, see FDA, ch. 5
# For more on this example, see AFDA, ch. 3
# or \demo\goodindexDemo.R
durtimefine <- seq(1964, 1967, length=181)
#fit = eval.fd(durtimefine, lognondursmth);
logNondurSm1964.67 = eval.fd(durtimefine, logNondurSm$fd);
lines(durtimefine, logNondurSm1964.67)
abline(v=1965:1966, lty=2)
# Figure 2.7. Phase-plane plot for a simple harmonic function
sin. <- expression(sin(2*pi*x))
D.sin <- D(sin., "x")
D2.sin <- D(D.sin, "x")
with(data.frame(x=seq(0, 1, length=46)),
plot(eval(D.sin), eval(D2.sin), type="l",
xlim=c(-10, 10), ylim=c(-50, 50),
xlab="Velocity", ylab="Acceleration") )
pi.2 <- (2*pi)
#lines(x=c(-pi.2,pi.2), y=c(0,0), lty=3)
abline(h=0, lty="longdash")
pi.2.2 <- pi.2^2
lines(x=c(0,0), y=c(-pi.2.2, pi.2.2), lty="longdash")
text(c(0,0), c(-47, 47), rep("no kinetic, max potential", 2))
text(c(-8.5,8.5), c(0,0), rep("max kinetic\nno potential", 2))
# Figure 2.8. Phase-plane plot of nondurable goods index 1964
goodsbasis <- create.bspline.basis(rangeval=c(1919,2000),
nbasis=979, norder=8)
LfdobjNonDur <- int2Lfd(4)
library(zoo)
logNondurSm <- smooth.basisPar(argvals=index(nondurables),
y=log10(coredata(nondurables)), fdobj=goodsbasis,
Lfdobj=LfdobjNonDur, lambda=1e-11)
phaseplanePlot(1964, logNondurSm$fd)
![swh spinner](/static/img/swh-spinner.gif)
Computing file changes ...