# inst/Devel.Examples/pomp.sir.R In CollocInfer: Collocation Inference for Dynamic Systems

```require(pomp)
library(CollocInfer)

## euler.sim pomp object from the pomp library

data(euler.sir)

## Parameters

pars = coef(euler.sir)
parnames = names(pars)

## Observations and observation times; we'll start by assuming that the
## observations are given by a scaled version of the infectives.

Y = t(data.array(euler.sir))
times = time(euler.sir)

Data = matrix(NA,nrow(Y),5)   ## Easiest way to specify that we only observe
Data[,2] = log(Y) + 1.217     ## the infectives.

## Solve the ODE to get a trajectory and state names

x <- trajectory(euler.sir)
varnames =  dimnames(x)[]  ## I really wish I could extract this sort of thing
x = t(x[,1,])                 ## directly from the pomp object

x[,5] = 1

## Now set up a basis representation for the ODE and start by smoothing the
## solved trajectory

ebasis = create.bspline.basis(c(0,4),norder=4,breaks=seq(0,4,0.1))

trajfd = smooth.basis(times,log(x[-1,]),fdPar(ebasis,int2Lfd(2),0.01))

coefs = trajfd\$fd\$coefs
colnames(coefs) = varnames

## Now we'll set up the CollocInfer objects

objs = LS.setup(pars=pars, coefs=coefs, fn=euler.sir,
basisvals=ebasis, lambda = 1e1,
data = Data, times=times, posproc=TRUE)

## And try smoothing

res1 = inneropt(Data,times,pars,coefs,objs\$lik,objs\$proc,control.in=list(trace=6))

## And the outer optimization

res2 = outeropt(Data,times,pars,res1\$coefs,objs\$lik,objs\$proc,active=1:3)

##########################
## Try without logging
#############################

Data = exp(Data)

trajfd = smooth.basis(times,x[-1,],fdPar(ebasis,int2Lfd(2),0.00001))

coefs = trajfd\$fd\$coefs
colnames(coefs) = varnames

coefs[,2] = smooth.basis(times,Data[,2],fdPar(ebasis,int2Lfd(2),0.00001))\$fd\$coef

## Now we'll set up the CollocInfer objects

objs = LS.setup(pars=pars, coefs=coefs, fn=euler.sir,
basisvals=ebasis, lambda = 1e-3,
data = Data, times=times, posproc=FALSE)

## And try smoothing

res1 = inneropt(Data,times,pars,coefs,objs\$lik,objs\$proc,control.in=list(trace=6))

## And the outer optimization

res2 = outeropt(Data,times,pars,res1\$coefs,objs\$lik,objs\$proc,active=1:3)

##### Now something more curious on the observation process; we'll use the difference
##### in the cases state variable.

proc = objs\$proc
lik = objs\$lik

lik\$bvals = eval.basis(times,ebasis) - eval.basis(c(0,times[-208]),ebasis)

Data = matrix(0,208,5)
Data[,4] = Y

## And try smoothing

res1 = inneropt(Data,times,pars,coefs,lik,proc,control.in=list(trace=6))

## And the outer optimization

res2 = outeropt(Data,times,pars,res1\$coefs,lik,proc,active=1:3)
```

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