MagiSolver {magi}

MAnifold-constrained Gaussian process Inference (MAGI)

Description

Core function of the MAGI method for inferring the parameters and trajectories of dynamic systems governed by ordinary differential equations.

Usage

MagiSolver(y, odeModel, tvec, control = list())

Arguments

y

data matrix of observations

odeModel

list of ODE functions and inputs. See details.

tvec

vector of discretization time points corresponding to rows of y. If missing, MagiSolver will use the column named 'time' in y.

control

list of control variables, which may include 'sigma', 'phi', 'theta', 'xInit', 'mu', 'dotmu', 'priorTemperature', 'niterHmc', 'nstepsHmc', 'burninRatio', 'stepSizeFactor', 'bandSize', 'useFixedSigma', 'kerneltype', 'skipMissingComponentOptimization', 'positiveSystem', 'verbose'. See details.

Details

The data matrix y has a column for each system component, and optionally a column 'time' with the discretization time points. If the column 'time' is not provided in y, a vector of time points must be provided via the tvec argument. The rows of y correspond to the discretization set I at which the GP is constrained to the derivatives of the ODE system. To set the desired discretization level for inference, use setDiscretization to prepare the data matrix for input into MagiSolver. Missing observations are indicated with NA or NaN.

The list odeModel is used for specification of the ODE system and its parameters. It must include five elements:

fOde

function that computes the ODEs, specified with the form f(theta, x, tvec). fOde should return a matrix where columns correspond to the system components of x, see examples.

fOdeDx

function that computes the gradients of the ODEs with respect to the system components. fOdeDx should return a 3-D array, where the slice [, i, j] is the partial derivative of the ODE for the j-th system component with respect to the i-th system component, see examples.

fOdeDtheta

function that computes the gradients of the ODEs with respect to the parameters \theta. fOdeDtheta should return a 3-D array, where the slice [, i, j] is the partial derivative of the ODE for the j-th system component with respect to the i-th parameter in \theta, see examples.

thetaLowerBound

a vector indicating the lower bounds of each parameter in \theta.

thetaUpperBound

a vector indicating the upper bounds of each parameter in \theta.

Additional control variables can be supplied to MagiSolver via the optional list control, which may include the following:

sigma

a vector of noise levels (observation noise standard deviations) \sigma for each component, at which to initialize MCMC sampling. By default, MagiSolver computes starting values for sigma via Gaussian process (GP) smoothing. If the noise levels are known, specify sigma together with useFixedSigma = TRUE.

phi

a matrix of GP hyper-parameters for each component, with rows for the kernel hyper-parameters and columns for the system components. By default, MagiSolver estimates phi via an optimization routine.

theta

a vector of starting values for the parameters \theta, at which to initialize MCMC sampling. By default, MagiSolver uses an optimization routine to obtain starting values.

xInit

a matrix of values for the system trajectories of the same dimension as y, at which to initialize MCMC sampling. Default is linear interpolation between the observed (non-missing) values of y and an optimization routine for entirely unobserved components of y.

mu

a matrix of values for the mean function of the GP prior, of the same dimension as y. Default is a zero mean function.

dotmu

a matrix of values for the derivatives of the GP prior mean function, of the same dimension as y. Default is zero.

priorTemperature

the tempering factor by which to divide the contribution of the GP prior, to control the influence of the GP prior relative to the likelihood. Default is the total number of observations divided by the total number of discretization points.

niterHmc

MCMC sampling from the posterior is carried out via the Hamiltonian Monte Carlo (HMC) algorithm. niterHmc specifies the number of HMC iterations to run. Default is 20000 HMC iterations.

nstepsHmc

the number of leapfrog steps per HMC iteration. Default is 200.

burninRatio

the proportion of HMC iterations to be discarded as burn-in. Default is 0.5, which discards the first half of the MCMC samples.

stepSizeFactor

initial leapfrog step size factor for HMC. Can be a specified as a scalar (applied to all posterior dimensions) or a vector (with length corresponding to the dimension of the posterior). Default is 0.01, and the leapfrog step size is automatically tuned during burn-in to achieve an acceptance rate between 60-90%.

bandSize

a band matrix approximation is used to speed up matrix operations, with default band size 20. Can be increased if MagiSolver returns an error indicating numerical instability.

useFixedSigma

logical, set to TRUE if sigma is known. If useFixedSigma = TRUE, the known values of \sigma must be supplied via the sigma control variable. Default is FALSE.

kerneltype

the GP covariance kernel, generalMatern is the default and recommended choice. Other available choices are matern, rbf, compact1, periodicMatern. See calCov for their definitions.

skipMissingComponentOptimization

logical, set to TRUE to skip automatic optimization for missing components. If skipMissingComponentOptimization = TRUE, values for xInit and phi must be supplied for all system components. Default is FALSE.

positiveSystem

logical, set to TRUE if the system cannot be negative. Default is FALSE.

verbose

logical, set to TRUE to output diagnostic and progress messages to the console. Default is FALSE.

Value

MagiSolver returns an object of class magioutput which contains the following elements:

theta

matrix of MCMC samples for the system parameters \theta, after burn-in.

xsampled

array of MCMC samples for the system trajectories at each discretization time point, after burn-in.

sigma

matrix of MCMC samples for the observation noise SDs \sigma, after burn-in.

phi

matrix of estimated GP hyper-parameters, one column for each system component.

lp

vector of log-posterior values at each MCMC iteration, after burn-in.

y, tvec, odeModel

from the inputs to MagiSolver.

References

Wong, S. W. K., Yang, S., & Kou, S. C. (2024). 'magi': A Package for Inference of Dynamic Systems from Noisy and Sparse Data via Manifold-Constrained Gaussian Processes. *Journal of Statistical Software*, 109 (4), 1-47. doi:10.18637/jss.v109.i04

Yang, S., Wong, S. W. K., & Kou, S. C. (2021). Inference of Dynamic Systems from Noisy and Sparse Data via Manifold-constrained Gaussian Processes. *Proceedings of the National Academy of Sciences*, 118 (15), e2020397118. doi:10.1073/pnas.2020397118

Examples

# Set up odeModel list for the Fitzhugh-Nagumo equations
fnmodel <- list(
  fOde = fnmodelODE,
  fOdeDx = fnmodelDx,
  fOdeDtheta = fnmodelDtheta,
  thetaLowerBound = c(0, 0, 0),
  thetaUpperBound = c(Inf, Inf, Inf)
)

# Example noisy data observed from the FN system
data(FNdat)

# Set discretization for a total of 81 equally-spaced time points from 0 to 20
yinput <- setDiscretization(FNdat, by = 0.25)

# Run MagiSolver
# Short sampler run for demo only, more iterations needed for convergence
MagiSolver(yinput, fnmodel, control = list(nstepsHmc = 5, niterHmc = 101))

# Use 3000 HMC iterations with 100 leapfrog steps per iteration
FNres <- MagiSolver(yinput, fnmodel, control = list(nstepsHmc = 100, niterHmc = 3000))
# Summary of parameter estimates
summary(FNres)
# Plot of inferred trajectories
plot(FNres, comp.names = c("V", "R"), xlab = "Time", ylab = "Level")