episode (version 1.0.0)

rodeo.ode: Regularised Ordinary Differential Equation Optimisation (RODEO)

Description

Fit the parameters (and optionally initial states) for a ordinary differential equation model with data sampled across different contexts.

Usage

# S3 method for ode
rodeo(x, op, x0, params, trace = FALSE, ...)

Arguments

x

ode-object created via mak, plk etc.

op

opt-object created via opt-function (compatibility check with x is conducted).

x0

A vector (or matrix) of non-negative initialisation(s) for the initial state in the optimisation problem, see details.

params

A list of initialisations for the parameter arguments in the optimisation problem, see details. A list of matrices if multiple initialisations are desired. Sparse matrix "dgCMatrix" is allowed.

trace

Logical indicating if status messages should be printed during rodeo.

...

Additional arguments passed to rodeo.

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.

Details

For running a single initialisation of the optimisation procedure, supply x0 as vector (where the initial states for the contexts are concatinated) and params as a list with a vector entry for each parameter. The resulting estimates are stored in matrices, with each column representing a lambda-value. The convergence codes, steps and losses are stored as vectors, one entry for each value of lambda.

For running multiple initialisations, supply x0 as matrix with the individual initialisations as columns (each column vector as described above) and params as list of matrices with initialisations stored column-wise.

The initial state estimates (for the different lambda-values) are returned as a matrix (array) (row = coordinate, column = lambda-value, (slice = initialisation)) and the parameter estimates are returned as a list of (lists with) sparse matrices (row = coordinate, column = lambda-value). The convergence codes, steps and losses are stored as matrices (row = lambda, column = initialisation).

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

If explicitly setting x0 to NULL, the system uses the first observation of each context (throws error if non finites).

If explicitly setting params to NULL, 0 initialisations are used.

See Also

rodeo, rodeo.aim

Examples

# NOT RUN {
set.seed(123)
# Example: Michaelis-Menten system with two 0-rate reactions
A <- matrix(c(1, 1, 0, 0,
              0, 0, 1, 0,
              0, 0, 1, 0,
              0, 1, 0, 0,
              0, 0, 0, 1), ncol = 4, byrow = TRUE)
B <- matrix(c(0, 0, 1, 0,
              1, 1, 0, 0,
              1, 0, 0, 1,
              0, 0, 0, 1,
              0, 0, 1, 0), ncol = 4, byrow = TRUE)
k <- c(1, 2, 0.5, 0, 0); x0 <- c(E = 2, S = 8, ES = 0.5, P = 0.5)
Time <- seq(0, 10, by = 1)

# Simulate data, in second context the catalytic rate has been doubled
m <- mak(A, B)
contexts <- cbind(1, c(1, 1, 2, 1, 1))
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 = .1), nrow = nrow(y))

# Fit data using rodeo on mak-object
op <- opt(y)
fit <- rodeo(m, op, x0 = NULL, params = NULL)

# Example: fit data knowing doubled catalytic rate
m_w_doubled <- mak(A, B, r = reg(contexts = contexts))
fit <- rodeo(m_w_doubled, op, x0 = NULL, params = NULL)

# }