| Title: | Continuous Time Stochastic Modelling using Template Model Builder |
|---|---|
| Description: | Perform state and parameter inference, and forecasting, in stochastic state-space systems using the 'ctsmTMB' R6 class. This class provides a user-friendly interface for working with stochastic state space models. Inference is based on maximum likelihood estimation, with derivatives efficiently computed through automatic differentiation enabled by the 'TMB'/'RTMB' packages (Kristensen et al., 2016) <doi:10.18637/jss.v070.i05>. The available inference methods include Kalman filters, in addition to a Laplace approximation-based smoothing method. For further details of these methods refer to the documentation of the 'CTSMR' package <https://ctsm.info/ctsmr-reference.pdf> and Thygesen (2025) <doi:10.48550/arXiv.2503.21358>. Forecasting capabilities include moment predictions and stochastic path simulations implemented in 'C++' using 'Rcpp' (Eddelbuettel et al., 2018) <doi:10.1080/00031305.2017.1375990> for computational efficiency. |
| Authors: | Phillip Vetter [aut, cre, cph], Jan Møller [ctb], Uffe Thygesen [ctb], Peder Bacher [ctb], Henrik Madsen [ctb] |
| Maintainer: | Phillip Vetter <[email protected]> |
| License: | GPL-3 |
| Version: | 1.1.0 |
| Built: | 2026-05-23 08:58:41 UTC |
| Source: | https://github.com/phillipbvetter/ctsmTMB |
The following public methods are used to construct a stochastic state space model system, consisting of a set of stochastic differential equations (SDEs), and one or more algebraic observation equations (AOEs). The AOEs are used to infer information about the value of the (latent) states governed by the SDEs, and thus must be functions of at least one state.
The function returns an object of class R6 and ctsmTMB,
which can be used to define a stochastic state space system. The model object has methods
likelihood and estimate for parameter estimation, filter for state reconstruction,
and predict / simulate for forecasting.
new()
Initialize private fields
ctsmTMB$new()
.private()
Extract the private fields of a ctsmTMB model object. Primarily used for debugging.
ctsmTMB$.private()
getPrivateFields()
Extract the private fields of a ctsmTMB model object. Primarily used for debugging.
ctsmTMB$getPrivateFields()
getPrivate()
Extract the private fields of a ctsmTMB model object. Primarily used for debugging.
ctsmTMB$getPrivate()
setTrainingMethod()
Set training method
ctsmTMB$setTrainingMethod(full.prediction)
full.predictionboolean whether or not to train against a full prediction over the data, rather than 1-step predictions
addSystem()
Define stochastic differential equation(s) on the form
d<state> ~ f(t,<states>, <inputs>) * dt + g(t, <states>, <inputs>) * dw
ctsmTMB$addSystem(form, ...)
forma formula specifying a stochastic differential equation
...additional formulas similar to form for specifying
multiple equations at once.
addObs()
Define algebraic observation equations on the form
<observation> ~ h(t, <states>, <inputs>) + e)
where h is the observation function, and e is normally
distributed noise with zero mean.
This function only specifies the observation name, and its mean through h.
ctsmTMB$addObs(form, ..., obsnames = NULL)
forma formula specifying an observation equation
...additional formulas similar to form for specifying
multiple equations at once.
obsnamescharacter vector specifying the name of the observation.
This is used when the left-hand side of form consists of more than just
a single variable (of class 'call').
setVariance()
Specify the variance of an observation equation.
A defined observation variable y in e.g. addObs(y ~
h(t,<states>,<inputs>) is perturbed by Gaussian noise with zero mean and
variance
to-be specified using setVariance(y ~ p(t,<states>,<inputs>).
We can for instance declare setVariance(y ~ sigma_x^2
where sigma_x is a fixed effect parameter to be declared through
setParameter.
ctsmTMB$setVariance(form, ...)
formformula class specifying the observation equation to be added to the system.
...additional formulas identical to form to specify multiple
observation equations at a time.
addInput()
Declare variables as data inputs
Declare whether a variable contained in system, observation or observation
variance equations is an input variable. If e.g. the system equation contains
an input variable u then it is declared using addInput(u).
The input u must be contained in the data.frame .data provided
when calling the estimate or predict methods.
ctsmTMB$addInput(...)
...variable names that specifies the name of input variables in the defined system.
setParameter()
Declare which variables that are (fixed effects) parameters in the specified model, and specify the initial optimizer guess, as well as lower / upper bounds during optimization. There are two ways to declare parameters:
You can declare parameters using formulas i.e. setParameter(
theta = c(1,0,10), mu = c(0,-10,10) ). The first value is the initial
value for the optimizer, the second value is the lower optimization
bound and the third value is the upper optimization bound.
You can provide a 3-column matrix where rows corresponds to different parameters, and the parameter names are provided as rownames of the matrix. The columns values corresponds to the description in the vector format above.
ctsmTMB$setParameter(...)
...a named vector or matrix as described above.
setAlgebraics()
Add algebraic relations.
Algebraic relations is a convenient way to transform parameters in your equations.
In the Ornstein-Uhlenbeck process the rate parameter theta is always positive, so
estimation in the log-domain is a good idea. Instead of writing exp(theta) directly
in the system equation one can transform into the log domain using the algebraic relation
setAlgebraics(theta ~ exp(logtheta)). All instances of theta is replaced
by exp(logtheta) when compiling the C++ function. Note that you must provide values
for logtheta now instead of theta when declaring parameters through
setParameter
ctsmTMB$setAlgebraics(form, ...)
formalgebraic formula
...additional formulas
setInitialState()
Declare the initial state values i.e. mean and covariance for the system states.
ctsmTMB$setInitialState(initial.state)
initial.statea named list of two entries 'x0' and 'p0' containing the initial state and covariance of the state.
setInitialVarianceScaling()
A scalar value that is multiplied onto the estimated initial state covariance matrix. The scaling is only applied when the initial state/cov is estimated, not when it is set by the user.
ctsmTMB$setInitialVarianceScaling(scaling)
scalinga numeric scalar value.
setLamperti()
Set a Lamperti Transformation
If the provided system equations have state dependent diffusion in a few available ways then it is advantageous to perform a transformation to remove the state dependence. This comes at the cost of a more complicated drift function. The following types of state-dependence is currently supported
'identity' - when the diffusion is state-independent (default)
'log' - when the diffusion is proportional to to x * dw
'logit' - when the diffusion is proportional to x * (1-x) * dw
'sqrt-logit' - when the diffusion is proportional to sqrt(x * (1-x)) * dw
ctsmTMB$setLamperti(transforms, states = NULL)
transformscharacter vector - one of either "identity, "log", "logit", "sqrt-logit"
statesa vector of the state names for which the specified transformations should be applied to.
setModelname()
Set modelname used to create the C++ file for TMB
When calling TMB::MakeADFun the (negative log) likelihood function
is created in the directory specified by the setCppfilesDirectory
method with name <modelname>.cpp
ctsmTMB$setModelname(name)
namestring defining the model name.
setMAP()
Enable maximum a posterior (MAP) estimation.
Adds a maximum a posterior contribution to the (negative log) likelihood
function by evaluating the fixed effects parameters in a multivariate Gaussian
with mean and covariance as provided.
ctsmTMB$setMAP(mean, cov)
meanmean vector of the Gaussian prior parameter distribution
covcovariance matrix of the Gaussian prior parameter distribution
setAdvancedSettings()
Enable maximum a posterior (MAP) estimation.
Adds a maximum a posterior contribution to the (negative log) likelihood
function by evaluating the fixed effects parameters in a multivariate Gaussian
with mean and covariance as provided.
ctsmTMB$setAdvancedSettings( force.ad = TRUE, rtmb.tapeconfig = NULL, tmb.tapeconfig = NULL )
force.ada boolean indicating whether to use state space functions that take advantage of the RTMB::AD(...,force=TRUE) hack which reduces compilation time call to MakeADFun by 20%. This breaks some functionalities such as REPORT.
rtmb.tapeconfigoptions to be passed to TapeConfig.
tmb.tapeconfigoptions to be passed to config.
getSystems()
Retrieve system equations.
ctsmTMB$getSystems()
getObservations()
Retrieve observation equations.
ctsmTMB$getObservations()
getVariances()
Retrieve observation variances
ctsmTMB$getVariances()
getAlgebraics()
Retrieve algebraic relations
ctsmTMB$getAlgebraics()
getInitialState()
Retrieve initially set state and covariance
ctsmTMB$getInitialState()
getParameters()
Get initial (and estimated) parameters.
ctsmTMB$getParameters(type = "all", value = "all")
typeone of "all", "free" or "fixed" parameters.
valueone of "all", "initial", "estimate", "lower" or "upper". When not "all", a named numeric vector is returned.
getTimers()
Retrieve initially timers
ctsmTMB$getTimers()
getEstimate()
Retrieve initially set state and covariance
ctsmTMB$getEstimate()
getLikelihood()
Retrieve initially set state and covariance
ctsmTMB$getLikelihood()
getPrediction()
Retrieve initially set state and covariance
ctsmTMB$getPrediction()
getSimulation()
Retrieve initially set state and covariance
ctsmTMB$getSimulation()
filter()
Perform state filtering (or smoothing for the 'laplace' method)
ctsmTMB$filter( data, pars = NULL, method = "ekf", ode.solver = "rk4", ode.timestep = diff(data$t), loss = "quadratic", loss_c = NULL, ukf.hyperpars = c(1, 0, 3), initial.state = self$getInitialState(), laplace.residuals = FALSE, estimate.initial.state = FALSE, use.cpp = TRUE, silent = FALSE, ... )
datadata.frame containing time-vector 't', observations and inputs. The observations
can take NA-values.
parsfixed parameter vector parsed to the objective function for prediction/filtration. The default
parameter values used are the initial parameters provided through setParameter, unless the estimate
methodcharacter vector specifying the filtering method used for state/likelihood calculations. Must be one of either "lkf", "ekf", "laplace".
ode.solverSets the ODE solver used in the Kalman Filter methods for solving the moment differential equations. The default "euler" is the Forward Euler method, alternatively the classical 4th order Runge Kutta method is available via "rk4".
ode.timestepnumeric value. Sets the time step-size in numerical filtering schemes.
The defined step-size is used to calculate the number of steps between observation time-points as
defined by the provided data. If the calculated number of steps is larger than N.01 where N
is an integer, then the time-step is reduced such that exactly N+1 steps is taken between observations
The step-size is used in the two following ways depending on the
chosen method:
Kalman filters: The time-step is used as the step-size in the numerical Forward-Euler scheme to compute the prior state mean and covariance estimate as the final time solution to the first and second order moment differential equations.
TMB method: The time-step is used as the step-size in the Euler-Maruyama scheme for simulating a sample path of the stochastic differential equation, which serves to link together the latent (random effects) states.
losscharacter vector. Sets the loss function type (only implemented for the kalman filter methods). The loss function is per default quadratic in the one-step residuals as is natural when the Gaussian (negative log) likelihood is evaluated, but if the tails of the distribution is considered too small i.e. outliers are weighted too much, then one can choose loss functions that accounts for this. The three available types available:
Quadratic loss (quadratic).
Quadratic-Linear (huber)
Quadratic-Constant (tukey)
The cutoff for the Huber and Tukey loss functions are determined from a provided cutoff
parameter loss_c. The implementations of these losses are approximations (pseudo-huber and sigmoid
approximation respectively) for smooth derivatives.
loss_ccutoff value for huber and tukey loss functions. Defaults to c=3
ukf.hyperparsThe hyperparameters alpha, beta, and kappa used for sigma points and weights construction in the Unscented Kalman Filter.
initial.statea named list of two entries 'x0' and 'p0' containing the initial state and covariance of the state
laplace.residualsboolean - whether or not to calculate one-step ahead residuals using the method of oneStepPredict.
estimate.initial.stateboolean value. When TRUE the initial state and covariance matrices are estimated as the stationary solution of the linearized mean and covariance differential equations. When the system contains time-varying inputs, the first element of these is used.
use.cppa boolean to indicate whether to use C++ to perform calculations
silentlogical value whether or not to suppress printed messages such as 'Checking Data', 'Building Model', etc. Default behaviour (FALSE) is to print the messages.
...additional arguments
smoother()
Perform state filtering (or smoothing for the 'laplace' method)
ctsmTMB$smoother( data, pars = NULL, method = "laplace", ode.solver = "euler", ode.timestep = diff(data$t), loss = "quadratic", loss_c = NULL, initial.state = self$getInitialState(), laplace.residuals = FALSE, estimate.initial.state = FALSE, silent = FALSE, ... )
datadata.frame containing time-vector 't', observations and inputs. The observations
can take NA-values.
parsfixed parameter vector parsed to the objective function for prediction/filtration. The default
parameter values used are the initial parameters provided through setParameter, unless the estimate
methodcharacter vector specifying the filtering method used for state/likelihood calculations. Must be one of either "lkf", "ekf", "laplace".
ode.solverSets the ODE solver used in the Kalman Filter methods for solving the moment differential equations. The default "euler" is the Forward Euler method, alternatively the classical 4th order Runge Kutta method is available via "rk4".
ode.timestepnumeric value. Sets the time step-size in numerical filtering schemes.
The defined step-size is used to calculate the number of steps between observation time-points as
defined by the provided data. If the calculated number of steps is larger than N.01 where N
is an integer, then the time-step is reduced such that exactly N+1 steps is taken between observations
The step-size is used in the two following ways depending on the
chosen method:
Kalman filters: The time-step is used as the step-size in the numerical Forward-Euler scheme to compute the prior state mean and covariance estimate as the final time solution to the first and second order moment differential equations.
TMB method: The time-step is used as the step-size in the Euler-Maruyama scheme for simulating a sample path of the stochastic differential equation, which serves to link together the latent (random effects) states.
losscharacter vector. Sets the loss function type (only implemented for the kalman filter methods). The loss function is per default quadratic in the one-step residuals as is natural when the Gaussian (negative log) likelihood is evaluated, but if the tails of the distribution is considered too small i.e. outliers are weighted too much, then one can choose loss functions that accounts for this. The three available types available:
Quadratic loss (quadratic).
Quadratic-Linear (huber)
Quadratic-Constant (tukey)
The cutoff for the Huber and Tukey loss functions are determined from a provided cutoff
parameter loss_c. The implementations of these losses are approximations (pseudo-huber and sigmoid
approximation respectively) for smooth derivatives.
loss_ccutoff value for huber and tukey loss functions. Defaults to c=3
initial.statea named list of two entries 'x0' and 'p0' containing the initial state and covariance of the state
laplace.residualsboolean - whether or not to calculate one-step ahead residuals using the method of oneStepPredict.
estimate.initial.stateboolean value. When TRUE the initial state and covariance matrices are estimated as the stationary solution of the linearized mean and covariance differential equations. When the system contains time-varying inputs, the first element of these is used.
silentlogical value whether or not to suppress printed messages such as 'Checking Data', 'Building Model', etc. Default behaviour (FALSE) is to print the messages.
...additional arguments
estimate()
Estimate the fixed effects parameters in the specified model.
ctsmTMB$estimate( data, method = "ekf", ode.solver = "rk4", ode.timestep = diff(data$t), loss = "quadratic", loss_c = NULL, ukf.hyperpars = c(1, 0, 3), initial.state = self$getInitialState(), trace = 10, control = list(trace = trace, iter.max = 1e+05, eval.max = 1e+05), use.hessian = FALSE, report = TRUE, laplace.residuals = FALSE, unconstrained.optim = FALSE, estimate.initial.state = FALSE, silent = FALSE, compile = FALSE, ... )
dataa data.frame with a time-column (must be named "t"), observations and
inputs used maximum-likelihood parameter and state estimation. The observations can take NA-values.
methoda character string specifying the likelihood method used for parameter and state estimation.
ode.solvera character string to determine the ODE solver scheme for the Kalman methods when solving the
moment differential equations. The methods are either Forward Euler ("euler") or 4th order Runge-Kutta ("rk4")
(default).
ode.timestepa numeric value that determines the ODE solver time step-size. The passed value
is used to calculate the number of steps between time-points in data$t such that an integer number of
steps are taken. For details see the "Estimation" vignette on the webpage. The step-size is used for solving
the moment differential equations (Kalman methods) or for simulating the stochastic SDE path (Laplace methods).
losscharacter vector. Sets the loss function type (only implemented for the kalman filter methods). The loss function is per default quadratic in the one-step residuals as is natural when the Gaussian (negative log) likelihood is evaluated, but if the tails of the distribution is considered too small i.e. outliers are weighted too much, then one can choose loss functions that accounts for this. The three available types available:
Quadratic loss (quadratic).
Quadratic-Linear (huber)
Quadratic-Constant (tukey)
The cutoff for the Huber and Tukey loss functions are determined from a provided cutoff
parameter loss_c. The implementations of these losses are approximations (pseudo-huber and sigmoid
approximation respectively) for smooth derivatives.
loss_ccutoff value for huber and tukey loss functions. Defaults to c=3
ukf.hyperparsThe hyperparameters alpha, beta, and kappa used for sigma points and weights construction in the Unscented Kalman Filter.
initial.statea named list of two entries 'x0' and 'p0' containing the initial state and covariance of the state
traceinteger passed to control which determines number of steps between each print-out
during optimization (use 0 to disable tracing print-outs).
controllist of control parameters parsed to nlminb as its control argument.
See ?stats::nlminb for more information
use.hessianboolean value. The default (TRUE) causes the optimization algorithm
stats::nlminb to use the fixed effects hessian of the (negative log) likelihood when
performing the optimization. This feature is only available for the kalman filter methods
without any random effects.
reportboolean - whether or not to report filtered states, observations and residuals.
laplace.residualsboolean - whether or not to calculate one-step ahead residuals using the method of oneStepPredict.
unconstrained.optimboolean value. When TRUE then the optimization is carried out unconstrained i.e.
without any of the parameter bounds specified during setParameter.
estimate.initial.stateboolean value. When TRUE the initial state and covariance matrices are estimated as the stationary solution of the linearized mean and covariance differential equations. When the system contains time-varying inputs, the first element of these is used.
silentlogical value whether or not to suppress printed messages such as 'Checking Data', 'Building Model', etc. Default behaviour (FALSE) is to print the messages.
compileboolean for (re)compiling the objective C++ file, used for methods ending with _cpp.
...additional arguments
likelihood()
Construct and extract function handlers for the negative log likelihood function.
The handlers from TMB's MakeADFun are constructed and returned.
This enables the user to e.g. choose their own optimization algorithm, or just
have more control of the optimization workflow.
ctsmTMB$likelihood( data, method = "ekf", ode.solver = "rk4", ode.timestep = diff(data$t), loss = "quadratic", loss_c = NULL, ukf.hyperpars = c(1, 0, 3), initial.state = self$getInitialState(), estimate.initial.state = FALSE, silent = FALSE, compile = FALSE, ... )
dataa data.frame containing time-vector 't', observations and inputs.
The observations can take NA-values.
methodcharacter vector specifying the filtering method used for state/likelihood calculations. Must be one of either "lkf", "ekf", "laplace".
ode.solverSets the ODE solver used in the Kalman Filter methods for solving the moment differential equations. The default "euler" is the Forward Euler method, alternatively the classical 4th order Runge Kutta method is available via "rk4".
ode.timestepthe time-step used in the filtering schemes. The time-step has two different uses depending on the chosen method.
Kalman Filters: The time-step is used when numerically solving the moment differential equations.
Laplace Approximation: The time-step is used in the Euler-Maruyama simulation scheme for simulating a sample path of the stochastic differential equation, which serves to link together the latent (random effects) states.
The defined step-size is used to calculate the number of steps between
observation time-points as
defined by the provided data. If the calculated number of steps is larger than N.01 where N
is an integer, then the time-step is reduced such that exactly N+1 steps is taken between observations
The step-size is used in the two following ways depending on the
chosen method:
Kalman filters: The time-step is used as the step-size in the numerical Forward-Euler scheme to compute the prior state mean and covariance estimate as the final time solution to the first and second order moment differential equations.
TMB method: The time-step is used as the step-size in the Euler-Maruyama scheme for simulating a sample path of the stochastic differential equation, which serves to link together the latent (random effects) states.
losscharacter vector. Sets the loss function type (only implemented for the kalman filter methods). The loss function is per default quadratic in the one-step residuals as is natural when the Gaussian (negative log) likelihood is evaluated, but if the tails of the distribution is considered too small i.e. outliers are weighted too much, then one can choose loss functions that accounts for this. The three available types available:
Quadratic loss (quadratic).
Quadratic-Linear (huber)
Quadratic-Constant (tukey)
The cutoff for the Huber and Tukey loss functions are determined from a provided cutoff
parameter loss_c. The implementations of these losses are approximations (pseudo-huber and sigmoid
approximation respectively) for smooth derivatives.
loss_ccutoff value for huber and tukey loss functions. Defaults to c=3
ukf.hyperparsThe hyperparameters alpha, beta, and kappa used for sigma points and weights construction in the Unscented Kalman Filter.
initial.statea named list of two entries 'x0' and 'p0' containing the initial state and covariance of the state
estimate.initial.stateboolean value. When TRUE the initial state and covariance matrices are estimated as the stationary solution of the linearized mean and covariance differential equations. When the system contains time-varying inputs, the first element of these is used.
silentlogical value whether or not to suppress printed messages such as 'Checking Data', 'Building Model', etc. Default behaviour (FALSE) is to print the messages.
compileboolean for (re)compiling the objective C++ file, used for methods ending with _cpp.
...additional arguments
The method retunrns a list which includes function handles for the likelihood function.
A vector of initial parameter values. This is the expected input size to fn, gr and he.
Negative log-likelihood function.
Negative log-likelihood gradient.
Negative log-likelihood hessian.
predict()
Perform prediction/filtration to obtain state mean and covariance estimates. The predictions are
obtained by solving the moment equations k.ahead steps forward in time when using the current step posterior
state estimate as the initial condition.
ctsmTMB$predict(
data,
pars = NULL,
method = "ekf",
ode.solver = "rk4",
ode.timestep = diff(data$t),
k.ahead = nrow(data) - 1,
return.k.ahead = 0:min(k.ahead, nrow(data) - 1),
return.variance = c("marginal", "none", "covariance", "correlation"),
ukf.hyperpars = c(1, 0, 3),
initial.state = self$getInitialState(),
estimate.initial.state = private$estimate.initial,
use.cpp = TRUE,
silent = FALSE,
...
)datadata.frame containing time-vector 't', observations and inputs. The observations
can take NA-values.
parsfixed parameter vector parsed to the objective function for prediction/filtration. The default
parameter values used are the initial parameters provided through setParameter, unless the estimate
function has been run, then the default values will be those at the found optimum.
methodThe prediction method
ode.solverSets the ODE solver used in the Kalman Filter methods for solving the moment differential equations. The default "euler" is the Forward Euler method, alternatively the classical 4th order Runge Kutta method is available via "rk4".
ode.timestepnumeric value. Sets the time step-size in numerical filtering schemes.
The defined step-size is used to calculate the number of steps between observation time-points as
defined by the provided data. If the calculated number of steps is larger than N.01 where N
is an integer, then the time-step is reduced such that exactly N+1 steps is taken between observations
The step-size is used in the two following ways depending on the
chosen method:
Kalman filters: The time-step is used as the step-size in the numerical Forward-Euler scheme to compute the prior state mean and covariance estimate as the final time solution to the first and second order moment differential equations.
TMB method: The time-step is used as the step-size in the Euler-Maruyama scheme for simulating a sample path of the stochastic differential equation, which serves to link together the latent (random effects) states.
k.aheadinteger specifying the desired number of time-steps (as determined by the provided data time-vector) for which predictions are made (integrating the moment ODEs forward in time without data updates).
return.k.aheadnumeric vector of integers specifying which k.ahead predictions to that should be returned.
return.variancea string to indicate what kind of dispersions should be reported.
ukf.hyperparsThe hyperparameters alpha, beta, and kappa used for sigma points and weights construction in the Unscented Kalman Filter.
initial.statea named list of two entries 'x0' and 'p0' containing the initial state and covariance of the state
estimate.initial.statebool - stationary estimation of initial mean and covariance
use.cppa boolean to indicate whether to use C++ to perform calculations
silentlogical value whether or not to suppress printed messages such as 'Checking Data', 'Building Model', etc. Default behaviour (FALSE) is to print the messages.
...additional arguments
A data.frame that contains for each time step the posterior state estimate at that time.step (k = 0), and the
prior state predictions (k = 1,...,k.ahead). If return.dispersion = TRUE then the state covariance/correlation
matrix is returned, otherwise only the marginal variances are returned.
simulate()
Perform prediction/filtration to obtain state mean and covariance estimates. The predictions are
obtained by solving the moment equations k.ahead steps forward in time when using the current step posterior
state estimate as the initial condition.
ctsmTMB$simulate(
data,
pars = NULL,
use.cpp = TRUE,
cpp.seeds = NULL,
method = c("ekf", "lkf", "ukf", "laplace", "laplace.thygesen"),
ode.solver = "rk4",
ode.timestep = diff(data$t),
simulation.timestep = diff(data$t),
k.ahead = nrow(data) - 1,
return.k.ahead = 0:min(k.ahead, nrow(data) - 1),
n.sims = 100,
ukf.hyperpars = c(1, 0, 3),
initial.state = self$getInitialState(),
estimate.initial.state = private$estimate.initial,
silent = FALSE,
...
)datadata.frame containing time-vector 't', observations and inputs. The observations
can take NA-values.
parsfixed parameter vector parsed to the objective function for prediction/filtration. The default
parameter values used are the initial parameters provided through setParameter, unless the estimate
function has been run, then the default values will be those at the found optimum.
use.cppa boolean to indicate whether to use C++ to perform calculations
cpp.seedsan integer seed value to control RNG normal draws on the C++ side.
methodThe natural TMB-style formulation where latent states are considered random effects and are integrated out using the Laplace approximation. This method only yields the gradient of the (negative log) likelihood function with respect to the fixed effects for optimization. The method is slower although probably has some precision advantages, and allows for non-Gaussian observation noise (not yet implemented). One-step / K-step residuals are not yet available in the package.
(Continuous-Discrete) Extended Kalman Filter where the system dynamics are linearized to handle potential non-linearities. This is computationally the fastest method.
(Continuous-Discrete) Unscented Kalman Filter. This is a higher order non-linear Kalman Filter which improves the mean and covariance estimates when the system display high nonlinearity, and circumvents the necessity to compute the Jacobian of the drift and observation functions.
All package features are currently available for the kalman filters, while TMB is limited to
parameter estimation. In particular, it is straight-forward to obtain k-step-ahead predictions
with these methods (use the predict S3 method), and stochastic simulation is also available
in the cases where long prediction horizons are sought, where the normality assumption will be
inaccurate
ode.solverSets the ODE solver used in the Kalman Filter methods for solving the moment differential equations. The default "euler" is the Forward Euler method, alternatively the classical 4th order Runge Kutta method is available via "rk4".
ode.timestepnumeric value. Sets the time step-size in numerical filtering schemes.
The defined step-size is used to calculate the number of steps between observation time-points as
defined by the provided data. If the calculated number of steps is larger than N.01 where N
is an integer, then the time-step is reduced such that exactly N+1 steps is taken between observations
The step-size is used in the two following ways depending on the
chosen method:
Kalman filters: The time-step is used as the step-size in the numerical Forward-Euler scheme to compute the prior state mean and covariance estimate as the final time solution to the first and second order moment differential equations.
TMB method: The time-step is used as the step-size in the Euler-Maruyama scheme for simulating a sample path of the stochastic differential equation, which serves to link together the latent (random effects) states.
simulation.timesteptimestep used in the euler-maruyama scheme
k.aheadinteger specifying the desired number of time-steps (as determined by the provided data time-vector) for which predictions are made (integrating the moment ODEs forward in time without data updates).
return.k.aheadnumeric vector of integers specifying which k.ahead predictions to that should be returned.
n.simsnumber of simulations
ukf.hyperparsThe hyperparameters alpha, beta, and kappa used for sigma points and weights construction in the Unscented Kalman Filter.
initial.statea named list of two entries 'x0' and 'p0' containing the initial state and covariance of the state
estimate.initial.statebool - stationary estimation of initial mean and covariance
silentlogical value whether or not to suppress printed messages such as 'Checking Data', 'Building Model', etc. Default behaviour (FALSE) is to print the messages.
...additional arguments
return.dispersionboolean value to indicate whether the covariance (instead of the correlation) should be returned.
A data.frame that contains for each time step the posterior state estimate at that time.step (k = 0), and the
prior state predictions (k = 1,...,k.ahead). If return.dispersion = TRUE then the state covariance/correlation
matrix is returned, otherwise only the marginal variances are returned.
print()
Function to print the model object
ctsmTMB$print()
clone()
The objects of this class are cloneable with this method.
ctsmTMB$clone(deep = FALSE)
deepWhether to make a deep clone.
library(ctsmTMB) model <- ctsmTMB$new() # adding a single system equations model$addSystem(dx ~ theta * (mu+u-x) * dt + sigma_x*dw) # adding an observation equation and setting variance model$addObs(y ~ x) model$setVariance(y ~ sigma_y^2) # add model input model$addInput(u) # add parameters model$setParameter( theta = c(initial = 1, lower=1e-5, upper=50), mu = c(initial=1.5, lower=0, upper=5), sigma_x = c(initial=1, lower=1e-10, upper=30), sigma_y = 1e-2 ) # set the model initial state model$setInitialState(list(1,1e-1)) # extract the likelihood handlers nll <- model$likelihood(data=Ornstein) # calculate likelihood, gradient and hessian w.r.t parameters nll$fn(nll$par) nll$gr(nll$par) nll$he(nll$par) # estimate the parameters using an extended kalman filter fit <- model$estimate(data=Ornstein) # perform moment predictions pred <- model$predict(data=Ornstein) # perform stochatic simulations sim <- model$simulate(data=Ornstein, n.sims=10)library(ctsmTMB) model <- ctsmTMB$new() # adding a single system equations model$addSystem(dx ~ theta * (mu+u-x) * dt + sigma_x*dw) # adding an observation equation and setting variance model$addObs(y ~ x) model$setVariance(y ~ sigma_y^2) # add model input model$addInput(u) # add parameters model$setParameter( theta = c(initial = 1, lower=1e-5, upper=50), mu = c(initial=1.5, lower=0, upper=5), sigma_x = c(initial=1, lower=1e-10, upper=30), sigma_y = 1e-2 ) # set the model initial state model$setInitialState(list(1,1e-1)) # extract the likelihood handlers nll <- model$likelihood(data=Ornstein) # calculate likelihood, gradient and hessian w.r.t parameters nll$fn(nll$par) nll$gr(nll$par) nll$he(nll$par) # estimate the parameters using an extended kalman filter fit <- model$estimate(data=Ornstein) # perform moment predictions pred <- model$predict(data=Ornstein) # perform stochatic simulations sim <- model$simulate(data=Ornstein, n.sims=10)
Create a ctsmTMB model faster avoiding $...
newModel()newModel()
Print of ctsmTMB model object
The data was simulated using a standard Euler-Maruyama method.
The simulated process is governed by the SDE #' dx ~ theta * (mu + u - x) * dt + sigma_x * dw
The parameters used for simulation were theta = 5, mu = 3, sigma_x = 1, sigma_y = 0.1
The simulation time-step was 1e-3, and observation time-step 1e-1. The simulation was taken from t = 0..20
The simulated input was given by u.sim = cumsum(rnorm(length(t.sim),sd=0.05))
where t.sim is the simulated time vector.
OrnsteinOrnstein
A data frame of 201 rows and 3 columns. The columns represent the
variables: t (time), y (observation) and u (input).
The data was simulated using a standard Euler-Maruyama method.
The simulated process is governed by the coupled SDEs
dx1 ~ theta * (mu + u - x1) * dt + sigma_x1 * dw1 dx2 ~ alpha * (x1 - x2) * dt + sigma_x2 * dw2
with observations
y1 ~ x1 y2 ~ x2
x2 mean-reverts toward x1, acting as a lagged/smoothed version of x1. The dataset is intended for testing multi-dimensional (> 1 state) inference.
The parameters used for simulation were theta = 5, mu = 3, alpha = 2, sigma_x1 = 1, sigma_x2 = 0.5, sigma_y1 = 0.1, sigma_y2 = 0.1
The simulation time-step was 1e-3, and observation time-step 1e-1. The simulation was taken from t = 0..20
The simulated input was given by u.sim = cumsum(rnorm(length(t.sim), sd=0.05))
where t.sim is the simulated time vector.
Ornstein_augmentedOrnstein_augmented
A data frame of 201 rows and 4 columns. The columns represent the
variables: t (time), y1 (observation of x1), y2 (observation of x2),
and u (input).
This function creates residual plots for an estimated ctsmTMB object
## S3 method for class 'ctsmTMB.fit' plot( x, print.plot = 1, type = "residuals", state.type = "prior", against.obs = NULL, ggtheme = getggplot2theme(), ylims = c(NA, NA), residual.burnin = 0L, residual.vs.obs.and.inputs = FALSE, ... )## S3 method for class 'ctsmTMB.fit' plot( x, print.plot = 1, type = "residuals", state.type = "prior", against.obs = NULL, ggtheme = getggplot2theme(), ylims = c(NA, NA), residual.burnin = 0L, residual.vs.obs.and.inputs = FALSE, ... )
x |
A R6 ctsmTMB fit object |
print.plot |
a single integer determining which element out of all
states/observations (depending on the argument to |
type |
a character vector either 'residuals' or 'states' determining what to plot. |
state.type |
a character vector either 'prior', 'posterior' or 'smoothed' determining what kind of states to plot. |
against.obs |
name of an observation to plot state predictions against. |
ggtheme |
ggplot2 theme to use for creating the ggplot. |
ylims |
limits on the y-axis for residual time-series plot |
residual.burnin |
integer N to remove the first N residuals |
residual.vs.obs.and.inputs |
the residual plots also include a new window with time-series plots of residuals, associated observations and inputs |
... |
additional arguments |
a (list of) ggplot residual plot(s)
library(ctsmTMB) model <- ctsmTMB$new() # create model model$addSystem(dx ~ theta * (mu+u-x) * dt + sigma_x*dw) model$addObs(y ~ x) model$setVariance(y ~ sigma_y^2) model$addInput(u) model$setParameter( theta = c(initial = 1, lower=1e-5, upper=50), mu = c(initial=1.5, lower=0, upper=5), sigma_x = c(initial=1, lower=1e-10, upper=30), sigma_y = 1e-2 ) model$setInitialState(list(1,1e-1)) # fit model to data fit <- model$estimate(Ornstein) # plot residuals ## Not run: plot(fit) plot(fit) # plot filtered states ## Not run: plot(fit, type="states")library(ctsmTMB) model <- ctsmTMB$new() # create model model$addSystem(dx ~ theta * (mu+u-x) * dt + sigma_x*dw) model$addObs(y ~ x) model$setVariance(y ~ sigma_y^2) model$addInput(u) model$setParameter( theta = c(initial = 1, lower=1e-5, upper=50), mu = c(initial=1.5, lower=0, upper=5), sigma_x = c(initial=1, lower=1e-10, upper=30), sigma_y = 1e-2 ) model$setInitialState(list(1,1e-1)) # fit model to data fit <- model$estimate(Ornstein) # plot residuals ## Not run: plot(fit) plot(fit) # plot filtered states ## Not run: plot(fit, type="states")
Plot of k-step predictions from a ctsmTMB prediction object
## S3 method for class 'ctsmTMB.pred' plot( x, y, k.ahead = unique(x$states[, "k.ahead"]), state.name = NULL, type = "states", against = NULL, ... )## S3 method for class 'ctsmTMB.pred' plot( x, y, k.ahead = unique(x$states[, "k.ahead"]), state.name = NULL, type = "states", against = NULL, ... )
x |
a ctsmTMB.pred object |
y |
not used |
k.ahead |
an integer indicating which k-ahead predictions to plot |
state.name |
a string indicating which states to plot |
type |
one of 'states' or 'observations', to plot |
against |
name of an observations to plot predictions against |
... |
additional arguments |
A plot of predicted states
library(ctsmTMB) model <- ctsmTMB$new() # create model model$addSystem(dx ~ theta * (mu+u-x) * dt + sigma_x*dw) model$addObs(y ~ x) model$setVariance(y ~ sigma_y^2) model$addInput(u) model$setParameter( theta = c(initial = 1, lower=1e-5, upper=50), mu = c(initial=1.5, lower=0, upper=5), sigma_x = c(initial=1, lower=1e-10, upper=30), sigma_y = 1e-2 ) model$setInitialState(list(1,1e-1)) # fit model to data fit <- model$estimate(Ornstein) # perform moment predictions pred <- model$predict(Ornstein) # plot the k.ahead=10 predictions plot(pred, against="y.data") # plot filtered states plot(fit, type="states", against="y")library(ctsmTMB) model <- ctsmTMB$new() # create model model$addSystem(dx ~ theta * (mu+u-x) * dt + sigma_x*dw) model$addObs(y ~ x) model$setVariance(y ~ sigma_y^2) model$addInput(u) model$setParameter( theta = c(initial = 1, lower=1e-5, upper=50), mu = c(initial=1.5, lower=0, upper=5), sigma_x = c(initial=1, lower=1e-10, upper=30), sigma_y = 1e-2 ) model$setInitialState(list(1,1e-1)) # fit model to data fit <- model$estimate(Ornstein) # perform moment predictions pred <- model$predict(Ornstein) # plot the k.ahead=10 predictions plot(pred, against="y.data") # plot filtered states plot(fit, type="states", against="y")
Plot a profile likelihood ctsmTMB object
## S3 method for class 'ctsmTMB.profile' plot(x, y, include.opt = TRUE, ...)## S3 method for class 'ctsmTMB.profile' plot(x, y, include.opt = TRUE, ...)
x |
a profile.ctsmTMB object |
y |
not in use |
include.opt |
boolean which indicates whether or not to include the total likelihood optimizer in the plot. |
... |
additional arguments |
library(ctsmTMB) model <- ctsmTMB$new() # create model model$addSystem(dx ~ theta * (mu+u-x) * dt + sigma_x*dw) model$addObs(y ~ x) model$setVariance(y ~ sigma_y^2) model$addInput(u) model$setParameter( theta = c(initial = 1, lower=1e-5, upper=50), mu = c(initial=1.5, lower=0, upper=5), sigma_x = c(initial=1, lower=1e-10, upper=30), sigma_y = 1e-2 ) model$setInitialState(list(1,1e-1)) # fit model to data fit <- model$estimate(Ornstein) # calculate profile likelihood # out <- profile(fit,parlist=list(theta=NULL)) out <- profile(fit,parlist=list(theta=NULL, mu=NULL)) # plot profile plot(out)library(ctsmTMB) model <- ctsmTMB$new() # create model model$addSystem(dx ~ theta * (mu+u-x) * dt + sigma_x*dw) model$addObs(y ~ x) model$setVariance(y ~ sigma_y^2) model$addInput(u) model$setParameter( theta = c(initial = 1, lower=1e-5, upper=50), mu = c(initial=1.5, lower=0, upper=5), sigma_x = c(initial=1, lower=1e-10, upper=30), sigma_y = 1e-2 ) model$setInitialState(list(1,1e-1)) # fit model to data fit <- model$estimate(Ornstein) # calculate profile likelihood # out <- profile(fit,parlist=list(theta=NULL)) out <- profile(fit,parlist=list(theta=NULL, mu=NULL)) # plot profile plot(out)
Basic print of ctsmTMB objects
## S3 method for class 'ctsmTMB' print(x, ...)## S3 method for class 'ctsmTMB' print(x, ...)
x |
an object of class 'ctsmTMB' |
... |
additional arguments (not in use) |
Print of ctsmTMB model object
library(ctsmTMB) model <- ctsmTMB$new() # print empty model print(model) # add elements to model and see new print model$addSystem(dx ~ theta * (mu+u-x) * dt + sigma_x*dw) model$addObs(y ~ x) model$setVariance(y ~ sigma_y^2) model$addInput(u) model$setParameter( theta = c(initial = 1, lower=1e-5, upper=50), mu = c(initial=1.5, lower=0, upper=5), sigma_x = c(initial=1, lower=1e-10, upper=30), sigma_y = 1e-2 ) print(model)library(ctsmTMB) model <- ctsmTMB$new() # print empty model print(model) # add elements to model and see new print model$addSystem(dx ~ theta * (mu+u-x) * dt + sigma_x*dw) model$addObs(y ~ x) model$setVariance(y ~ sigma_y^2) model$addInput(u) model$setParameter( theta = c(initial = 1, lower=1e-5, upper=50), mu = c(initial=1.5, lower=0, upper=5), sigma_x = c(initial=1, lower=1e-10, upper=30), sigma_y = 1e-2 ) print(model)
Basic print of objects ctsmTMB fit objects
## S3 method for class 'ctsmTMB.fit' print(x, ...)## S3 method for class 'ctsmTMB.fit' print(x, ...)
x |
a ctsmTMB fit object |
... |
additional arguments |
Print of ctsmTMB fit object
library(ctsmTMB) model <- ctsmTMB$new() # create model model$addSystem(dx ~ theta * (mu+u-x) * dt + sigma_x*dw) model$addObs(y ~ x) model$setVariance(y ~ sigma_y^2) model$addInput(u) model$setParameter( theta = c(initial = 1, lower=1e-5, upper=50), mu = c(initial=1.5, lower=0, upper=5), sigma_x = c(initial=1, lower=1e-10, upper=30), sigma_y = 1e-2 ) model$setInitialState(list(1,1e-1)) # fit model to data fit <- model$estimate(Ornstein) # print fit print(fit)library(ctsmTMB) model <- ctsmTMB$new() # create model model$addSystem(dx ~ theta * (mu+u-x) * dt + sigma_x*dw) model$addObs(y ~ x) model$setVariance(y ~ sigma_y^2) model$addInput(u) model$setParameter( theta = c(initial = 1, lower=1e-5, upper=50), mu = c(initial=1.5, lower=0, upper=5), sigma_x = c(initial=1, lower=1e-10, upper=30), sigma_y = 1e-2 ) model$setInitialState(list(1,1e-1)) # fit model to data fit <- model$estimate(Ornstein) # print fit print(fit)
Performs full multi-dimensional profile likelihood calculations
## S3 method for class 'ctsmTMB.fit' profile( fitted, parlist, grid.size = rep(10, length(parlist)), grid.qnt = rep(3, length(parlist)), hessian = FALSE, silent = FALSE, control = list(trace = 0, iter.max = 1000, eval.max = 1000), ... )## S3 method for class 'ctsmTMB.fit' profile( fitted, parlist, grid.size = rep(10, length(parlist)), grid.qnt = rep(3, length(parlist)), hessian = FALSE, silent = FALSE, control = list(trace = 0, iter.max = 1000, eval.max = 1000), ... )
fitted |
a ctmsTMB fit object |
parlist |
a named-list of parameters to profile over. The user can either supply grid-values in the list or leave it empty. If the any one list is empty then grid-values will be calculated using the estimated parameter mean value and standard deviation. |
grid.size |
a vector of |
grid.qnt |
a vector of |
hessian |
a boolean indicating whether to use the hessian or not during the profile optimization. |
silent |
boolean whether or not to mute current iteration number
the |
control |
a list of optimization output controls (see nlminb) |
... |
various arguments (not in use) |
The implementation was modified from that of https://github.com/kaskr/adcomp/blob/master/TMB/R/tmbprofile.R
library(ctsmTMB) model <- ctsmTMB$new() # create model model$addSystem(dx ~ theta * (mu+u-x) * dt + sigma_x*dw) model$addObs(y ~ x) model$setVariance(y ~ sigma_y^2) model$addInput(u) model$setParameter( theta = c(initial = 1, lower=1e-5, upper=50), mu = c(initial=1.5, lower=0, upper=5), sigma_x = c(initial=1, lower=1e-10, upper=30), sigma_y = 1e-2 ) model$setInitialState(list(1,1e-1)) # fit model to data fit <- model$estimate(Ornstein) # calculate profile likelihood out <- profile(fit,parlist=list(theta=NULL))library(ctsmTMB) model <- ctsmTMB$new() # create model model$addSystem(dx ~ theta * (mu+u-x) * dt + sigma_x*dw) model$addObs(y ~ x) model$setVariance(y ~ sigma_y^2) model$addInput(u) model$setParameter( theta = c(initial = 1, lower=1e-5, upper=50), mu = c(initial=1.5, lower=0, upper=5), sigma_x = c(initial=1, lower=1e-10, upper=30), sigma_y = 1e-2 ) model$setInitialState(list(1,1e-1)) # fit model to data fit <- model$estimate(Ornstein) # calculate profile likelihood out <- profile(fit,parlist=list(theta=NULL))
Basic summary of ctsmTMB fit object
## S3 method for class 'ctsmTMB.fit' summary(object, correlation = FALSE, ...)## S3 method for class 'ctsmTMB.fit' summary(object, correlation = FALSE, ...)
object |
a ctsmTMB fit object |
correlation |
boolean indicating whether or not to display the parameter correlation structure |
... |
additional arguments |
a summary of the estimated ctsmTMB model fit
library(ctsmTMB) model <- ctsmTMB$new() # create model model$addSystem(dx ~ theta * (mu+u-x) * dt + sigma_x*dw) model$addObs(y ~ x) model$setVariance(y ~ sigma_y^2) model$addInput(u) model$setParameter( theta = c(initial = 1, lower=1e-5, upper=50), mu = c(initial=1.5, lower=0, upper=5), sigma_x = c(initial=1, lower=1e-10, upper=30), sigma_y = 1e-2 ) model$setInitialState(list(1,1e-1)) # fit model to data fit <- model$estimate(Ornstein) # print model summary summary(fit, correlation=TRUE)library(ctsmTMB) model <- ctsmTMB$new() # create model model$addSystem(dx ~ theta * (mu+u-x) * dt + sigma_x*dw) model$addObs(y ~ x) model$setVariance(y ~ sigma_y^2) model$addInput(u) model$setParameter( theta = c(initial = 1, lower=1e-5, upper=50), mu = c(initial=1.5, lower=0, upper=5), sigma_x = c(initial=1, lower=1e-10, upper=30), sigma_y = 1e-2 ) model$setInitialState(list(1,1e-1)) # fit model to data fit <- model$estimate(Ornstein) # print model summary summary(fit, correlation=TRUE)