rodeo: Regularised Ordinary Differential Equation Optimisation...

In episode: Estimation with Penalisation in Systems of Ordinary Differential Equations

Description

Fit the parameters (and optionally initial states) for an Ordinary Differential Equation model with data sampled across contexts. Two implementations exist:

rodeo.ode

The raw form of rodeo, based on ode (created via mak, plk, etc.) and opt-objects. If needed, everything can be user-specified through these two objects. However, parameter initialisations are always required and no action towards adjusting weights or scales of the parameters are taken.

rodeo.aim

The more automatic form of rodeo, where automatic parameter initialisations and (optional) adjustments of scales and weights are available through aim.

Usage

rodeo(x, ...)

Arguments

x	Object that specifies ode system.
...	Additional arguments passed to rodeo.

Details

For details on the loss function, optimisation, etc. See documentation of opt.

Value

An object with S3 class "rodeo":

o	Original ode-object.
op	Original opt-object with default values for lambda_min_ratio, lambda and (if needed) a inserted, if these were originally NULL.
params	Parameter estimates, stored as list of sparse column format matrices, "dgCMatrix" (or a list of those if multiple initialisations). Rows represent coordinates and columns represent the lambda value.
x0s	Initial state estimates stored in a matrix (or array). Rows represent coordinates, columns represent the lambda value and (if multiple initialisations) slices represent initialisations.
dfs	A matrix (or array, if multiple initialisations) of degrees of freedom. Row represents a parameter (the first is always the initial state parameter), columns represent lambda, slices represent initialisation, if multiple are provided.
codes	A matrix (or array) of convergence codes organised as dfs. 0: The convergence criteria is met (see details in opt). Current estimate is probably a local minimum. 1: Backtracking in the last iteration yields no numerical improvement, but no unsual behavior observed. Current estimate is probably a local minimum. However, if exact_gradient = FALSE in the reg-object in the ode-object, changing this may improve the code. Alternatively one can adjust backtracking via backtrack_max and tau_min in reg objects in ode object. 2: The optimisation procedure exceeded maximal number of steps (step_max in reg objects). 3: The last gradient was unsually large. Either the tolerances in reg objects are off or the ODE systems is very sensitive and runs over long time spans. In the latter case, initialisation(s) may have inappropriate zeros (change initialisation and/or make sure they start at smaller lambda value). 4: The numeric ODE solver exceeded maximal number of steps. Check if supplied initial states were out of bounds, if not increase step_max (or tol) in reg-objects in ode-object.
steps	A matrix (or array) holding number of steps used in optimisation procedure. Organised as dfs.
losses	A vector (or matrix) of unpenalised losses at optimum for each lambda value (stored row-wise if multiple are provided).
penalties	A matrix (or array) of penalties for each parameter, organised as dfs.
jerr	A matrix (or array) of summary codes (for internal debugging), organised as dfs.

See Also

rodeo.aim, rodeo.ode

Examples

set.seed(123)
# Michaelis-Menten system with two 0-rate reactions
A <- matrix(c(1, 1, 0, 0,
 0, 0, 1, 0,
 0, 0, 1, 0,
 0, 0, 0, 1,
 1, 0, 0, 0), ncol = 4, byrow = TRUE)
B <- matrix(c(0, 0, 1, 0,
 1, 1, 0, 0,
 1, 0, 0, 1,
 0, 0, 1, 0,
 1, 0, 1, 0), ncol = 4, byrow = TRUE)
k <- c(1.1, 2.5, 0.25, 0, 0); x0 <- c(E = 2, S = 2, ES = 8.5, P = 2.5)
Time <- seq(0, 10, by = 1)

# Simulate data, in second context the catalytic rate has been doubled
contexts <- cbind(1, c(1, 1, 2, 1, 1))
m <- mak(A, B, r = reg(contexts = contexts))
y <- numsolve(m, c(Time, Time), cbind(x0, x0 + c(2, -2, 0, 0)), contexts * k)
y[, -1] <- y[, -1] + matrix(rnorm(prod(dim(y[, -1])), sd = .5), nrow = nrow(y))

# Example: fit data using rodeo on mak-object
op <- opt(y, nlambda = 10)
fit <- rodeo(m, op, x0 = NULL, params = NULL)
fit$params$rate

# Example: fit dat using rodeo on aim-object
a <- aim(m, op)
a$params$rate
fit <- rodeo(a)
fit$params$rate

episode documentation built on May 1, 2019, 11:17 p.m.