pcode | R Documentation |
Obtain estimates of both structural and nuisance parameters of an ODE model by parameter cascade method.
pcode(data, time, ode.model, par.names, state.names, likelihood.fun, par.initial, basis.list,lambda,controls)
data |
A data frame or a matrix contain observations from each dimension of the ODE model. |
time |
A vector contain observation times or a matrix if time points are different between dimensions. |
ode.model |
An R function that computes the time derivative of the ODE model given observations of states variable and structural parameters. |
par.names |
The names of structural parameters defined in the 'ode.model'. |
state.names |
The names of state variables defined in the 'ode.model'. |
likelihood.fun |
A likelihood function passed to PCODE in case of that the error terms do not have a Normal distribution. |
par.initial |
Initial value of structural parameters to be optimized. |
basis.list |
A list of basis objects for smoothing each dimension's observations. Can be the same or different across dimensions. |
lambda |
Penalty parameter for controling the fidelity of interpolation. |
controls |
A list of control parameters. See Details. |
The controls
argument is a list providing addition inputs for the nonlinear least square optimizer or general optimizer optim
:
nquadpts
Determine the number of quadrature points for approximating an integral. Default is 101.
smooth.lambda
Determine the smoothness penalty for obtaining initial value of nuisance parameters.
tau
Initial value of Marquardt parameter. Small values indicate good initial values for structural parameters.
tolx
Tolerance for parameters of objective functions. Default is set at 1e-6.
tolg
Tolerance for the gradient of parameters of objective functions. Default is set at 1e-6.
maxeval
The maximum number of evaluation of the outter optimizer. Default is set at 20.
structural.par |
The structural parameters of the ODE model. |
nuisance.par |
The nuisance parameters or the basis coefficients for interpolating observations. |
library(fda) library(deSolve) library(MASS) library(pracma) #Simple ode model example #define model parameters model.par <- c(theta = c(0.1)) #define state initial value state <- c(X = 0.1) #Define model for function 'ode' to numerically solve the system ode.model <- function(t, state,parameters){ with(as.list(c(state,parameters)), { dX <- theta*X*(1-X/10) return(list(dX)) }) } #Observation time points times <- seq(0,100,length.out=101) #Solve the ode model desolve.mod <- ode(y=state,times=times,func=ode.model,parms = model.par) #Prepare for doing parameter cascading method #Generate basis object for interpolation and as argument of pcode #21 konts equally spaced within [0,100] knots <- seq(0,100,length.out=21) #order of basis functions norder <- 4 #number of basis funtions nbasis <- length(knots) + norder - 2 #creating Bspline basis basis <- create.bspline.basis(c(0,100),nbasis,norder,breaks = knots) #Add random noise to ode solution for simulating data nobs <- length(times) scale <- 0.1 noise <- scale*rnorm(n = nobs, mean = 0 , sd = 1) observation <- desolve.mod[,2] + noise #parameter estimation pcode(data = observation, time = times, ode.model = ode.model, par.initial = 0.1, par.names = 'theta',state.names = 'X', basis.list = basis, lambda = 1e2)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.