# demo/OldDemo/FhNEx.R In CollocInfer: Collocation Inference for Dynamic Systems

```##########################################
##       FitzHugh-Nagumo Example       ###
##########################################

library("CollocInfer")

#########################################################
####   Data generation from a solution of the ODE #######
#########################################################

data(FhNdata)

#  define the functions using function make.fhn

#  the functions that make.fhn defines:
#        fn       =  fhn.fun,
#        fn.ode  = fhn.fun.ode,
#        dfdx    = fhn.dfdx,
#        dfdp    = fhn.dfdp,
#        d2fdx2  = fhn.d2fdx2,
#        d2fdxdp = fhn.d2fdxdp

FhN = make.fhn()

# We can also remove the 'R' part of the data (ie, assume we did not measure
# it) and set the corresponding coefficients to zero

V.FhNdata = FhNdata
V.FhNdata[,2] = NA

# Now we want to set up a basis expansion; this makes use of
# routines in the fda library

knots  = seq(0,20,0.5)
norder = 3
nbasis = length(knots) + norder - 2
range  = c(0,20)

FhNbasis = create.bspline.basis(range=range,nbasis=nbasis,
norder=norder,breaks=knots)

# We'll start off by creating a smooth of each state variable to get initial
# values for the coefficients.

fdnames=list(NULL,c('V','R'),NULL)
bfdPar = fdPar(FhNbasis,lambda=1,int2Lfd(1))
DEfd0 = smooth.basis(FhNtimes,FhNdata,bfdPar,fdnames=fdnames)\$fd

#  plot the solution

plotfit.fd(FhNdata, FhNtimes, DEfd0)

#  extract the coefficients and assign variable names

coefs0 = DEfd0\$coefs
colnames(coefs0) = FhNvarnames

#  In the case of only V being measured, we can get coefficients for R,
#  so we set these to zero as starting values

V.coefs0 = coefs0
V.coefs0[,2] = 0

#  Set a value for lambda

lambda = 1000

# Optimize the coefficients  in coefs0, that is, run the lowest level or
# inner optimization loop, maximizing function J(c|theta,lambda)

Ires1	= Smooth.LS(FhN,FhNdata,FhNtimes,FhNpars,coefs0,FhNbasis,lambda,
in.meth='nlminb')

# We have some handy plotting routines

plotout1 = CollocInferPlots(Ires1\$coefs,FhNpars,Ires1\$lik,Ires1\$proc,FhNtimes,FhNdata)

# Now we can do the profiling, running the outer loop to optimize
#  parameter values.  You may wish to turn off output buffering in the
#  Misc menu in the command window at this point.

Ores1 = Profile.LS(FhN,FhNdata,FhNtimes,FhNpars,coefs=Ires1\$coefs,
basisvals=FhNbasis,lambda=lambda)

# And look at the result

Ores1\$pars

plotout2 = CollocInferPlots(Ores1\$coefs,Ores1\$pars,Ores1\$lik,Ores1\$proc,FhNtimes,FhNdata)

#  Repeat these analyses for the data with only V measured

Ores1.V = Profile.LS(FhN,V.FhNdata,FhNtimes,FhNpars,coefs=V.coefs0,
basisvals=FhNbasis,lambda=lambda)

Ores1.V\$pars

plotout2.V = CollocInferPlots(Ores1.V\$coefs,Ores1.V\$pars,Ores1.V\$lik,Ores1.V\$proc,FhNtimes,V.FhNdata)

# If we only coded the FitzHugh Nagumo function with no derivatives we would just
# have make.fhn()\$fn, in which case we default to estimating the derivatives by
# finite differencing.  The first argument provides access only to the
# function evaluating the right hand side.

fhn.fun <- function(times, y, p, more) {
r = y
r[, "V"] = p["c"] * (y[, "V"] - y[, "V"]^3/3 + y[, "R"])
r[, "R"] = -(y[, "V"] - p["a"] + p["b"] * y[, "R"])/p["c"]
return(r)
}

Ires1a = Smooth.LS(fhn.fun,FhNdata,FhNtimes,FhNpars,coefs0,FhNbasis,lambda,
in.meth='nlminb')

# Now we can do the profiling

Ores1a = Profile.LS(fhn.fun,FhNdata,FhNtimes,FhNpars,coefs=Ires1a\$coefs,
basisvals=FhNbasis,lambda=lambda,out.meth='nlminb')

Ores1a\$pars

plotout1a = CollocInferPlots((Ores1a\$coefs,Ores1a\$pars,Ores1a\$lik,Ores1a\$proc,FhNtimes,V.FhNdata)

#### Option 2:  set-up functions for proc and lik objects and then call
# optimizers.  This is a longer form of the analysis, in which the
#  lik and proc objects defining the first and second terms of
#  J(c|theta,lambda) are specified.  This is done in function LS.setup.
#  Then we can optimize initial parameter values using function
#  ParsMatchOpt

profile.obj = LS.setup(pars=FhNpars,fn=FhN,lambda=lambda,times=FhNtimes,
coefs=coefs1,basisvals=FhNbasis)
lik  = profile.obj\$lik
proc = profile.obj\$proc

# Now we can get initial parameter estimates from "gradient matching"
#  using function ParsMatchOpt

Pres  = ParsMatchOpt(FhNpars,coefs1,proc)
pars1 = Pres\$pars

pars1

# Smoothing can be done more specifically with 'inneropt'

Ires2 = inneropt(FhNdata,times=FhNtimes,pars1,coefs1,lik,proc)

# And we can also do profiling

Ores2 = outeropt(FhNdata,FhNtimes,pars1,coefs1,lik,proc)

Ores2\$pars

DEfd2 = fd(Ores2\$coefs,FhNbasis,fdnames)

plotfit.fd(FhNdata, FhNtimes, DEfd2)

# Option 3:  Set everything up manually.  This is the really long way to
# do the analysis, in which the error sum of squares fit is defined
# manually for both terms of J

## lik object

lik2 = make.SSElik()                      # we are using squared error
lik2\$bvals = eval.basis(FhNtimes,FhNbasis)  # values of the basis at times
lik2\$more = make.id()                     # identity transformation
lik2\$more\$weights = c(1,0)                # only use data for V

## proc object

qpts = knots[1:(length(knots)-1)]+diff(knots)/2  # Collocation points at
# midpointsbetween knots

proc2 = make.SSEproc()                    # we are using squared error
proc2\$bvals = list(bvals = eval.basis(qpts,FhNbasis),    # values and derivative of
dbvals = eval.basis(qpts,FhNbasis,1))  # basis at collocation points
proc2\$more = make.fhn()                   # FitzHugh-Nagumo right hand side
proc2\$more\$names = FhNvarnames            # State variable names
proc2\$more\$parnames = FhNparnames         # Parameter names
proc2\$more\$qpts = qpts                    # Collocation points
proc2\$more\$weights  = c(1000,1000)        # Weight relative to observations of
# each collocation point.

# We can also specify optimization methods

in.meth = 'nlminb'            # Inner Optimization
control.in = list(rel.tol=1e-12,iter.max=1000,eval.max=2000,trace=0)  # Optimization control

out.meth = 'optim'            # Outer Optimization
control.out = list(trace=6,maxit=100,reltol=1e-8,meth='BFGS') # Optimization control

# We can now fix the smooth or 'V' and try and choose 'R' to make the differential
# equation match as well as possible.

Vres0 = FitMatchOpt(coefs=V.coefs0,which=2,pars=FhNpars,proc2)
V.coefs1 = Vres0\$coefs

# And we can call the same inner optimization as above, this time with our
# own control parameters and method

Ires3 = inneropt(V.FhNdata,FhNtimes,FhNpars,V.coefs1,lik2,proc2,
in.meth=in.meth,control.in=control.in)

# And we can also do profiling with specified controls

Ores3.V = outeropt(V.FhNdata,FhNtimes,FhNpars,coefs=V.coefs1,lik=lik2,proc=proc2,
in.meth=in.meth,out.meth=out.meth,
control.in=control.in,control.out=control.out)

Ores3.V\$pars

DEfd3.V = fd(Ores3.V\$coefs,FhNbasis,fdnames)

plot(DEfd3.V['V'])
points(FhNtimes, V.FhNdata[,1])

plot(DEfd3.V['R'])

# We can also use finite differencing to calculate derivatives of the right hand side of the
# FitzHugh-Nagumo equations, this is defined by modifying the proc object.

proc2a = make.SSEproc()                    # we are using squared error
proc2a\$bvals = list(bvals = eval.basis(qpts,FhNbasis),    # values and derivative of
dbvals = eval.basis(qpts,FhNbasis,1))  # basis at collocation points
proc2a\$more = make.findif.ode()                   # FitzHugh-Nagumo right hand side
proc2a\$more\$names = FhNvarnames            # State variable names
proc2a\$more\$parnames = FhNparnames         # Parameter names
proc2a\$more\$qpts = qpts                    # Collocation points
proc2a\$more\$weights  = c(1000,1000)        # Weight relative to observations of
# each collocation point.
proc2a\$more\$more = list(fn = make.fhn()\$fn, eps=1e-6) # Tell findif that the function
# to difference is fhn and the
# difference stepsize

Ires3a = inneropt(V.FhNdata,FhNtimes,FhNpars,Ores3.V\$coefs,lik2,proc2a,
in.meth=in.meth,control.in=control.in)

# And we can also do profiling with specified controls

Ores3a = outeropt(V.FhNdata, FhNtimes, FhNpars, coefs=Ores3.V\$coefs,
lik=lik2, proc=proc2a,
in.meth=in.meth, out.meth=out.meth,
control.in=control.in, control.out=control.out)
```

## 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.