# demo/Fdascript.R In Data2LD: Functional Data Analysis with Linear Differential Equations

#  The data are analyzed in the Ramsay and Silverman books, and consist of
#  of the script version of "fda" written by Jim Ramsay.  Because script
#  is cursive, there are two variables, "X" being the motion across the
#  front of the body, and "Y" being fore-and-aft motion.

#  The two second order linear differential equations analyzed here are:
#  $D^2X(t) = -b_X X(t) + a_X(t)$
#  $D^2Y(t) = -b_Y Y(t) + a_Y(t)$
#  where $b_X$ and $b_Y$ are scalar constants, called the "reaction speed"
#  for the coefficient.  The coefficient functions $a_X(t)$ and $a_Y(t)$
#  are step functions, each having 20 steps.  They can be thught of as
#  coefficient functions multiplying a dummy forcing function U(t) which
#  only takes the value one.
#  Thus, each equation there are 21 parameters to estimate from the data,
#  the coordinate's reaction speed and the 20 heights of the steps.
#  Although there are no shared terms in these two equations, so that
#  they could be estimated indendently, for the purposes of neatness and
#  clarity, we shall treat them as a system of two equations with a
#  total of 42 parameters to be estimated.
#  Each estimated step function is a spline function of order one.

#  --------------------------------------------------------------------
#        Load the fda script data and define the time values
#  --------------------------------------------------------------------

#  select out the two 1401 by 20 matrices of coordinate values

fdascriptX <- as.matrix(fdascript[,1,1])
fdascriptY <- as.matrix(fdascript[,1,2])

#  Define the observation times in 100ths of a second

centisec <- seq(0,230,len=1401)
fdarange <- range(centisec)

yList = vector("list", 2)
yList[[1]]$argvals <- centisec yList[[1]]$y       <- fdascriptX
yList[[2]]$argvals <- centisec yList[[2]]$y       <- fdascriptY

#  --------------------------------------------------------------------
#   Set up the basis object for representing the X and Y coordinate
#   functions
#  --------------------------------------------------------------------

#  The basis functions for each coordinate will be order 6 B-splines,
#  32 basis functions per second, 4 per 8 herz cycle
#  2.3 seconds 2.3*32=73.6.
#  Order 6 is used because we want to estimate a smooth second
#  derivative.
#  This implies norder + no. of interior knots <- 74 - 1 + 6 <- 79
#  basis functions.

nbasis   <- 79
norder   <-  6 #  order 6 for a smooth 2nd deriv.
fdabasis <- create.bspline.basis(fdarange, nbasis, norder)

XbasisList <- vector("list",2)
XbasisList[[1]] <- fdabasis
XbasisList[[2]] <- fdabasis

#  --------------------------------------------------------------------
#   Use the basis to estimate the curve functions and their
#   second derivatives
#  --------------------------------------------------------------------

ResultX <- smooth.basis(centisec, fdascriptX, fdabasis)
fdascriptfdX <- ResultX$fd D2fdascriptX <- eval.fd(centisec, fdascriptfdX, 2) ResultY <- smooth.basis(centisec, fdascriptY, fdabasis) fdascriptfdY <- ResultY$fd
D2fdascriptY <- eval.fd(centisec, fdascriptfdY, 2)

#  --------------------------------------------------------------------
#   Set up the basis objects for representing the coefficient
#   functions
#  --------------------------------------------------------------------

#  Define a constant basis and function over the 230 centisecond range

conbasis <- create.constant.basis(fdarange)
confd    <- fd(1,conbasis)
confdPar <- fdPar(confd)

#  Define the order one Bspline basis for the coefficient
#  of the forcing coefficients

nevent    <- 20
nAorder   <- 1
nAbasis   <- nevent
nAknots   <- nAbasis + 1
Aknots    <- seq(fdarange[1],fdarange[2],len=nAknots)
Abasisobj <- create.bspline.basis(fdarange,nAbasis,nAorder,Aknots)

#  --------------------------------------------------------------------
#  Set up single homogeneous terms in the homogeneous portion of the
#  equations  D^2x = -beta_x x and D^2y = -beta_y y.
#  Each term involves a scalar constant multiplier -1.
#  --------------------------------------------------------------------

Xterm <- make.Xterm(confdPar, 0.4, TRUE, variable=1, derivative=0, factor=-1)
Yterm <- make.Xterm(confdPar, 0.4, TRUE, variable=2, derivative=0, factor=-1)

#  Set up to list arrays to hold these term specifications

XListX <- vector("list",1)
XListX[[1]] <- Xterm

XListY <- vector("list",1)
XListY[[1]] <- Yterm

#  --------------------------------------------------------------------
#  Set up coefficients for the forcing terms
#  \alpha_x(t)*1 and \alpha_y(t)*1.
#  The forcing functions are the unit constant function.
#  --------------------------------------------------------------------

parvecA <- matrix(1,nAbasis,1)
FtermX <- make.Fterm(Abasisobj, parvecA, TRUE, Ufd=confd, factor=1)
FtermY <- make.Fterm(Abasisobj, parvecA, TRUE, Ufd=confd, factor=1)

#  Set up list arrays to hold forcing terms specs

FListX <- vector("list", 1)
FListX[[1]] <- FtermX

FListY <- vector("list", 1)
FListY[[1]] <- FtermY

#  --------------------------------------------------------------------
#  make variables for X and Y coordinates and system object fdaVariableList
#  --------------------------------------------------------------------

Xvariable <- make.Variable(name="X", order=2, XList=XListX, FList=FListX)
Yvariable <- make.Variable(name="Y", order=2, XList=XListY, FList=FListY)

#  List array for the whole system

fdaVariableList <- vector("list",2)
fdaVariableList[[1]] <- Xvariable
fdaVariableList[[2]] <- Yvariable

#  check the system specification for consistency,  This is a mandatory
#  step because objects are set up in the output list that are needed
#  in the analysis.  A summary is displayed.

checkfdaModel <- checkModel(XbasisList, fdaVariableList, TRUE)

fdaModel <- checkfdaModel$modelList nparam <- checkfdaModel$nparam

print(paste("Number of parameters =",nparam))

rhoVec <- matrix(0.5,2,1)

Data2LDResult <- Data2LD(yList, XbasisList, fdaModel, rhoVec, summary=TRUE)

MSE0      <- Data2LDResult$MSE # Mean squared error for fit to data DpMSE0 <- Data2LDResult$DpMSE      #  gradient with respect to parameter values
D2ppMSE0  <- Data2LDResult$D2ppMSE # Hessian matrix XfdList <- Data2LDResult$XfdParList
df        <- Data2LDResult$df gcv <- Data2LDResult$gcv
ISE       <- Data2LDResult$ISE Var.theta <- Data2LDResult$Var.theta
Rmat      <- Data2LDResult$Rmat Smat <- Data2LDResult$Smat
DRarray   <- Data2LDResult$DRarray DSmat <- Data2LDResult$DSmat

MSE0
DpMSE0

#  --------------------------------------------------------------------
#  Set up sequence of rho values to be used in the optimization
#  and the display level, iteration limit and convergence criterion
#  --------------------------------------------------------------------

gamvec <- 0:7
rhoMat <- t(matrix(exp(gamvec)/(1+exp(gamvec)),8,2))
nrho   <- dim(rhoMat)[2]

#  values controlling optimization

dbglev   <- 1     #  display level
iterlim  <- 50    #  maximum number of iterations
convrg   <- 1e-6  #  convergence criterion

#  --------------------------------------------------------------------
#  Set up arrays to hold results over rho and
#  Initialize coefficient list
#  --------------------------------------------------------------------

dfesave <- matrix(0,nrho,1)
gcvsave <- matrix(0,nrho,1)
MSEsave <- matrix(0,nrho,1)
ISEsave <- matrix(0,nrho,1)
thesave <- matrix(0,nrho,nparam)

fdaModel.opt <- fdaModel

#  --------------------------------------------------------------------
#  Step through rho values, optimizing at each step, and saving
#  results along the way.
#  --------------------------------------------------------------------

for (irho in 1:nrho) {
rhoVeci <- matrix(rhoMat[,irho],2,1)
print(paste("rho <- ",round(rhoVeci[1],5)))
Result <- Data2LD.opt(yList, XbasisList, fdaModel.opt, rhoVeci,
convrg, iterlim, dbglev)
theta.opti       <- Result$thetastore fdaModel.opt <- modelVec2List(fdaModel.opt, theta.opti) Data2LDResult <- Data2LD(yList, XbasisList, fdaModel.opt, rhoVeci, summary=TRUE) thesave[irho,] <- theta.opti dfesave[irho,1] <- Data2LDResult$df
gcvsave[irho,1]  <- Data2LDResult$gcv MSEsave[irho,1] <- Data2LDResult$MSE
ISEsave[irho,1]  <- Data2LDResult$ISE } # Each of the optimizations is started from the converged value for # previous optimization, and they all converge smoothly and rapidly. # If you want see how this can fail, change fdaModel.opt to just # fdaModel, so that each iteration starts from the initial parameter # values. You will see a catastrophic failure for the last value of # rho. # -------------------------------------------------------------------- # Display degrees of freedom and gcv values for each step # -------------------------------------------------------------------- print(" rho df gcv") print(round(cbind(rhoMat[1,], dfesave, gcvsave),4)) # plot flow across rho values for the two homogeneous coefficients matplot(rhoMat[1,], thesave[,1:2], type="b", pch="o",lwd=2, xlab = 'rho', ylabel = 'beta coefficients') legend(x=0.6, y=0.04, c("X", "Y"), lty=c(1,2), lwd=2, col=c(1,2)) # plot flow of MSE and ISE par(mfrow=c(2,1)) plot(rhoMat[1,], MSEsave, type="b", pch="o",lwd=2, xlab = "rho", ylab = "MSE") subplot(2,1,2) plot(rhoMat[1,], ISEsave, type="b", pch="o",lwd=2, xlab = "rho", ylab = "ISE") # -------------------------------------------------------------------- # Evaluate the fit for parameter values at the highest rho value # and display the numerical results # -------------------------------------------------------------------- irho <- 8 theta <- thesave[irho,] fdaModel <- modelVec2List(fdaModel, theta) rho <- rhoMat[,irho] Data2LDResult <- Data2LD(yList, XbasisList, fdaModel, rho, summary=TRUE) MSE <- Data2LDResult$MSE        #  Mean squared error for fit to data
DpMSE      <- Data2LDResult$DpMSE # gradient with respect to parameter values D2ppMSE <- Data2LDResult$D2ppMSE    #  Hessian matrix
XfdParList <- Data2LDResult$XfdParList # List of fdPar objects for variable values df <- Data2LDResult$df         #  Degrees of freedom for fit
gcv        <- Data2LDResult$gcv # Generalized cross-validation coefficient ISE <- Data2LDResult$ISE        #  Size of second term, integrated squared error
Var.theta  <- Data2LDResult$Var.theta # Estimate sampling variance for parameters print(paste("MSE = ", round(MSE,4))) print(paste("df = ", round(df, 4))) print(paste("gcv = ", round(gcv,4))) print(paste("X reaction speed = ",round(theta[1],3))) print(paste("Y reaction speed = ",round(theta[2],3))) # -------------------------------------------------------------------- # Display the graphical results # -------------------------------------------------------------------- fdascriptfdX <- Data2LDResult$XfdParList[[1]]$fd fdascriptfdY <- Data2LDResult$XfdParList[[2]]$fd fdascriptTimefine <- seq(fdarange[1],fdarange[2],len=201) fdascriptfineX <- eval.fd(fdascriptTimefine, fdascriptfdX) fdascriptfineY <- eval.fd(fdascriptTimefine, fdascriptfdY) fdascriptKnotsX <- eval.fd(Aknots, fdascriptfdX) fdascriptKnotsY <- eval.fd(Aknots, fdascriptfdY) fdascriptKnotsX <- approx(centisec, fdascriptX[,1], Aknots) fdascriptKnotsY <- approx(centisec, fdascriptY[,1], Aknots) # plot fit to the data as a script par(mfrow=c(1,1)) plot(fdascriptfineX, fdascriptfineY, type="l", xlab="X coordinate", ylab="Y coordinate") points(fdascriptX[,1], fdascriptY[,1], pch="o") points(fdascriptKnotsX$y, fdascriptKnotsY\$y, pch="X", col=2)

#  plot fit to the data in two-panels

par(mfrow=c(2,1))
plot(centisec, fdascriptX[,1], type="p",
xlab="Time (centisec)", ylab="X-coordinate")
lines(fdascriptTimefine, fdascriptfineX)

plot(centisec, fdascriptY[,1], type="p",
xlab="Time (centisec)", ylab="Y-coordinate")
lines(fdascriptTimefine, fdascriptfineY)

#  plot the step heights

fdascriptKnotsX    = eval.fd(Aknots, fdascriptfdX);
fdascriptKnotsY    = eval.fd(Aknots, fdascriptfdY);

par(mfrow=c(2,1))

plot( c(  0,  0), c(0,fdascriptKnotsX[1]), type="l", lwd=2,
xlim=c(0,230), ylim=c(-40,40), xlab = "", ylab = "X")
lines(c(230,230), c(fdascriptKnotsX[20],0), type="l", lwd=2)
lines(c(  0,230), c(0,0), type="l", lwd = 2, lty=2)
for (i in 1:20) {
lines(c(Aknots[i],Aknots[i+1]), c(fdascriptKnotsX[i],fdascriptKnotsX[i]),
type="l", lwd=2)
}
for (i in 1:19) {
lines(c(Aknots[i+1],Aknots[i+1]),
c(fdascriptKnotsX[i],fdascriptKnotsX[i+1]), type="l", lwd=2)
}

plot( c(  0,  0), c(0,fdascriptKnotsY[1]), type="l", lwd=2,
xlim=c(0,230), ylim=c(-40,40), xlab = "Time (centiseconds", ylab = "Y")
lines(c(230,230), c(fdascriptKnotsY[20],0), type="l", lwd=2)
lines(c(  0,230), c(0,0), type="l", lwd = 2, lty=2)
for (i in 1:20) {
lines(c(Aknots[i],Aknots[i+1]), c(fdascriptKnotsY[i],fdascriptKnotsY[i]),
type="l", lwd=2)
}
for (i in 1:19) {
lines(c(Aknots[i+1],Aknots[i+1]),
c(fdascriptKnotsY[i],fdascriptKnotsY[i+1]), type="l", lwd=2)
}

#  plot step heights as a script

indX <-  3:22;
indY <- 23:42;

StepX <- theta[indX]
StepY <- theta[indY]

plot(StepX, StepY, type="b",
xlab="X coordinate", ylab="Y coordinate")


## Try the Data2LD package in your browser

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

Data2LD documentation built on Aug. 6, 2020, 1:08 a.m.