Package: RKGx_ode Type: Package Title: Reproducing kernel Hilbert Space and gradient matching based estimation of parameters of Differential equations Version: 1.0 Date: 2017-06-13 Author: Mu Niu , Simon Rogers , Maurizio Filippone , Dirk Husmeier ## Working Example ## user need to provide the differential equations. If they want to generate simulation data ## by solving the ode, they also need to provide the true value of ode parameters, ## initial value of system states and the time span of the simulation. ## ## ode model : dx/dt = f(x, theta) x is system state theta is ode parameter. ## example: lotka volterra model : 2 system states and 4 ode parameters ## dx1/dt = alpha * x1 - beta * x2 * x1 ## dx2/dt = -gamma * x2 + delta * x1*x2 ## ## true value of ode paraemter alpha = 1 beta =1 gamma =4 delta =1 ## initial states are x1 = 0.5 x2 = 1 ## time span of simulation c(0,6) ## if user do not want to use simulation data, they need to provdie observations for each system ## states, the time points for each observation and the ode equations. ## the current setting also require user provide a function df(x,theta)/dtheta which is ## LV_grlNODE in the example. library(pspline) ### define ode 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] ) ) } 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 } ### generate data by solving odes source('ode.r') kkk = ode$new(2,fun=LV_fun,grfun=LV_grlNODE) xinit = as.matrix(c(0.5,1)) tinterv = c(0,6) kkk$solve_ode(c(1,1,4,1),xinit,tinterv) ## add noise to the true solution y_no is the noisy data set.seed(19573) n_o = max( dim( kkk$y_ode) ) noise = 0.01 ## variance of noise y_no = t(kkk$y_ode) + rmvnorm(n_o,c(0,0),noise*diag(2)) ################################ standard gradient matching ############ source('kernel1.r') source('rkhs1.r') ## 1st system state interpolation ann1 = RBF$new(1) bbb1 = rkhs$new(t(y_no)[1,],kkk$t,rep(1,n_o),1,ann1) bbb1$skcross(c(5) ) resb1 = bbb1$predict() plot(kkk$t,resb1$pred) ## 2nd system state interpolation ann2 = RBF$new(1) bbb2 = rkhs$new(t(y_no)[2,],kkk$t,rep(1,n_o),1,ann2) bbb2$skcross(c(5)) resb2 = bbb2$predict() plot(kkk$t,resb2$pred) ## parameter estimation intp= rbind(resb1$pred,resb2$pred) grad= rbind(resb1$grad,resb2$grad) kkk$optim_par( c(0.1,0.1,0.1,0.1), intp, grad ) ################################## use warping to improve parameter estimation ############### source('WarpSin.r') ################################ warp 1st state ############################### p0=6 ## a guess of the period of warped signal eps= 1 ## the standard deviation of period lambda_t= 50 ## the weight of fixing the end of interval y_c = resb1$pred ##y_use[1,] ################# learn wapring function ############### fixlen = 3 ### set the starting value of hyper par. ## you can also use bbb$slowWarp(2,6,p0,eps) to find good starting value. wsigm = Sigmoid$new(1) bbbs = Warp$new( y_c,kkk$t,rep(1,n_o),lambda_t,wsigm) ppp = bbbs$warpSin( fixlen, 10 ) ############### smooth version of the warping function ### t_me= bbbs$tw ben = MLP$new(c(5,5)) rkben = rkhs$new(t(t_me),kkk$t,rep(1,n_o),1,ben) rkben$mkcross(c(5,5)) resm = rkben$predict() ################# relearn interpolation in warped time domain ### ann1w = RBF$new(1) bbb1w = rkhs$new(t(y_no)[1,],resm$pred,rep(1,n_o),1,ann1w) bbb1w$skcross(5) resb1 = bbb1w$predict() ################################# warp 2nd state ############################## p0=5.3 ## the guessing period eps= 1 ## the standard deviation of period lambda_t= 50 ## the weight of fixing the end of interval y_c2 = resb2$pred ####################### learn wapring function ################# fixlen2 = 3 wsigm2 = Sigmoid$new(1) bbbs2 = Warp$new( y_c2,kkk$t,rep(1,n_o),lambda_t,wsigm2) ppp2 = bbbs2$warpSin( fixlen2, 10 ) t_me2= bbbs2$tw ben2 = MLP$new(c(5,5)) rkben2 = rkhs$new(t(t_me2),kkk$t,rep(1,n_o),1,ben2) rkben2$mkcross(c(5,5)) resm2 = rkben2$predict() ann2w = RBF$new(1) bbb2w = rkhs$new(t(y_no)[2,],resm2$pred,rep(1,n_o),1,ann2w) bbb2w$skcross(5) resb2w = bbb2w$predict() ############### gradient Matching for warped signal and estimating ode parameters #### intp= rbind(resb1$pred,resb2w$pred) grad= rbind(resb1$grad*resm$grad,resb2w$grad*resm2$grad) kkk$optim_par( c(0.1,0.1,0.1,0.1), intp, grad )