rkg: The 'rkg' function

In mu2013/KGode: Kernel Based Gradient Matching for Parameter Inference in Ordinary Differential Equations

Description

This function is used to create 'rkhs' class object and estimate ode parameters using standard gradient matching.

Usage

rkg(kkk, y_no, ktype)

Arguments

kkk	ode class object.
y_no	matrix(of size n_s*n_o) containing noisy observations. The row(of length n_s) represent the ode states and the column(of length n_o) represents the time points.
ktype	character containing kernel type. User can choose 'rbf' or 'mlp' kernel.

Details

Arguments of the 'rkg' function are 'ode' class, noisy observation, and kernel type. It return the interpolation for each of the ode states. The Ode parameters are estimated using gradient matching, and the results are stored in the 'ode' class as the ode_par attribute.

Value

return list containing :

• intp - list containing interpolation for each ode state.

• bbb - rkhs class objects for each ode state.

Author(s)

Mu Niu [e-mail]

Examples

## Not run: 
require(mvtnorm)
noise = 0.1  ## set the variance of noise
SEED = 19537
set.seed(SEED)
## Define ode function, we use lotka-volterra model in this example. 
## we have two ode states x[1], x[2] and four ode parameters alpha, beta, gamma and delta.
LV_fun = function(t,x,par_ode){
  alpha=par_ode[1]
  beta=par_ode[2]
  gamma=par_ode[3]
  delta=par_ode[4]
  as.matrix( c( alpha*x[1]-beta*x[2]*x[1] , -gamma*x[2]+delta*x[1]*x[2] ) )
}
## Define the gradient of ode function against ode parameters 
## df/dalpha,  df/dbeta, df/dgamma, df/ddelta where f is the differential equation.
LV_grlNODE= function(par,grad_ode,y_p,z_p) { 
alpha = par[1]; beta= par[2]; gamma = par[3]; delta = par[4]
dres= c(0)
dres[1] = sum( -2*( z_p[1,]-grad_ode[1,])*y_p[1,]*alpha ) 
dres[2] = sum( 2*( z_p[1,]-grad_ode[1,])*y_p[2,]*y_p[1,]*beta)
dres[3] = sum( 2*( z_p[2,]-grad_ode[2,])*gamma*y_p[2,] )
dres[4] = sum( -2*( z_p[2,]-grad_ode[2,])*y_p[2,]*y_p[1,]*delta)
dres
}

## create a ode class object
kkk0 = ode$new(2,fun=LV_fun,grfun=LV_grlNODE)
## set the initial values for each state at time zero.
xinit = as.matrix(c(0.5,1))
## set the time interval for the ode numerical solver.
tinterv = c(0,6)
## solve the ode numerically using predefined ode parameters. alpha=1, beta=1, gamma=4, delta=1.
kkk0$solve_ode(c(1,1,4,1),xinit,tinterv) 

## Add noise to the numerical solution of the ode model and use it as the noisy observation.
n_o = max( dim( kkk0$y_ode) )
t_no = kkk0$t
y_no =  t(kkk0$y_ode) + rmvnorm(n_o,c(0,0),noise*diag(2))

## Create a ode class object by using the simulation data we created from the ode numerical solver.
## If users have experiment data, they can replace the simulation data with the experiment data.
## Set initial value of ode parameters.
init_par = rep(c(0.1),4)
init_yode = t(y_no)
init_t = t_no
kkk = ode$new(1,fun=LV_fun,grfun=LV_grlNODE,t=init_t,ode_par= init_par, y_ode=init_yode )

## The following examples with CPU or elapsed time > 10s

##Use function 'rkg' to estimate the ode parameters. The standard gradient matching method is coded
##in the the 'rkg' function. The parameter estimations are stored in the returned vector of 'rkg'.
## Choose a kernel type for 'rkhs' interpolation. Two options are provided 'rbf' and 'mlp'. 
ktype ='rbf'
rkgres = rkg(kkk,y_no,ktype)
## show the results of ode parameter estimation using the standard gradient matching
kkk$ode_par

## End(Not run)

mu2013/KGode documentation built on June 23, 2020, 8:04 p.m.