Description Usage Arguments Value References See Also Examples
This function calls function Data2LD
in each iteration of a quasi-Newton
type optimization algorithm, allowing for a set of linear restrictions on the
parameter vector.
1 2 | Data2LD.opt(yList, XbasisList, modelList, rhoMat,
convcrit = 1e-4, iterlim = 20, dbglev = 1, parMap = diag(rep(1,nparam)))
|
yList |
A list of length NVAR. Each list contains in turn a struct object with fields:
|
XbasisList |
A list array of length d. Each member contains in turn a functional data object or a functional basis object. |
modelList |
A list each member of which specifies one of linear differential
equations in the system. These can be
constructed using function |
rhoMat |
A matrix with the number of rows equal to the number of values of the smoothing parameter per variable to be used in the optimization, and number of columns equal to the number of variables. Each entry in the matrix is a value rho in the interval [0,1) so that 0 <= rho < 1. For each variable ivar in the system of equations and each optimization iopt, the data are weighted by rho(iopt,ivar) and the roughness penalty by 1-rho(iopt,ivar). It is expected that the values of rho within each column will be in ascending order, and the estimated parameters for each row are passed along as initial values to be used for the optimization defined by the values of rho in the next row. |
convcrit |
A one- or two-part convergence criterion. The first part is for testing the convergence of the criterion values, and the second part, if present in a vector of length 2, is for testing the norm of the gradient. |
iterlim |
Maximum number of iterations allowed. |
dbglev |
Control the output on each iteration. If 0, no output is displayed. If 1, the criterion and gradient norm at each iteration. If 2, the stepsize, slope and criterion at each step in the line search. if 3, displays, in addition to 2, a plot of the criterion values and slope after each line search followed by a pause. |
parMap |
A rectangular matrix with number of rows equal to the
number of parameters to be estimated, and number of columns equal to the number of
parameters less the number of linear constraints on the estimated parameters.
The columns of |
A named list object containing these results of the analysis:\
theta.opt: A vector containing the optized parameter values.
modelList.opt: A list object defining the optimized linear equation
system. This may be displayed using function printModel
.
J. O. Ramsay and G. Hooker (2017) Dynamic Data Analysis. Springer.
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 | #
# Example 1: The refinery data
#
N <- 194
TimeData <- RefineryData[,1]
TrayData <- RefineryData[,2]
ValvData <- RefineryData[,3]
# Define the List array containing the tray data
TrayList <- list(argvals=RefineryData[,1], y=RefineryData[,2])
TrayDataList <- vector("list",1)
TrayDataList[[1]] <- TrayList
# construct the basis object for tray variable
# Order 5 spline basis with four knots at the 67
# to allow discontinuity in the first derivative
# and 15 knots between 67 and 193
Traynorder <- 5
Traybreaks <- c(0, rep(67,3), seq(67, 193, len=15))
Traynbasis <- 22
TrayBasis <- create.bspline.basis(c(0,193), Traynbasis, Traynorder,
Traybreaks)
# Set up the basis list for the tray variable
TrayBasisList <- vector("list",1)
TrayBasisList[[1]] <- TrayBasis
# Both alpha and beta coefficient functions will be constant.
# Define the constant basis
conbasis <- create.constant.basis(c(0,193))
betafdPar <- fdPar(conbasis)
alphafdPar <- fdPar(conbasis)
# Construct a step basis (order 1) and valve level
Valvbreaks <- c(0,67,193)
Valvnbasis <- 2
Valvnorder <- 1
Valvbasis <- create.bspline.basis(c(0,193), Valvnbasis, Valvnorder, Valvbreaks)
# smooth the valve data
Valvfd <- smooth.basis(TimeData, ValvData, Valvbasis)$fd
# Set up the model list for the tray variable
# Define single homogeneous term
# Xterm Fields: funobj parvec estimate variable deriv. factor
XTerm <- make.Xterm(betafdPar, 0.04, TRUE, 1, 0, -1)
XList <- vector("list", 1)
XList[[1]] <- XTerm
# Define the single forcing term
# Fterm Fields: funobj parvec estimate Ufd factor
FTerm <- make.Fterm(alphafdPar, 1.0, TRUE, Valvfd, 1)
FList <- vector("list", 1)
FList[[1]] <- FTerm
# Define the single differential equation in the model
TrayVariable <- make.Variable(XList=XList, FList=FList,
name="Tray level", order=1)
# Set up the model List array
TrayModelList <- vector("list",1)
TrayModelList[[1]] <- TrayVariable
# Run a check on TrayVariableList and make the modelList object
checkModelList <- checkModel(TrayBasisList, TrayModelList, summarywrd=TRUE)
TrayModelList <- checkModelList$model
nparam <- checkModelList$nparam
# Evaluate the fit to the data given the initial parameter estimates (0 and 0)
# This also initializes the four-way tensors so that they are not re-computed
# for subsquent analyses.
rhoVec <- 0.5
Data2LDList <- Data2LD(TrayDataList, TrayBasisList, TrayModelList, rhoVec)
MSE <- Data2LDList$MSE # Mean squared error for fit to data
DMSE <- Data2LDList$DpMSE # gradient with respect to parameter values
# These parameters control the optimization.
dbglev <- 1 # debugging level
iterlim <- 50 # maximum number of iterations
convrg <- c(1e-8, 1e-4) # convergence criterion
# This sets up an increasing sequence of rho values
gammavec <- 0:7
rhoMat <- as.matrix(exp(gammavec)/(1+exp(gammavec)),length(gammavec),1)
nrho <- length(rhoMat)
dfesave <- matrix(0,nrho,1)
gcvsave <- matrix(0,nrho,1)
MSEsave <- matrix(0,nrho,1)
thesave <- matrix(0,nrho,nparam)
# This initializes the list object containing coefficient estimates
TrayModelList.opt <- TrayModelList
# Loop through values of rho.
# for test purposes, only the first value rho = 0.5 is used here
for (irho in 1:1) {
rhoi <- rhoMat[irho,]
print(paste("rho <- ",round(rhoi,5)))
Data2LDResult <-
Data2LD.opt(TrayDataList, TrayBasisList, TrayModelList.opt,
rhoi, convrg, iterlim, dbglev)
theta.opti <- Data2LDResult$thetastore
TrayModelList.opt <- modelVec2List(TrayModelList.opt, theta.opti)
Data2LDList <- Data2LD(TrayDataList, TrayBasisList, TrayModelList.opt, rhoi)
MSE <- Data2LDList$MSE
df <- Data2LDList$df
gcv <- Data2LDList$gcv
ISE <- Data2LDList$ISE
Var.theta <- Data2LDList$Var.theta
thesave[irho,] <- theta.opti
dfesave[irho] <- df
gcvsave[irho] <- gcv
MSEsave[irho] <- MSE
}
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.