# demo/RMEx.r In CollocInfer: Collocation Inference for Dynamic Systems

```# This demonstration file is intended to accompany the CollocInfer manual.

library(CollocInfer)

# To ensure reproducibility
set.seed((2004*2007)/2014)

# This file carrys out the analysis for the Rosenzweig-MacArthur model that is
# given in the manual. We assume that the reader is familiar with the FitzHugh-Nagumo
# example that is first developed in the manual, so we will focus on the elaborations
# that this model requires.

# First, we write down the three-species Rosenzweig-Macarthur model in a form
# suitable for CollocInfer

RosMac2 = function(t,x,p,more){

p = exp(p)
dx = x

dx[,'C1'] = p['rho1']*x[,'C1']*(1- x[,'C1']/p['kappaC1']- x[,'C2']/p['kappaC2']) - p['pi']*p['gamma']*x[,'C1']*x[,'B']/(p['kappaB']+p['pi']*x[,'C1']+x[,'C2'])
dx[,'C2'] = p['rho2']*x[,'C2']*(1- x[,'C1']/p['kappaC1']- x[,'C2']/p['kappaC2']) - p['gamma']*x[,'C2']*x[,'B']/(p['kappaB']+p['pi']*x[,'C1']+x[,'C2'])
dx[,'B'] = p['chi']*p['gamma']*(p['pi']*x[,'C1']+x[,'C2'])*x[,'B']/(p['kappaB']+p['pi']*x[,'C1']+x[,'C2']) - p['delta']*x[,'B']

return(dx)
}

# We will define some interesting parameters

RMpars = c(0.2,0.025,0.125,2.2e4,1e5,5e6,1,1e9,0.3)
RMParnames = c('pi','rho1','rho2','kappaC1','kappaC2','gamma','chi','kappaB','delta')

# Which we represent on the log scale

logpars = log(RMpars)
names(logpars) = RMParnames

# And we also need some initial conditions (with named entries)

RMVarnames = c('C1','C2','B')

x0 = c(50,50,2)
names(x0) = RMVarnames

# The following version of RosMac2 is suitable for use with lsoda

RosMac2ODE = function(t,z,p){
p = exp(p)
x = exp(z)
dx = x

dx['C1'] = p['rho1']*x['C1']*(1- x['C1']/p['kappaC1']-x['C2']/p['kappaC2']) - p['pi']*p['gamma']*x['C1']*x['B']/(p['kappaB']+p['pi']*x['C1']+x['C2'])
dx['C2'] = p['rho2']*x['C2']*(1- x['C2']/p['kappaC2']- x['C1']/p['kappaC1']) - p['gamma']*x['C2']*x['B']/(p['kappaB']+p['pi']*x['C1']+x['C2'])
dx['B'] = p['chi']*p['gamma']*(p['pi']*x['C1']+x['C2'])*x['B']/(p['kappaB']+p['pi']*x['C1']+x['C2']) - p['delta']*x['B']

return(list(dx/x))
}

# With this we can solve the ODE at 200 successive days
time = 0:200
res0 = lsoda(log(x0),time,RosMac2ODE,p = logpars)

# and plot the solutions
matplot(exp(res0[,2:4]),type='l')

# We'll obtain data by adding noise
data = res0[,2:4] + 0.2*matrix(rnorm(603),201,3)

# and name the columns
colnames(data) = RMVarnames

# Giving us the following plot
matplot(data,cex.lab=2.5,cex.axis=2.5,cex=1.5,xlab='days',pch = c('1','2','B'))

# Now we need to set up the CollocInfer machinery

# First we'll define a basis with knots each time point

rr = range(time)
knots = seq(rr,rr,by=1)

bbasis = create.bspline.basis(rr,norder=4,breaks=knots)

# And obtain an initial set of parameters and coefficients from smoothing

coef0 = smooth.basis(time,data,fdPar(bbasis,int2Lfd(2),10))\$fd\$coef
colnames(coef0) = RMVarnames

# We will now create the profiling objects, but use the log transformation by
# setting posproc=TRUE

out = LS.setup(pars=logpars,coefs=coef0,basisvals=bbasis,fn=RosMac2,lambda=1e5,
times=time,posproc=TRUE)

lik = out\$lik
proc = out\$proc

# We'll do gradient matching to get parameter estimates corresponding to this
# smooth

res1 = ParsMatchOpt(logpars,coef0,proc)

# And now profiling
res3 = outeropt(data,time,res1\$pars,coef0,lik,proc)

# Let's have a look at the parameters that we got
exp(res3\$pars)

# and the agreement with data and model
out3 = CollocInferPlots(res3\$coefs,res3\$pars,lik,proc,times=time,data=data,
cex.lab=2.5,cex.axis=2.5,cex=1.5,lwd=3)

## Now we will compliate the model. In reality, we only observe the sum of C1
# and C2.

data2 = cbind( log( exp(data[,'C1'])+exp(data[,'C2'])), data[,'B'])
matplot(data2,cex.lab=2.5,cex.axis=2.5,pch=c('C','B'),cex=1.5,xlab='days')

# To deal with this, we need to define a transformation function from our
# state variables to the expected observations. In this case (since the states
# are on the log scale) we exponentiate to get back to the original scale, add
# C1 and C2 and then take the log again.
RMobsfn = function(t,x,p,more)
{
x = exp(x)
y = cbind( x[,'C1']+x[,'C2'],x[,'B'])
return(log(y))
}

# We can now create new profiling objects that incorporate this transformation
# function by specfying likfn

out = LS.setup(pars=logpars,coefs=coef0,basisvals=bbasis,fn=RosMac2,lambda=1e5,
times=time,posproc=TRUE,likfn=RMobsfn)

lik2 = out\$lik
proc2 = out\$proc

# To see how we perform in this situation, we'll start by setting the two
# columns of the data smooth to zero

coef02 = coef0
coef02[,1:2] = 0

# We'll now pull these columns into line with the rotifer column (which we can
# still smooth.

Fres3 = FitMatchOpt(coef02,1:2,res1\$pars,proc2)

# And we'll run profiling and have a look at what we get.
res32 = outeropt(data2,time,res1\$pars,Fres3\$coefs,lik2,proc2)
exp(res32\$pars)
out32 = CollocInferPlots(res32\$coefs,res32\$pars,lik2,proc2,times=time,data=data2,
datanames=c('C','B'),cex.lab=2.5,cex.axis=2.5,cex=1.5,lwd=3)

# The section at the end of this demo goes through setting up lik2 and proc2
# manually rather than through LS.setup.

### In this framework is is not unreasonable to expect that we have
# repeated experiments. When these are very well structured and all have common
# observation times, this can be easily accommodated in CollocInfer (it can
# accommodate less regular replicated experiments, but requires work to set things
# up manuall.

# We'll create a second experiment starting from new initial conditions

x03 = c(15,25,4)
names(x03) = RMVarnames

res03 = lsoda(log(x03),time,RosMac2ODE,p = logpars)

data03 =  res03[,2:4] + 0.2*matrix(rnorm(603),201,3)

# and set up a three dimensional array in which the experiment number is the
# second dimension and the third dimension is the variable being measured (this
# is to agree with conventions in the fda package)

alldat = array(0,c(201,2,3))
alldat[,1,] = data
alldat[,2,] = data03

# Nowe we'll smooth the second experiment

coef3 = smooth.basis(time,data03,fdPar(bbasis,int2Lfd(2),10))\$fd\$coef

# and create a three-dimensional array with all the coefficients together

coefs = array(0,c(dim(coef3),2,3))
coefs[,1,] = coef0
coefs[,2,] = coef3

# These three dimensional arrays can be given to LSsetup which understands
# thi structure and knows what to do with it.

out = LS.setup(pars=logpars,coefs=coefs,basisvals=bbasis,fn=RosMac2,lambda=1e5,
times=time,data=alldat,posproc=TRUE,names=RMVarnames)

lik3 = out\$lik
proc3 = out\$proc

# We can now call the inner optimisation to use a model-based smooth
res13 = inneropt(data=out\$data,times=out\$times,pars=res1\$pars,coefs=out\$coefs,lik=lik3,proc=proc3)

# And use profiling to estimate parameters
res33 = outeropt(data=out\$data,times=out\$times,res1\$pars,res13\$coefs,lik3,proc3)

# These parameters should hopefully be closer to the truth than with only one experimental run
exp(res33\$pars)

# And we can also examine the fit to the data and model. In this case, the times vector
# wraps around creating a few unpleasant graphical effects.

out3 = CollocInferPlots(res33\$coefs,res33\$pars,lik3,proc3,times=out\$times,data=out\$data,datanames=c('B','C'),
cex.lab=2.5,cex.axis=2.5,cex=1.5,lwd=3)

###### Manual set-up

# Here we create lik2 and proc2 (corresponding to the indirectly observed
# single-run experiment above in order to demonstrate their structure.

# First we need to specify the matrices of basis values that we will use

# at observation time points
bvals.obs = eval.basis(time,bbasis)

# and quadrature times, these are midpoints between knots plus the end points
# The quadrature weights are all equal, but we have multiplied them by the lambda
# that we are using, in this case 1e5.

qpts = c(knots,knots[1:(length(knots)-1)]+diff(knots)/2,knots[length(knots)])
qwts = 1e5*matrix(1,length(qpts),3)/length(qpts)

# basis values for proc is a list

dbvals = eval.basis(qpts,bbasis,1))

# Now we create the lik object

# make.SSElik() sets up the squared error criteria
lik.m = make.SSElik()

# Attach the values of the basis expansion at the observation time points
lik.m\$bvals = bvals.obs

# We need to specify the transformation of the state variables that is to be
# compared with the data. In this case, we will use finite differencing to
# to compute the derivatives that we need; we can achieve this by first
# employing make.findif.ode()
lik.m\$more = make.findif.ode()

# and then telling findif.ode that the function it is finite differencing is
# RMobsfn
lik.m\$more\$more = list(fn = RMobsfn,eps=1e-6,more=NULL)

# We also need to give lik.m a set of weights. This has to occur in
# the more element of lik.m because it is used inside lik.m\$fn.
lik.m\$more\$weights = array(1,dim(data))

# We can also create the proc object manually. First we call make.SSEproc()
# in order to set up the squared error functions
proc.m = make.SSEproc()

# We also specify the basis values and their derivatives at the quadrature points

# Now we need to tell SSEproc about the right hand side of the ODE and its
# derivatives. Here we will use finite differencing again, as in lik.m:
proc.m\$more = make.findif.ode()

# In fact, we are going to finite difference the right hand side of the
# the ODE for the log transformed data. Here we specify the log transform
proc.m\$more\$more = list(fn = make.logtrans()\$fn,eps=1e-6)

# and then give it the (non log-transform) Rosenzweig-MacArthur equations
proc.m\$more\$more\$more\$fn = RosMac2

# proc.m\$more also needs to include some elements for internal processing. In
# particular the following specify the quadrature points and weights
proc.m\$more\$weights = qwts
proc.m\$more\$qpts = qpts

# We will also tell it about the parameter and variable names.
proc.m\$more\$parnames = RMParnames
proc.m\$more\$names = RMVarnames

# And let's check that this all works

Fres.m3 = FitMatchOpt(coef02,1:2,res1\$pars,proc.m)
res.m32 = outeropt(data2,time,res1\$pars,Fres3\$coefs,lik2,proc.m)
exp(res.m32\$pars)
out.m32 = CollocInferPlots(res.m32\$coefs,res.m32\$pars,lik.m,proc.m,times=time,data=data2)

# If you compare this to out32, you should have the same estimates.
```

## Try the CollocInfer package in your browser

Any scripts or data that you put into this service are public.

CollocInfer documentation built on May 2, 2019, 4:03 a.m.