Collocation Inference in R

Description

Functions carry out collocation inference method for nonlinear continuous-time dynamic systems. These are based on basis-expansion representations for the state of the system. Gradient-matching, profiling and EM algorithms are supported.

Details

Author(s)

Giles Hooker, Luo Xiao

Maintainer: Giles Hooker <giles.hooker@cornell.edu>

References

Ramsay, James O., Giles Hooker, Jiguo Cao and David Campbell (2007), "Parameter Estimation in Ordinary Differential Equations: A Generalized Smoothing Approach", Journal of the Royal Statistical Society, 69

Ramsay, James O., and Silverman, Bernard W. (2006), Functional Data Analysis, 2nd ed., Springer, New York.

Chemostat Example Data

Description

Five-species Chemostat Model

Usage

ChemoData

Format

ChemoData

A 61 by 2 matrix of data observed in a chemostat.

ChemoTime

A vector of 61 observation times corresponding to ChemoData.

ChemoPars

Named parameter vector as a starting point for estimation ChemoData.

ChemoVarnames

c('N','C1','C2','B','S'): the state variable names for the chemostat system.

ChemoParnames

parameter names for the chemostat system.

Source

Yoshida, T., L. E. Jones, S. P. Ellner, G. F. Fussmann and N. G. Hairston, 2003, "Rapid evolution drives ecological dynamics in a predator-prey system", Nature, 424, pp. 303-306.

Rosenzweig-MacArthur Model Applied to Chemostat Data

Description

Two-Species Rosenzweig-MacArthur Model

Usage

ChemoRMData

Format

ChemoRMData

A 108 by 2 matrix of data observed in a chemostat.

ChemoRMPars

Named parameter vector as a starting point for estimation in ChemoRMData.

ChemoRMTime

A vector of 108 observation times corresponding to ChemoData.

RMparnames

parameter names for the Rosenzweig-MacArthur system.

RMvarnames

the state variable names for the Rosenzweig-MacArthur system.

Source

Becks, L., S. P. Ellner, L. E. Jones, and N. G. Hairston, 2010, "Reduction of adaptive genetic diversity radically alters eco-evolutionary community dynamics", Ecology Letters, 13, pp. 989-997.

CollocInfer internal functions

Description

Internal undocumentation functions

Usage

blocks2mat(H)

Diagnostic PLots for CollocInfer

Description

Diagnostic Plots on the Results of CollocInfer

Usage

CollocInferPlots(coefs,pars,lik,proc,times=NULL,data=NULL,
      cols=NULL,datacols=NULL,datanames=NULL,ObsPlot=TRUE,DerivPlot=TRUE,
      cex.axis=1.5,cex.lab=1.5,cex=1.5,lwd=2)

Arguments

Details

Timevec is taken to be the quadrature values. Three plots can be produced:

If ObsPlot=TRUE a plot is given of the predicted values of the observations along with the observations themselves (if given).

If DerivPlot=TRUE two plots are produced. The first gives the value of the derivative of the estimated trajectory (dashed) and the value of the right-hand-side of the ordinary differential equation in proc (hence the predicted derivative) (solid). The second plot gives their difference in the first panel as well as the estimated trajectory in the second panel.

Value

A list containing elements used in plotting:

FitzHugh-Nagumo data

Description

Data generated for FitzHugh-Nagumo Examples

Usage

FhNdata

Format

FhNdata

A 41 by 2 matrix of data generated from the FitzHugh Nagumo equations.

FhNtimes

A vector of 41 observation times corresponding to FhNdata.

FhNpars

Named parameter vector used to generate FhNdata.

FhNvarnames

c('V','R'): the state variable names for the FitzHugh Nagumo system.

FhNparnames

c('a','b','c') parameter names for the FitzHugh Nagumo system.

Source

James Ramsay, Giles Hooker David Campbell and Jiguo Cao, 2007. "Parameter Estimation for Differential Equations: A Generalized Smoothing Approach". Journal of the Royal Statistical Society Vol 69 No 5.

Estimated Parameters for FitzHugh-Nagumo data

Description

Parameters Estimated for FhN Data – used to speed up examples

Usage

FhNestPars

Format

FhNestPars

Estimated parameters for the FhN Data example.

FhNestCoefs

Estimated coefficients for the FhN Data example.

Source

James Ramsay, Giles Hooker David Campbell and Jiguo Cao, 2007. "Parameter Estimation for Differential Equations: A Generalized Smoothing Approach". Journal of the Royal Statistical Society Vol 69 No 5.

Estimating Hidden States

Description

Estimating hidden states to maximize agreement with the process.

Usage

FitMatchOpt(coefs,which,pars,proc,meth='nlminb',control=list())

FitMatchErr(coefs,allcoefs,which,pars,proc,sgn=1)

FitMatchDC(coefs,allcoefs,which,pars,proc,sgn=1)

FitMatchDC2(coefs,allcoefs,which,pars,proc,sgn=1)

FitMatchList(coefs,allcoefs,which,pars,proc,sgn=1)

Arguments

Details

These routines allow the values of coefficients for some states to be optimized relative to the others. That is, the objective defined by proc is minimized over those states specified in which leaving the others constant. This would be typically done, for example, a smooth is taken to estimate some states non-parametrically, but data is not available on all of them.

A number of optimization routines have been implemented in FitMatchOpt, some experimentation is advised.

Value

See Also

ParsMatchErr, SplineCoefsErr, inneropt

Examples

###############################
#### Some Data            #####
###############################

data(FhNdata)

# And parameter estimates

data(FhNest)


###############################
####  Basis Object      #######
###############################

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

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


# Initial values for coefficients will be obtained by smoothing

fd.data = FhNdata[,1]

DEfd = smooth.basis(FhNtimes,fd.data,fdPar(bbasis,1,0.5))   

coefs = cbind(DEfd$fd$coefs,rep(0,nbasis))
colnames(coefs) = FhNvarnames

#############################################################
### If We Only Observe One State, We Can Re-Smooth Others ### 
#############################################################

profile.obj = LS.setup(pars=FhNpars,coefs=coefs,fn=make.fhn(),
                      basisvals=bbasis,lambda=1000,times=FhNtimes)
lik = profile.obj$lik
proc= profile.obj$proc

#  DD = Matrix(diag(1,200),sparse=TRUE)
#  tDD = t(DD)

fres = FitMatchOpt(coefs=coefs,which=2,pars=FhNpars,proc)

plot(fd(fres$coefs,bbasis))

IntegrateForward

Description

Solves a differential equation going forward based on a proc object.

Usage

IntegrateForward(y0,ts,pars,proc,more)

Arguments

Value

Returns the output from solving the differential equation using the lsoda routines. Specifically, it returns a list with elements

times

The output times.

states

The output states.

See Also

Profile.LS, Profile.multinorm

Examples

proc = make.SSEproc()
proc$more = make.fhn()
proc$more$names = c('V','R')

y0 = c(-1,1)
names(y0) = c('V','R')

pars = c(0.2,0.2,3)
names(pars) = c('a','b','c')

ts = seq(0,20,0.5)

value = IntegrateForward(y0,ts,pars,proc)

matplot(value$times,value$states)

North Shore data

Description

Groundwater Data from Vancouver's North Shore

Usage

NSgroundwater

Format

NSgroundwater

A 315 by 1 matrix of data on groundwater level collected in vancouver.

NStimes

A vector of 315 observation times corresponding to NSgroundwater.

NSrainfall

Rainfall as a covariate to NSgroundwater; this quantity is lagged by 3 days.

Estimate of Parameters from Smooth

Description

Objective function and derivatives to estimate parameters with a fixed smooth.

Usage

ParsMatchOpt(pars,coefs,proc,active=1:length(pars),meth='nlminb',control=list())

ParsMatchErr(pars,coefs,proc,active=1:length(pars),allpars,sgn=1)

ParsMatchDP(pars,coefs,proc,active=1:length(pars),allpars,sgn=1)

ParsMatchList(pars,coefs,proc,active=1:length(pars),allpars,sgn=1)

Arguments

Details

These routines fix the estimated states at the value given by coefs and estimate pars to maximize agreement between the fixed state and the objective given by the proc object.

A number of optimization routines have been implemented in FitMatchOpt, some experimentation is advised.

Value

See Also

FitMatchErr, SplineCoefsErr, inneropt

Examples

data(FhNdata)

###############################
####  Basis Object      #######
###############################

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

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


# Initial values for coefficients will be obtained by smoothing

DEfd = smooth.basis(FhNtimes,FhNdata,fdPar(bbasis,1,0.5))   # Smooth to estimate
                                                            # coefficients first
coefs = DEfd$fd$coefs
colnames(coefs) = FhNvarnames



#################################
### Initial Parameter Guesses ###
#################################

profile.obj = LS.setup(pars=FhNpars,coefs=coefs,fn=make.fhn(),basisvals=bbasis,
    lambda=1000,times=FhNtimes)
lik = profile.obj$lik
proc= profile.obj$proc

pres = ParsMatchOpt(FhNpars,coefs,proc)

npars = pres$pars

Profile.covariance

Description

Newey-West estimate of covariance of parameter estimates from profiling.

Usage

Profile.covariance(pars,active=NULL,times,data,coefs,lik,proc,
          in.meth='nlminb',control.in=NULL,eps=1e-6,GN=FALSE)

Arguments

Details

Currently assumes a lag-5 auto-correlation among observation vectors.

Value

Returns a Newey-West estimate of the covariance matrix of the parameter estimates.

See Also

ProfileErr, ProfileSSE, Profile.LS, Profile.multinorm

Examples

# See example in Profile.LS

Profile Estimation with Collocation Inference

Description

Profile estimation and objective functions for collocation estimation of parameters in continuous-time stochastic processes.

Usage

Profile.GausNewt(pars,times,data,coefs,lik,proc,in.meth="nlminb",
      control.in=NULL,active=1:length(pars),
      control=list(reltol=1e-6,maxit=50,maxtry=15,trace=1))

ProfileSSE(pars,allpars,times,data,coefs,lik,proc,in.meth='nlminb',
      control.in=NULL,active=1:length(pars),dcdp=NULL,oldpars=NULL,
      use.nls=TRUE,sgn=1)
      
ProfileErr(pars,allpars,times,data,coefs,lik,proc,in.meth = "house",
      control.in=NULL,sgn=1,active=1:length(allpars))

ProfileDP(pars,allpars,times,data,coefs,lik,proc,in.meth = "house",
      control.in=NULL,sgn=1,sumlik=1,active=1:length(allpars))

ProfileList(pars,allpars,times,data,coefs,lik,proc,in.meth = "house",
      control.in=NULL,sgn=1,active=1:length(allpars))

Arguments

Details

Profile.GausNewt provides a simple implementation of Gauss-Newton optimization and may not result in optimized values that are as good as other optimizers in R.

When Profile.GausNewt is not being used, ProfileSEE and ProfileERR create the following temporary files:

counter.tmp

The number of halving-steps taken on the current optimization step.

curcoefs.tmp

The current estimates of the coefficients.

optcoefs.tmp

The optimal estimates of the coefficients in the iteration history.

These need to be removed manually when the optimization is finished. optcoefs.tmp will contain the optimal value of coefs for plotting the estimated trajectories.

Value

See Also

outeropt, Profile.LS, Profile.multinorm, LS.setup, multinorm.setup

Profile Estimation Functions

Description

These functions are wrappers that create lik and proc functions and run generalized profiling.

Usage

Profile.LS(fn,data,times,pars,coefs=NULL,basisvals=NULL,lambda,
                        fd.obj=NULL,more=NULL,weights=NULL,quadrature=NULL,
                        likfn = make.id(), likmore = NULL,
                        in.meth='nlminb',out.meth='nls',
                        control.in,control.out,eps=1e-6,active=NULL,posproc=FALSE,
                        poslik=FALSE,discrete=FALSE,names=NULL,sparse=FALSE)

Profile.multinorm(fn,data,times,pars,coefs=NULL,basisvals=NULL,var=c(1,0.01),
                        fd.obj=NULL,more=NULL,quadrature=NULL,
                        in.meth='nlminb',out.meth='optim',
                        control.in,control.out,eps=1e-6,active=NULL,
                        posproc=FALSE,poslik=FALSE,discrete=FALSE,names=NULL,sparse=FALSE)

Arguments

Details

These functional all carry out the profiled optimization method of Ramsay et al 2007. Profile.LS uses a sum of squared errors criteria for both fit to data and the fit of the derivatives to a differential equation. Profile.multinorm uses multivariate normal approximations. discrete changes the system to a discrete-time difference equation with the right hand side function giving the transition function.

Note that these all call outeropt, which creates temporary files 'curcoefs.tmp' and 'optcoefs.tmp' to update coefficients as pars evolves. These overwrite existing files of those names and are deleted before the function terminates.

Value

A list with elements

See Also

outeropt, ProfileErr, ProfileSSE, LS.setup, multinorm.setup

Examples

###############################
####   Data             #######
###############################

data(FhNdata)

###############################
####  Basis Object      #######
###############################

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

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


#### Start from pre-estimated values to speed up optimization

data(FhNest)

spars = FhNestPars
coefs = FhNestCoefs

lambda = 10000

res1 = Profile.LS(make.fhn(),data=FhNdata,times=FhNtimes,pars=spars,coefs=coefs,
  basisvals=bbasis,lambda=lambda,in.meth='nlminb',out.meth='nls')

Covar = Profile.covariance(pars=res1$pars,times=FhNtimes,data=FhNdata,
  coefs=res1$coefs,lik=res1$lik,proc=res1$proc) 


## Not run: 
## Alternative, starting from perturbed coefficients -- takes too long for 
# automatic checks in CRAN

# Initial values for coefficients will be obtained by smoothing

DEfd = smooth.basis(FhNtimes,FhNdata,fdPar(bbasis,1,0.5))   # Smooth to estimate
                                                            # coefficients first
coefs = DEfd$fd$coefs
colnames(coefs) = FhNvarnames


###############################
####     Optimization       ###
###############################

spars = c(0.25,0.15,2.5)          # Perturbed parameters
names(spars)=FhNparnames
lambda = 10000

res1 = Profile.LS(make.fhn(),data=FhNdata,times=FhNtimes,pars=spars,coefs=coefs,
  basisvals=bbasis,lambda=lambda,in.meth='nlminb',out.meth='nls')

par(mfrow=c(2,1))
plotfit.fd(FhNdata,FhNtimes,fd(res1$coefs,bbasis))

## End(Not run)  
  
 
 
  
## Not run: 
####################################################
###  An Explicitly Multivariate Normal Formation ### 
####################################################

var = c(1,0.0001)

res2 = Profile.multinorm(make.fhn(),FhNdata,FhNtimes,pars=res1$pars,
          res1$coefs,bbasis,var=var,out.meth='nlminb', in.meth='nlminb')

## End(Not run)

SEIR data

Description

Data generated for SEIR Examples

Usage

SEIRdata

Format

SEIRdata

A 261 by 1 matrix of data generated from the SEIR Gillespie process run over 5 years.

SEIRtimes

A vector of 261 observation times corresponding to SEIRdata.

SEIRpars

Named parameter vector used to generate SEIRdata.

SEIRvarnames

c('V','R'): the state variable names for the SEIR system.

SEIRparnames

parameter names for the SEIR system.

Source

Giles Hooker, Stephen P. Ellner, David Earn and Laura Roditi, 2010. "Parameterizing State-space Models for Infectious Disease Dynamics by Generalized Profiling: Measles in Ontario", Technical Report, Cornell University.

Model-Based Smoothing Functions

Description

Perform the inner optimization to estimate coefficients given parameters.

Usage

Smooth.LS(fn,data,times,pars,coefs=NULL,basisvals=NULL,lambda,fd.obj=NULL,
        more=NULL,weights=NULL,quadrature=NULL,likfn = make.id(), 
        likmore = NULL,in.meth='nlminb',control.in,eps=1e-6,
        posproc=FALSE,poslik=FALSE,discrete=FALSE,names=NULL,
        sparse=FALSE)
        
Smooth.multinorm(fn,data,times,pars,coefs=NULL,basisvals=NULL,var=c(1,0.01),
        fd.obj=NULL,more=NULL,quadrature=NULL,in.meth='nlminb',
        control.in,eps=1e-6,posproc=FALSE,poslik=FALSE,discrete=FALSE,
        names=NULL,sparse=FALSE)

Arguments

Details

These routines create lik and proc objects and call inneropt.

Value

A list with elements

See Also

inneropt, LS.setup, multinorm.setup, SplineCoefsErr

Examples

###############################
####   Data             #######
###############################

data(FhNdata)

###############################
####  Basis Object      #######
###############################

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

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


#### Start from pre-estimated values to speed up optimization

data(FhNest)

spars = FhNestPars
coefs = FhNestCoefs

lambda = 10000

res1 = Smooth.LS(make.fhn(),data=FhNdata,times=FhNtimes,pars=spars,coefs=coefs,
  basisvals=bbasis,lambda=lambda,in.meth='nlminb')


## Not run: 
# Henon system

hpars = c(1.4,0.3)              # Parameters
t = 1:200

x = c(-1,1)                     # Create some dataa
X = matrix(0,200+20,2)
X[1,] = x

for(i in 2:(200+20)){ X[i,] = make.Henon()$ode(i,X[i-1,],hpars,NULL) }

X = X[20+1:200,]

Y = X + 0.05*matrix(rnorm(200*2),200,2)

basisvals = diag(rep(1,200))    # Basis is just identiy
coefs = matrix(0,200,2)


# For sum of squared errors

lambda = 10000

res1	= Smooth.LS(make.Henon(),data=Y,times=t,pars=hpars,coefs,basisvals=basisvals,
  lambda=lambda,in.meth='nlminb',discrete=TRUE)

## End(Not run)


## Not run: 
# For multinormal transitions

var = c(1,0.01)

res2 = Smooth.multinorm(make.Henon(),data=Y,t,pars=hpars,coefs,basisvals=NULL,
  var=var,in.meth='nlminb',discrete=TRUE)

## End(Not run)

Spline Estimation Functions

Description

Model-based smoothing; estimation, objective criterion and derivatives.

Usage

SplineEst.NewtRaph(coefs,times,data,lik,proc,pars,
      control=list(reltol=1e-12,maxit=1000,maxtry=10,trace=0))
      
SplineCoefsList(coefs,times,data,lik,proc,pars,sgn=1)

SplineCoefsErr(coefs,times,data,lik,proc,pars,sgn=1)

SplineCoefsDC(coefs,times,data,lik,proc,pars,sgn=1)

SplineCoefsDP(coefs,times,data,lik,proc,pars,sgn=1)

SplineCoefsDC2(coefs,times,data,lik,proc,pars,sgn=1)

SplineCoefsDCDP(coefs,times,data,lik,proc,pars,sgn=1)

Arguments

Details

SplineEst.NewtRaph performs a simple Newton-Raphson estimate for the optimal value of the coefficients. This estimate lacks the convergence checks of other estimation packages, but may yeild a fast solution when needed.

Value

The output of gradients is in terms of an array with dimensions corresponding to derivatives. Derivatives with with respect to coefficients are given in dimensions before those that give derivatives with respect to parameters.

See Also

inneropt, Smooth.LS

forward.prediction.error

Description

Forward prediction error objective for choice of lambda in square error criteria.

Usage

forward.prediction.error(times,data,coefs,lik,proc,pars,whichtimes=NULL)

Arguments

Details

Forward prediction error can be used to choose values of lambda in the profiled estimation routines. The ordinary differential equation is solved starting from the starting times specified in whichtimes and measured at the corresponding measurement times. The error is then recorded. This should then be minimized by a grid search.

Value

The forwards prediction error from the estimates.

See Also

ProfileSSE, outeropt

Inner Optimization Functions

Description

Estmates coefficients given parameters.

Usage

inneropt(data,times,pars,coefs,lik,proc,in.meth='nlminb',control.in=list())

Arguments

Details

This minimizes the objective function defined by the addition of the lik and proc objectives with respect to the coefficients. A number of generic optimization routines can be used and some experimentation is recommended.

Value

A list with elements

See Also

outeropt, Smooth.LS,LS.setup, multinorm.setup, SplineCoefsErr

Examples

## Not run: 
# FitzHugh-Nagumo Equations

data(FhNdata)   # Some data
data(FhNest)    # with some parameter estimates

knots = seq(0,20,0.2)         # Create a basis
norder = 3
nbasis = length(knots) + norder - 2
range = c(0,20)

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

lambda = 10000               # Penalty value

DEfd = smooth.basis(FhNtimes,FhNdata,fdPar(bbasis,1,0.5))   # Smooth to estimate
                                                            # coefficients first
coefs = DEfd$fd$coefs
colnames(coefs) = FhNvarnames

profile.obj = LS.setup(pars=FhNpars,coefs=coefs,fn=make.fhn(),
                        basisvals=bbasis,lambda=lambda,times=FhNtimes)

lik = profile.obj$lik
proc= profile.obj$proc

res = inneropt(FhNdata,times=FhNtimes,FhNpars,coefs,lik,proc,in.meth='nlminb')

plot(fd(res$coefs,bbasis))

## End(Not run)

Finite Difference Functions

Description

Returns a list of functions that calculate finite difference derivatives.

Usage

make.findif.ode()

make.findif.loglik()

make.findif.var()

Details

All these functions require the sepcification of more$eps to give the size of the finite differencing step. They also require more to specify the original object (ODE right hand side functions, definitions of lik and proc objects).

Value

A list of functions that calculate the derivatives via finite differencing schemes.

See Also

LS.setup, multinorm.setup

Examples

# Sum of squared errors with finite differencing to get right-hand-side derivatives

proc = make.SSEproc()
proc$more = make.findif.ode()


# Finite differencing for the log likelihood

lik = make.findif.loglik()
lik$more = make.SSElik()


# Multivariate normal transitions with finite differencing for mean and variance functions

lik = make.multinorm()
lik$more = c(make.findif.ode,make.findif.var)

# Finite differencing for transition density of a discrete time system

proc = make.Dproc()
proc$more = make.findif.loglik()

Observation Process Distribution Function

Description

Returns a list of functions that calculate the observation process distribution and its derivatives; designed to be used with the collocation inference functions.

Usage

make.SSElik()

make.multinorm()

Details

These functions require more to be a list with elements:

fn

The transform function of the states to observations, or to their derivatives.

dfdx

The derivative of fn with respect to states.

dfdp

The derivative of fn with respect to parameters.

d2fdx2

The second derivative of fn with respect to states.

d2fdxdp

The cross derivative of fn with respect to states and parameters.

make.Multinorm further requires:

var.fn

The variance given in terms of states and parameters.

var.dfdx

The derivative of var.fn with respect to states.

var.dfdp

The derivative of var.fn with respect to parameters.

var.d2fdx2

The second derivative of var.fn with respect to states.

var.d2fdxdp

The cross derivative of var.fn with respect to states and parameters.

make.SSElik further requres weights giving weights to each observation.

Value

A list of functions that calculate the log observation distribution and its derivatives.

See Also

LS.setup, multinorm.setup

Examples

# Straightforward sum of squares:

lik = make.SSElik()
lik$more = make.id()

# Multivariate normal about an exponentiated state with constant variance

lik = make.multinorm()
lik$more = c(make.exp(),make.cvar())

Log Transforms

Description

Functions to modify liklihood, transform, lik and proc objects so that the operate with the state defined on a log scale.

Usage

make.logtrans()

make.exptrans()

make.logstate.lik()

make.exp.Cproc()

make.exp.Dproc()

Details

All functions require more to specify the original object (ODE right hand side functions, definitions of lik and proc objects).

Value

A list of functions that calculate log transforms and derivatives in various contexts.

See Also

LS.setup, make.Cproc, make.Dproc

Examples

# Model the log of an SEIR process

proc = make.SSEproc()
proc$more = make.logtrans()
proc$more$more = make.SEIR()

# Observe a linear combination  of

lik = make.logstate.lik()
lik$more = make.SSElik()
lik$more$more = make.genlin()

# SEIR Model with multivariate transition densities

proc = make.exp.Cproc()
proc$more = make.multinorm()
proc$more$more = c(make.SEIR(),make.cvar())

Process Distributions

Description

Functions to define process distributions in the collocation inference package.

Usage

make.Dproc()

make.Cproc()

make.SSEproc()

Details

All functions require more to specify this distribution. This should be a list containing

fn

The distribution specified.

dfdx

The derivative of fn with respect to states.

dfdp

The derivative of fn with respect to parameters.

d2fdx2

The second derivative of fn with respect to states.

d2fdxdp

The cross derivative of fn with respect to states and parameters.

For Cproc and Dproc this should specify the distribution; for SSEproc it should specify the right hand side of a differential equation.

Value

A list of functions that the process distribution

See Also

LS.setup, multinorm.setup

Examples

# FitzHugh-Nagumo Equations

proc = make.SSEproc()
proc$more = make.fhn()

# Henon Map

proc = make.Dproc()
proc$more = make.Henon


# SEIR with multivariate normal transitions

proc = make.Cproc()
proc$more = make.multinorm()
proc$more$more = c(make.SEIR(),make.var.SEIR())

Transfer Functions

Description

Returns a list of functions that calculate the transform and its derivatives.

Usage

make.id()

make.exp()

make.genlin()

make.fhn()

make.Henon()

make.SEIR()

make.NS()

chemo.fun(times,y,p,more=NULL)

Arguments

All the functions created by make... functions, require the arguments needed by chemo.fun

Details

make.genlin requires the specification of further elements in the list. In particular the element more should be a list containing

mat

a matrix defining the linear transform before any parameters are added. This may be all zero, but it may also specify fixed elements, if desired.

sub

a k-by-3 matrix indicating which parameters should be entered into which elements of mat. Each row is a triple giving the row and colum of mat to be specified and the element of the parameter vector that should be substituted. sub over-rides any values in mat.

force

if input functions are given, these are given as a list.

force.mat

specifying the influence of the elements of force on the state variables. Defined as in mat.

force.sub

defined as in sub, over-rides the elements of force.mat with parameter values.

make.diagnostics estimates forcing-function diagnostics as in Hooker, 2009 for goodness-of-fit assessment. It requires

psi

Values of a basis expansion for forcing functions at the quadrature points.

which

Which states are to be forced?

fn, dfdx, d2fdx2

Functions and derivatives as would be used to estimate parameters for the original equations.

pars

Parameters to go into more$fn.

make.SEIR estimates parameters and a seasonal variation in the infection rate in an SEIR model. It requires the specification of the seasonal change rate in more by a list with objects

beta.fun

A function to calculate beta, it should have arguments t, p and betadef and return a matrix giving the value of beta at times t with parameters p.

beta.dfdp

Should calculate the derivative of beta.fun with respect to p, at times t returning a matrix. The matrix should be of size length(t) by length(p) where p is the entire parameter vector.

betadef

Additional inputs (eg bases) to beta.fun and beta.dfdp.

make.NS provides functions for the North Shore example. This is a possibly time-varying forced linear system of one dimension. It requires more to specify betabasis to describe the autoregressive coefficient, and alphabasis to provide a contant in front of the functional data object rainfd.

chemo.fun Is a five-state predator-prey-resources model used as an example. It stands alone as a function and should be used with the findif.ode functions.

Value

A list of functions that calculate the transform and its derivatives, in a form compatible with the collocation inference functions.

See Also

LS.setup, multinorm.setup

Examples

# Observe the FitzHugh-Nagumo equations

proc = make.SSEproc()
proc$more = make.fhn()

lik = make.SSElik()
lik$more = make.id()

# Observe an unknown scalar transform of each component of a Henon map, given
# in the first two elements of the parameter vector:

proc = make.Dproc()
proc$more = make.multinorm()
proc$more$more = c(make.Henon,make.cvar)

lik = make.multinorm()
lik$more = c(make.genlin,make.cvar)
lik$more$more = list(mat = matrix(0,2,2),sub=matrix(c(1,1,1,2,2,2),2,3,byrow=TRUE))

# Model SEIR equations on the log scale and then exponentiate

lik = make.SSElik()
lik$more = make.exp()

proc = make.SSEproc()
proc$more = make.logtrans()
proc$more$more = make.SEIR()

Variance Functions

Description

Returns a list of functions that calculate a (possibly state and parameter dependent) variance.

Usage

make.cvar()

make.var.SEIR()

Details

make.cvar requires the specification of further elements in the list. In particular the element more should be a list containing

Value

A list of functions that calculate a variance function and its derivatives, in a form compatible with the collocation inference functions.

See Also

make.multinorm

Examples

# Multivariate normal observation of the state vector.

lik = make.multinorm()
lik$more = c(make.id(),make.cvar())

Outer Optimization Functions

Description

Outer optimization; performs profiled estimation.

Usage

outeropt(data,times,pars,coefs,lik,proc,
      in.meth='nlminb',out.meth='nlminb',
      control.in=list(),control.out=list(),active=1:length(pars))

Arguments

Details

The outer optimization for parameters looks only at the objective defined by the lik object. For every parameter value, coefs are optimized by inneropt and then the value of lik for these coefficients is computed.

A number of optimization routines can be used here, some experimentation is recommended. Libraries for these optimization routines are not pre-loaded. Where these functions take options as explicit arguments instead of a list, they should be listed in control.out and will be called by their names.

The routine creates temporary files 'curcoefs.tmp' and 'optcoefs.tmp' to update coefficients as pars evolves. These overwrite existing files of those names and are deleted before the function terminates.

Value

A list containing

See Also

inneropt, Profile.LS, ProfileSSE, ProfileErr, LS.setup, multinorm.setup

Examples

## Not run: 
data(FhNdata)
data(FhNest)

knots = seq(0,20,0.2)         # Create a basis
norder = 3
nbasis = length(knots) + norder - 2
range = c(0,20)

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

lambda = 10000               # Penalty value

DEfd = smooth.basis(FhNtimes,FhNdata,fdPar(bbasis,1,0.5))   # Smooth to estimate
                                                            # coefficients first
coefs = DEfd$fd$coefs
colnames(coefs) = FhNvarnames

profile.obj = LS.setup(pars=FhNpars,coefs=coefs,fn=make.fhn(),basisvals=bbasis,
      lambda=lambda,times=FhNtimes)

lik = profile.obj$lik
proc= profile.obj$proc

res = outeropt(data=FhNdata,times=FhNtimes,pars=FhNpars,coefs=coefs,lik=lik,proc=proc,
    in.meth="nlminb",out.meth="nlminb",control.in=NULL,control.out=NULL)


plot(res$coefs,main='outeropt')
print(blah)

## End(Not run)

Setup Functions for proc and lik objects

Description

These functions set up lik and proc objects of squared error and multinormal processes.

Usage

LS.setup(pars,coefs=NULL,fn,basisvals=NULL,lambda,fd.obj=NULL,
        more=NULL,data=NULL,weights=NULL,times=NULL,quadrature=NULL,
        likfn = make.id(), likmore = NULL,eps=1e-6,
        posproc=FALSE,poslik=FALSE,discrete=FALSE,names=NULL,sparse=FALSE)

multinorm.setup(pars,coefs=NULL,fn,basisvals=NULL,var=c(1,0.01),fd.obj=NULL,
        more=NULL,data=NULL,times=NULL,quadrature=NULL,eps=1e-6,posproc=FALSE,
        poslik=FALSE,discrete=FALSE,names=NULL,sparse=FALSE)

Arguments

Details

These functions provide basic setup utilities for the collocation inference methods. They define lik and proc objects for sum of squared errors and multivariate normal log likelihoods with nonlinear transfer functions describing the evolution of the state vector.

LS.setup

Creates sum of squares functions

multinorm.setup

Creates multinormal log likelihoods for a continuous-time system.

Value

A list with elements

See Also

inneropt, outeropt, Profile.LS, Profile.multinorm, Smooth.LS, Smooth.multinorm

Examples

# FitzHugh-Nagumo

t = seq(0,20,0.05)            # Observation times

pars = c(0.2,0.2,3)           # Parameter vector
names(pars) = c('a','b','c')

knots = seq(0,20,0.2)         # Create a basis
norder = 3
nbasis = length(knots) + norder - 2
range = c(0,20)

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

lambda = 10000               # Penalty value

coefs = matrix(0,nbasis,2)   # Coefficient matrix

profile.obj = LS.setup(pars=pars,coefs=coefs,fn=make.fhn(),basisvals=bbasis,
                       lambda=lambda,times=t)


# Using multinorm

var = c(1,0.01)

profile.obj = multinorm.setup(pars=pars,coefs=coefs,fn=make.fhn(),
                                        basisvals=bbasis,var=var,times=t)


# Henon - discrete

hpars = c(1.4,0.3)
t = 1:200

coefs = matrix(0,200,2)
lambda = 10000

profile.obj = LS.setup(pars=hpars,coefs=coefs,fn=make.Henon(),basisvals=NULL,
                             lambda=lambda,times=t,discrete=TRUE)