imd: Integral Matching Design

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

Description

Get the design matrix, formated data and weights used in integral matching.

Usage

imd(o, op, x = NULL, params = NULL)

Arguments

o	An object of class ode (created via mak, plk, etc.).
op	An object of class opt (created via opt). Compatibility check with o is conducted.
x	Matrix of dimension m-x-(d+1) containing custom time and smoothed values of the processes used for the integral matching, see details. If NULL (default) a linear-interpolation smoother is employed.
params	List of vectors of initial parameter values. Default (NULL) translates to the zero-vector. If the ODE system is linear in the parameter vector and the boundary constraints on are not unusual, then params is ignored.

Details

The design matrix is as follows:

X = (\int_{t_{i-1}}^{t_{i}}{\frac{df}{dparam}(x(s), context_l * param) * context_l ds})_{i=1}^n

Here f is the vector field of the ODE-system, x is the smoothed process and params is (internally) scaled with scales in reg objects in o.

Similiarly the observations and weights are concatinateed as follows:

Y = (x_{t_i} - x_{t_{i-1}})_{i=1}^n

W = ((w_{t_i} + w_{t_{i-1}}) / 2)_{i=1}^n

The number of decreases in time in x must match the number of decreases in time in y in op. The process x is simply a linear interpolation of x, hence finer discretisations give more refined integral matching. Each context is handled seperately. Moreover, the compatibility checks in rodeo are also conducted here.

Value

A list with:

X	List of design matrices, one for each parameter argument.
Y	A vector of observations.
W	A vector of weights.
X0	A matrix of initial states, one for each context. Interpolated from x.
xout	A matrix containing the values of the process at the time points in y, interpolated from x. Use it to check that x works as intended.

See Also

aim, rodeo.aim, numint

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, 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(2.1, 2.25, 1.5, 0, 0); x0 <- c(E = 8, S = 10, ES = 1.5, P = 1.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))

# Get the design matrix used in integral matching
d <- imd(m, opt(y))
head(d$X[[1]])

# Compare with glmnet
lasso <- glmnet::glmnet(x = d$X[[1]], y = d$Y, intercept = FALSE, lower.limits = 0)
a <- aim(m, opt(y, nlambda = 100), adapts = "scales")
all.equal(lasso$beta, a$params$rate)

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