# inst/OldDemo/SEIREX.r In CollocInfer: Collocation Inference for Dynamic Systems

```##########
# This demonstration illustrates the use of log transformations with
# CollocInfer. We use the SEIR equations with a seasonally varying
# infection rate for this purpose.
#
# The equations are given as
#
# Sdot = mu - beta(t)*S*(I+i) - nu*S        (Susceptibles)
# Edot = beta(t)*S*(I+i) - (sigma+nu)*E     (Exposed)
# Idot = sigma*E - (gamma+nu)*I              (Infectious)
#
#
# Here beta(t) - the infection rate - is parameterized by a sinusoidal function
# plus a constant.
#
# Other parameters are
# i - a visiting process
# nu - death rate
# sigma - the rate of movement from Exposed to Infectious.
# gamma - the rate of recovery from infection.
#
# It is generally more stable to solve the equations for the log states rather
# than the states themselves. CollocInfer contains a number of useful tools
# that let you make this transition without needing to re-code all of your
# differential equations.

library('CollocInfer')

#### labels for the state process

SEIRnames = c('S','E','I')

#### define the right side evaluation function

SEIRodefn = make.SEIR()\$fn.ode

N = 5e6

####  define true parameter values and names

pars = c(7/5,0.02*N,0,0.02,365/8,365/5,17*(365/5)/N, 0, 0.08*17*(365/5)/N)
parnames = c('rho','mu','i','nu','sigma','gamma','b0','b1','b2')
names(pars) = parnames

#### A couple of functions to define the infection rate

beta.fun = function(t,p,more){
return( p['b0'] + p['b1']*sin(2*pi*t) + p['b2']*cos(2*pi*t) )
}

beta.dfdp = function(t,p,more){
dfdp =  cbind(rep(1,length(t)), sin(2*pi*t), cos(2*pi*t))
colnames(dfdp) = c('b0','b1','b2')
return(dfdp)
}

betamore = list(beta.fun=beta.fun,
beta.dfdp=beta.dfdp,
beta.ind=c('b0','b1','b2'))

#### Now create ODE paths

y0 = N*c(0.9,0,0.1)

y0 =N*c(.056,.00014,.000085)
names(y0) = SEIRnames
times = seq(0,3,1/52)

y = lsoda(y0,times=times,func=make.SEIR()\$ode.fn,parms=list(p=pars,more=betamore))

data = y[,4,drop=FALSE]*7/5

# and we're going to need some basis values

rr = range(times)
knots = seq(rr,rr,1/52)
norder = 3
nbasis = length(knots)+norder-2

bbasis = create.bspline.basis(range=rr,norder=norder,nbasis=nbasis,breaks=knots)
qpts = 0.5*(knots[1:(length(knots)-1)] + knots[2:length(knots)])

bvals.obs = Matrix(eval.basis(times,bbasis),sparse=TRUE)

bvals = list(bvals = Matrix(eval.basis(qpts,bbasis),sparse=TRUE),
dbvals = Matrix(eval.basis(qpts,bbasis,1),sparse=TRUE))

### Proc object, the first of these is just the standard squared error deviation
# from the right hand side of a differential equation

sproc = make.SSEproc()
sproc\$bvals = bvals
sproc\$more = make.SEIR()
sproc\$more\$more = betamore
sproc\$more\$qpts = qpts
sproc\$more\$weights = matrix(1,length(qpts),3)%*%diag(c(1e2,1e0,1e0))
sproc\$more\$names = SEIRnames
sproc\$more\$parnames = parnames

### However, ODEs are often much more numerically stable if represented on a
# log scale. The make.logtrans() function will take the right hand side functions
# and derivatives defined for any differential equation and provide the equivalent
# system for the log state. Note that this does affect the way  you represent
# your state when considering the observations.

lsproc = sproc
lsproc\$more = make.logtrans()
lsproc\$more\$more = make.SEIR()
lsproc\$more\$more\$more = betamore
lsproc\$more\$qpts = qpts
lsproc\$more\$weights = matrix(1,length(qpts),3)%*%diag(c(1e2,1e0,1e0))
lsproc\$more\$names = SEIRnames
lsproc\$more\$parnames = parnames

### Lik objects, this is the standard squared error.

slik = make.SSElik()
slik\$bvals = eval.basis(times,bbasis)
slik\$more = make.genlin()
slik\$more\$weights = rep(1,length(times))
slik\$more\$more = list(mat = matrix(0,1,3),sub=matrix(c(1,3,1),1,3))
slik\$more\$names = SEIRnames
slik\$more\$parnames = parnames

# Log transform transformation. For this, we note that we have represented
# the trajectory on the log scale and will need to transform back.

lslik = make.logstate.lik()
lslik\$bvals = slik\$bvals
lslik\$more\$weights = slik\$more\$weights
lslik\$more = slik
lslik\$more\$parnames = parnames

# Numerically things work much better on the log scale

dfd = data2fd(log(data/pars),times,bbasis)

coefs = matrix(0,nbasis,3)
coefs[,3] = dfd\$coefs

res = FitMatchOpt(coefs=coefs,which=1:2,proc=lsproc,pars=pars,meth='nlminb')

res2 = inneropt(data=data,times=times,pars=pars,proc=lsproc,lik=lslik,coefs=res\$coefs)

res3 = outeropt(data=data,times=times,pars=pars,proc=lsproc,lik=lslik,coefs=res2\$coefs,
active=c('rho','i','b0','b1','b2'))

## Or just log the observations, this is faster.

spars = pars
spars = 1

tdata = log(data/pars)

dfd = data2fd(tdata,times,bbasis)

coefs = matrix(0,nbasis,3)
coefs[,3] = dfd\$coefs

res = FitMatchOpt(coefs=coefs,which=1:2,proc=lsproc,pars=spars)

res2 = inneropt(data=tdata,times=times,pars=spars,proc=lsproc,lik=slik,coefs=res\$coefs)

res3 = outeropt(data=tdata,times=times,pars=spars,proc=lsproc,lik=slik,coefs=res2\$coefs,
active=c('i','b0','b1','b2'),out.meth='nlminb')
```

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