# pda.fd: Principal Differential Analysis In fda: Functional Data Analysis

## Description

Principal differential analysis (PDA) estimates a system of n linear differential equations that define functions that fit the data and their derivatives. There is an equation in the system for each variable.

Each equation has on its right side the highest order derivative that is used, and the order of this derivative, m_j, j=1,...,n can vary over equations.

On the left side of equation is a linear combination of all the variables and all the derivatives of these variables up to order one less than the order m_j of the highest derivative.

In addition, the right side may contain linear combinations of forcing functions as well, with the number of forcing functions varying over equations.

The linear combinations are defined by weighting functions multiplying each variable, derivative, and forcing function in the equation. These weighting functions may be constant or vary over time. They are each represented by a functional parameter object, specifying a basis for an expansion of a coefficient, a linear differential operator for smoothing purposes, a smoothing parameter value, and a logical variable indicating whether the function is to be estimated, or kept fixed.

## Usage

 ```1 2``` ```pda.fd(xfdlist, bwtlist=NULL, awtlist=NULL, ufdlist=NULL, nfine=501, returnMatrix=FALSE) ```

## Arguments

 `xfdlist` a list whose members are functional data objects representing each variable in the system of differential equations. Each of these objects contain one or more curves to be represented by the corresponding differential equation. The length of the list is equal to the number of differential equations. The number N of replications must be the same for each member functional data object. `bwtlist` this argument contains the weight coefficients that multiply, in the right side of each equation, all the variables in the system, and all their derivatives, where the derivatives are used up to one less than the order of the variable. This argument has, in general, a three-level structure, defined by a three-level hierarchy of list objects. At the top level, the argument is a single list of length equal to the number of variables. Each component of this list is itself a list At the second level, each component of the top level list is itself a list, also of length equal to the number of variables. At the third and bottom level, each component of a second level list is a list of length equal to the number of orders of derivatives appearing on the right side of the equation, including the variable itself, a derivative of order 0. If m indicates the order of the equation, that is the order of the derivative on the left side, then this list is length m. The components in the third level lists are functional parameter objects defining estimates for weight functions. For a first order equation, for example, m = 1 and the single component in each list contains a weight function for the variable. Since each equation has a term involving each variable in the system, a system of first order equations will have n^2 at the third level of this structure. There MUST be a component for each weight function, even if the corresponding term does not appear in the equation. In the case of a missing term, the corresponding component can be `NULL`, and it will be treated as a coefficient fixed at 0. However, in the case of a single differential equation, `bwtlist` can be given a simpler structure, since in this case only m coefficients are required. Therefore, for a single equation, `bwtlist` can be a list of length m with each component containing a functional parameter object for the corresponding derivative. `awtlist` a two-level list containing weight functions for forcing functions. In addition to terms in each of the equations involving terms corresponding to each derivative of each variable in the system, each equation can also have a contribution from one or more exogenous variables, often called forcing functions. This argument defines the weights multiplying these forcing functions, and is a list of length n, the number of variables. Each component of this is is in turn a list, each component of which contains a functional parameter object defining a weight function for a forcing function. If there are no forcing functions for an equation, this list can be `NULL`. If none of the equations involve forcing functions, `awtlist` can be `NULL`, which is its default value if it is not in the argument list. `ufdlist` a two-level list containing forcing functions. This list structure is identical to that for `awtlist`, the only difference being that the components in the second level contain functional data objects for the forcing functions, rather than functional parameter objects. `nfine` a number of values for a fine mesh. The estimation of the differential equation involves discrete numerical quadrature estimates of integrals, and these require that functions be evaluated at a fine mesh of values of the argument. This argument defines the number to use. The default value of 501 is reset to five times the largest number of basis functions used to represent any variable in the system, if this number is larger. `returnMatrix` logical: If TRUE, a two-dimensional is returned using a special class from the Matrix package.

## Value

an object of class `pda.fd`, being a list with the following components:

 `bwtlist` a list array of the same dimensions as the corresponding argument, containing the estimated or fixed weight functions defining the system of linear differential equations. `resfdlist` a list of length equal to the number of variables or equations. Each members is a functional data object giving the residual functions or forcing functions defined as the left side of the equation (the derivative of order m of a variable) minus the linear fit on the right side. The number of replicates for each residual functional data object is the same as that for the variables. `awtlist` a list of the same dimensions as the corresponding argument. Each member is an estimated or fixed weighting function for a forcing function.

`pca.fd`, `cca.fd`
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231``` ```#See analyses of daily weather data for examples. ## ## set up objects for examples ## # constant basis for estimating weight functions cbasis = create.constant.basis(c(0,1)) # monomial basis: {1,t} for estimating weight functions mbasis = create.monomial.basis(c(0,1),2) # quartic spline basis with 54 basis functions for # defining functions to be analyzed xbasis = create.bspline.basis(c(0,1),24,5) # set up functional parameter objects for weight bases cfd0 = fd(0,cbasis) cfdPar = fdPar(cfd0) mfd0 = fd(matrix(0,2,1),mbasis) mfdPar = fdPar(mfd0) # sampling points over [0,1] tvec = seq(0,1,len=101) ## ## Example 1: a single first order constant coefficient unforced equation ## Dx = -4*x for x(t) = exp(-4t) beta = 4 xvec = exp(-beta*tvec) xfd = smooth.basis(tvec, xvec, xbasis)\$fd xfdlist = list(xfd) bwtlist = list(cfdPar) # perform the principal differential analysis result = pda.fd(xfdlist, bwtlist) # display weight coefficient for variable bwtlistout = result\$bwtlist bwtfd = bwtlistout[[1]]\$fd par(mfrow=c(1,1)) plot(bwtfd) title("Weight coefficient for variable") print(round(bwtfd\$coefs,3)) # display residual functions reslist = result\$resfdlist plot(reslist[[1]]) title("Residual function") ## ## Example 2: a single first order varying coefficient unforced equation ## Dx(t) = -t*x(t) or x(t) = exp(-t^2/2) bvec = tvec xvec = exp(-tvec^2/2) xfd = smooth.basis(tvec, xvec, xbasis)\$fd xfdlist = list(xfd) bwtlist = list(mfdPar) # perform the principal differential analysis result = pda.fd(xfdlist, bwtlist) # display weight coefficient for variable bwtlistout = result\$bwtlist bwtfd = bwtlistout[[1]]\$fd par(mfrow=c(1,1)) plot(bwtfd) title("Weight coefficient for variable") print(round(bwtfd\$coefs,3)) # display residual function reslist = result\$resfdlist plot(reslist[[1]]) title("Residual function") ## ## Example 3: a single second order constant coefficient unforced equation ## Dx(t) = -(2*pi)^2*x(t) or x(t) = sin(2*pi*t) ## xvec = sin(2*pi*tvec) xfd = smooth.basis(tvec, xvec, xbasis)\$fd xfdlist = list(xfd) bwtlist = list(cfdPar,cfdPar) # perform the principal differential analysis result = pda.fd(xfdlist, bwtlist) # display weight coefficients bwtlistout = result\$bwtlist bwtfd1 = bwtlistout[[1]]\$fd bwtfd2 = bwtlistout[[2]]\$fd par(mfrow=c(2,1)) plot(bwtfd1) title("Weight coefficient for variable") plot(bwtfd2) title("Weight coefficient for derivative of variable") print(round(c(bwtfd1\$coefs, bwtfd2\$coefs),3)) print(bwtfd2\$coefs) # display residual function reslist = result\$resfdlist par(mfrow=c(1,1)) plot(reslist[[1]]) title("Residual function") ## ## Example 4: two first order constant coefficient unforced equations ## Dx1(t) = x2(t) and Dx2(t) = -x1(t) ## equivalent to x1(t) = sin(2*pi*t) ## xvec1 = sin(2*pi*tvec) xvec2 = 2*pi*cos(2*pi*tvec) xfd1 = smooth.basis(tvec, xvec1, xbasis)\$fd xfd2 = smooth.basis(tvec, xvec2, xbasis)\$fd xfdlist = list(xfd1,xfd2) bwtlist = list( list( list(cfdPar), list(cfdPar) ), list( list(cfdPar), list(cfdPar) ) ) # perform the principal differential analysis result = pda.fd(xfdlist, bwtlist) # display weight coefficients bwtlistout = result\$bwtlist bwtfd11 = bwtlistout[[1]][[1]][[1]]\$fd bwtfd21 = bwtlistout[[2]][[1]][[1]]\$fd bwtfd12 = bwtlistout[[1]][[2]][[1]]\$fd bwtfd22 = bwtlistout[[2]][[2]][[1]]\$fd par(mfrow=c(2,2)) plot(bwtfd11) title("Weight for variable 1 in equation 1") plot(bwtfd21) title("Weight for variable 2 in equation 1") plot(bwtfd12) title("Weight for variable 1 in equation 2") plot(bwtfd22) title("Weight for variable 2 in equation 2") print(round(bwtfd11\$coefs,3)) print(round(bwtfd21\$coefs,3)) print(round(bwtfd12\$coefs,3)) print(round(bwtfd22\$coefs,3)) # display residual functions reslist = result\$resfdlist par(mfrow=c(2,1)) plot(reslist[[1]]) title("Residual function for variable 1") plot(reslist[[2]]) title("Residual function for variable 2") ## ## Example 5: a single first order constant coefficient equation ## Dx = -4*x for x(t) = exp(-4t) forced by u(t) = 2 ## beta = 4 alpha = 2 xvec0 = exp(-beta*tvec) intv = (exp(beta*tvec) - 1)/beta xvec = xvec0*(1 + alpha*intv) xfd = smooth.basis(tvec, xvec, xbasis)\$fd xfdlist = list(xfd) bwtlist = list(cfdPar) awtlist = list(cfdPar) ufdlist = list(fd(1,cbasis)) # perform the principal differential analysis result = pda.fd(xfdlist, bwtlist, awtlist, ufdlist) # display weight coefficients bwtlistout = result\$bwtlist bwtfd = bwtlistout[[1]]\$fd awtlistout = result\$awtlist awtfd = awtlistout[[1]]\$fd par(mfrow=c(2,1)) plot(bwtfd) title("Weight for variable") plot(awtfd) title("Weight for forcing function") # display residual function reslist = result\$resfdlist par(mfrow=c(1,1)) plot(reslist[[1]], ylab="residual") title("Residual function") ## ## Example 6: two first order constant coefficient equations ## Dx = -4*x for x(t) = exp(-4t) forced by u(t) = 2 ## Dx = -4*t*x for x(t) = exp(-4t^2/2) forced by u(t) = -1 ## beta = 4 xvec10 = exp(-beta*tvec) alpha1 = 2 alpha2 = -1 xvec1 = xvec0 + alpha1*(1-xvec10)/beta xvec20 = exp(-beta*tvec^2/2) vvec = exp(beta*tvec^2/2); intv = 0.01*(cumsum(vvec) - 0.5*vvec) xvec2 = xvec20*(1 + alpha2*intv) xfd1 = smooth.basis(tvec, xvec1, xbasis)\$fd xfd2 = smooth.basis(tvec, xvec2, xbasis)\$fd xfdlist = list(xfd1, xfd2) bwtlist = list( list( list(cfdPar), list(cfdPar) ), list( list(cfdPar), list(mfdPar) ) ) awtlist = list(list(cfdPar), list(cfdPar)) ufdlist = list(list(fd(1,cbasis)), list(fd(1,cbasis))) # perform the principal differential analysis result = pda.fd(xfdlist, bwtlist, awtlist, ufdlist) # display weight functions for variables bwtlistout = result\$bwtlist bwtfd11 = bwtlistout[[1]][[1]][[1]]\$fd bwtfd21 = bwtlistout[[2]][[1]][[1]]\$fd bwtfd12 = bwtlistout[[1]][[2]][[1]]\$fd bwtfd22 = bwtlistout[[2]][[2]][[1]]\$fd par(mfrow=c(2,2)) plot(bwtfd11) title("weight on variable 1 in equation 1") plot(bwtfd21) title("weight on variable 2 in equation 1") plot(bwtfd12) title("weight on variable 1 in equation 2") plot(bwtfd22) title("weight on variable 2 in equation 2") print(round(bwtfd11\$coefs,3)) print(round(bwtfd21\$coefs,3)) print(round(bwtfd12\$coefs,3)) print(round(bwtfd22\$coefs,3)) # display weight functions for forcing functions awtlistout = result\$awtlist awtfd1 = awtlistout[[1]][[1]]\$fd awtfd2 = awtlistout[[2]][[1]]\$fd par(mfrow=c(2,1)) plot(awtfd1) title("weight on forcing function in equation 1") plot(awtfd2) title("weight on forcing function in equation 2") # display residual functions reslist = result\$resfdlist par(mfrow=c(2,1)) plot(reslist[[1]]) title("residual function for equation 1") plot(reslist[[2]]) title("residual function for equation 2") ```