bootstrap {KGode}

The 'bootstrap' function

Description

This function is used to perform bootstrap procedure to estimate parameter uncertainty.

Usage

bootstrap(kkk, y_no, ktype, K, ode_par, intp_data, www = NULL)

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.

K

the number of bootstrap replicates to collect.

ode_par

a vector of ode parameters estimated using gradient matching.

intp_data

a list of interpolations produced by gradient matching for each ode state.

www

an optional warping object (if warping has been performed using warpfun).

Details

Arguments of the 'bootstrap' function are 'ode' class, noisy observation, kernel type, the set of parameters that have been estimated before using gradient matching, a list of interpolations for each of the ode state from gradient matching, and the warping object (if warping has been performed). It returns a vector of the median absolute standard deviations for each ode state, computed from the bootstrap replicates.

Value

return a vector of the median absolute deviation (MAD) 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

## Perform bootstrap procedure to estimate the median absolute deviations of ode parameters
# here we get the resulting interpolation from gradient matching using 'rkg' for each ode state
bbb = rkgres$bbb
nst = length(bbb)
intp_data = list()
for( i in 1:nst) {
    intp_data[[i]] = bbb[[i]]$predictT(bbb[[i]]$t)$pred
}
K = 12 # the number of bootstrap replicates
mads = bootstrap(kkk, y_no, ktype, K, ode_par, intp_data)

## show the results of ode parameter estimation and its uncertainty
## using the standard gradient matching
ode_par
mads

############# gradient matching + ODE regularisation
crtype='i'
lam=c(10,1,1e-1,1e-2,1e-4)
lamil1 = crossv(lam,kkk,bbb,crtype,y_no)
lambdai1=lamil1[[1]]
res = third(lambdai1,kkk,bbb,crtype)
oppar = res$oppar

### do bootstrap here for gradient matching + ODE regularisation
ode_par = oppar
K = 12
intp_data = list()
for( i in 1:nst) {
    intp_data[[i]] = res$rk3$rk[[i]]$predictT(bbb[[i]]$t)$pred
}
mads = bootstrap(kkk, y_no, ktype, K, ode_par, intp_data)
ode_par
mads

############# gradient matching + ODE regularisation + warping
###### warp state
peod = c(6,5.3) #8#9.7     ## the guessing period
eps= 1          ## the standard deviation of period
fixlens=warpInitLen(peod,eps,rkgres)
kkkrkg = kkk$clone()
www = warpfun(kkkrkg,bbb,peod,eps,fixlens,y_no,kkkrkg$t)

### do bootstrap here for gradient matching + ODE regularisation + warping
nst = length(bbb)
K = 12
ode_par = www$wkkk$ode_par
intp_data = list()
for( i in 1:nst) {
    intp_data[[i]] = www$bbbw[[i]]$predictT(www$wtime[i, ])$pred
}
mads = bootstrap(kkk, y_no, ktype, K, ode_par, intp_data,www)
ode_par
mads


## End(Not run)