traj-match: Parameter estimation by fitting the trajectory of a model's...

In pomp: Statistical Inference for Partially Observed Markov Processes

Description

This function attempts to match trajectories of a model's deterministic skeleton to data. Trajectory matching is equivalent to maximum likelihood estimation under the assumption that process noise is entirely absent, i.e., that all stochasticity is measurement error. Accordingly, this method uses only the skeleton and dmeasure components of a POMP model.

Usage

  ## S4 method for signature 'pomp'
traj.match(object, start, est = character(0),
           method = c("Nelder-Mead","subplex","SANN","BFGS",
                      "sannbox","nloptr"),
           transform = FALSE, ...)
  ## S4 method for signature 'traj.matched.pomp'
traj.match(object, est, transform, ...)
  ## S4 method for signature 'pomp'
traj.match.objfun(object, params, est, transform = FALSE, ...)

Arguments

object	A pomp object. If object has no skeleton slot, an error will be generated.
start	named numeric vector containing an initial guess for parameters. By default start=coef(object) if the latter exists.
params	optional named numeric vector of parameters. This should contain all parameters needed by the skeleton and dmeasure slots of object. In particular, any parameters that are to be treated as fixed should be present here. Parameter values given in params for parameters named in est will be ignored. By default, params=coef(object) if the latter exists.
est	character vector containing the names of parameters to be estimated. In the case of traj.match.objfun, the objective function that is constructed will assume that its argument contains the parameters in this order.
method	Optimization method. Choices are subplex, “sannbox”, and any of the methods used by optim.
transform	logical; if TRUE, optimization is performed on the transformed scale.
...	Extra arguments that will be passed either to the optimizer (optim, subplex, nloptr, or sannbox, via their control (optim, subplex, sannbox) or opts (nloptr) lists) or to the ODE integrator. In traj.match, extra arguments will be passed to the optimizer. In traj.match.objfun, extra arguments are passed to trajectory. If extra arguments are needed by both optimizer and trajectory, construct an objective function first using traj.match.objfun, then give this objective function to the optimizer.

Details

In pomp, trajectory matching is the term used for maximizing the likelihood of the data under the assumption that there is no process noise. Specifically, traj.match calls an optimizer (optim, subplex, and sannbox are the currently supported options) to minimize an objective function. For any value of the model parameters, this objective function is calculated by

1. computing the deterministic trajectory of the model given the parameters. This is the trajectory returned by trajectory, which relies on the model's deterministic skeleton as specified in the construction of the pomp object object.

2. evaluating the negative log likelihood of the data under the measurement model given the deterministic trajectory and the model parameters. This is accomplished via the model's dmeasure slot. The negative log likelihood is the objective function's value.

The objective function itself — in a form suitable for use with optim-like optimizers — is created by a call to traj.match.objfun. Specifically, traj.match.objfun will return a function that takes a single numeric-vector argument that is assumed to cotain the parameters named in est, in that order.

Value

traj.match returns an object of class traj.matched.pomp. This class inherits from class pomp and contains the following additional slots:

transform, est

the values of these arguments on the call to traj.match.

evals

number of function and gradient evaluations by the optimizer. See optim.

value

value of the objective function. Larger values indicate better fit (i.e., traj.match attempts to maximize this quantity.

convergence, msg

convergence code and message from the optimizer. See optim.

Available methods for objects of this type include summary and logLik. The other slots of this object can be accessed via the $ operator.

traj.match.objfun returns a function suitable for use as an objective function in an optim-like optimizer.

See Also

trajectory, pomp, optim, subplex

Examples

  pompExample(ou2)
  true.p <- c(
	      alpha.1=0.9,alpha.2=0,alpha.3=-0.4,alpha.4=0.99,
	      sigma.1=2,sigma.2=0.1,sigma.3=2,
	      tau=1,
              x1.0=50,x2.0=-50
	      )
  simdata <- simulate(ou2,nsim=1,params=true.p,seed=43553)
  guess.p <- true.p
  res <- traj.match(
		    simdata,
		    start=guess.p,
		    est=c('alpha.1','alpha.3','alpha.4','x1.0','x2.0','tau'),
		    maxit=2000,
		    method="Nelder-Mead",
		    reltol=1e-8
		    )

  summary(res)

  plot(range(time(res)),range(c(obs(res),states(res))),type='n',xlab="time",ylab="x,y")
  points(y1~time,data=as(res,"data.frame"),col='blue')
  points(y2~time,data=as(res,"data.frame"),col='red')
  lines(x1~time,data=as(res,"data.frame"),col='blue')
  lines(x2~time,data=as(res,"data.frame"),col='red')

  pompExample(ricker)
  ofun <- traj.match.objfun(ricker,est=c("r","phi"),transform=TRUE)
  optim(fn=ofun,par=c(2,0),method="BFGS")

  pompExample(bbs)
  ## some options are passed to the ODE integrator
  ofun <- traj.match.objfun(bbs,est=c("beta","gamma"),transform=TRUE,hmax=0.001,rtol=1e-6)
  optim(fn=ofun,par=c(0,-1),method="Nelder-Mead",control=list(reltol=1e-10))

Example output

newly created object(s):
 ou2 
$params
    alpha.1     alpha.2     alpha.3     alpha.4     sigma.1     sigma.2 
  0.9152842   0.0000000  -0.3440927   0.9846510   2.0000000   0.1000000 
    sigma.3         tau        x1.0        x2.0 
  2.0000000   7.7347208  39.2445881 -63.8254689 

$loglik
[1] -691.9203

$eval
[1] 657  NA

$convergence
[1] 0

newly created object(s):
 ricker 
$par
[1] 2.062053 2.985564

$value
[1] 2135.11

$counts
function gradient 
      77       13 

$convergence
[1] 0

$message
NULL

In 'pomp': the following unrecognized argument(s) will be stored for use by user-defined functions: 'nbasis','degree','period','logvar','logitvar'
newly created object(s):
 bbs 
$par
[1]  0.7277326 -0.2751144

$value
[1] 94.01269

$counts
function gradient 
      95       NA 

$convergence
[1] 0

$message
NULL

pomp documentation built on May 2, 2019, 4:09 p.m.