imd: Integral Matching Design

Description

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

Usage

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

Arguments

An object of class ode (created via mak, plk, etc.).

An object of class opt (created via opt). Compatibility check with o is conducted.

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.

Value

A list with:

List of design matrices, one for each parameter argument.

A vector of observations.

A vector of weights.

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.

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.

Examples

Run this code

# NOT RUN {
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)

# }