trajectory {pomp}

Trajectory of a deterministic model

Description

Compute trajectories of the deterministic skeleton of a Markov process.

Usage

## S4 method for signature 'missing'
trajectory(
  t0,
  times,
  params,
  skeleton,
  rinit,
  ...,
  ode_control = list(),
  format = c("pomps", "array", "data.frame"),
  verbose = getOption("verbose", FALSE)
)

## S4 method for signature 'data.frame'
trajectory(
  object,
  ...,
  t0,
  times,
  params,
  skeleton,
  rinit,
  ode_control = list(),
  format = c("pomps", "array", "data.frame"),
  verbose = getOption("verbose", FALSE)
)

## S4 method for signature 'pomp'
trajectory(
  object,
  params,
  ...,
  skeleton,
  rinit,
  ode_control = list(),
  format = c("pomps", "array", "data.frame"),
  verbose = getOption("verbose", FALSE)
)

## S4 method for signature 'traj_match_objfun'
trajectory(object, ..., verbose = getOption("verbose", FALSE))

Arguments

t0

The zero-time, i.e., the time of the initial state. This must be no later than the time of the first observation, i.e., t0 <= times[1].

times

the sequence of observation times. times must indicate the column of observation times by name or index. The time vector must be numeric and non-decreasing.

params

a named numeric vector or a matrix with rownames containing the parameters at which the simulations are to be performed.

skeleton

optional; the deterministic skeleton of the unobserved state process. Depending on whether the model operates in continuous or discrete time, this is either a vectorfield or a map. Accordingly, this is supplied using either the vectorfield or map fnctions. For more information, see skeleton specification. Setting skeleton=NULL removes the deterministic skeleton.

rinit

simulator of the initial-state distribution. This can be furnished either as a C snippet, an R function, or the name of a pre-compiled native routine available in a dynamically loaded library. Setting rinit=NULL sets the initial-state simulator to its default. For more information, see rinit specification.

...

additional arguments are passed to pomp.

ode_control

optional list; the elements of this list will be passed to ode if the skeleton is a vectorfield, and ignored if it is a map.

format

the format in which to return the results.

format = "pomps" causes the trajectories to be returned as a single ‘pomp’ object (if a single parameter vector has been furnished to trajectory) or as a ‘pompList’ object (if a matrix of parameters have been furnished). In each of these, the states slot will have been replaced by the computed trajectory. Use states to view these.

format = "array" causes the trajectories to be returned in a rank-3 array with dimensions nvar x ncol(params) x ntimes. Here, nvar is the number of state variables and ntimes the length of the argument times. Thus if x is the returned array, x[i,j,k] is the i-th component of the state vector at time times[k] given parameters params[,j].

format = "data.frame" causes the results to be returned as a single data frame containing the time and states. An ordered factor variable, ‘.id’, distinguishes the trajectories from one another.

verbose

logical; if TRUE, diagnostic messages will be printed to the console.

object

optional; if present, it should be a data frame or a ‘pomp’ object.

Details

In the case of a discrete-time system, the deterministic skeleton is a map and a trajectory is obtained by iterating the map. In the case of a continuous-time system, the deterministic skeleton is a vector-field; trajectory uses the numerical solvers in deSolve to integrate the vectorfield.

Value

The format option controls the nature of the return value of trajectory. See above for details.

See Also

More on pomp elementary algorithms: elementary_algorithms, kalman, pfilter(), pomp-package, probe(), simulate(), spect(), wpfilter()

More on methods for deterministic process models: flow(), skeleton(), skeleton_spec, traj_match

Examples

  ## The basic components needed to compute trajectories
  ## of a deterministic dynamical system are
  ## rinit and skeleton.

  ## The following specifies these for a simple continuous-time
  ## model: dx/dt = r (1+e cos(t)) x

  trajectory(
    t0 = 0, times = seq(1,30,by=0.1),
    rinit = function (x0, ...) {
      c(x = x0)
    },
    skeleton = vectorfield(
      function (r, e, t, x, ...) {
        c(x=r*(1+e*cos(t))*x)
      }
    ),
    params = c(r=1,e=3,x0=1)
  ) -> po

  plot(po,log='y')

  ## In the case of a discrete-time skeleton,
  ## we use the 'map' function.  For example,
  ## the following computes a trajectory from
  ## the dynamical system with skeleton
  ## x -> x exp(r sin(omega t)).

  trajectory(
    t0 = 0, times=seq(1,100),
    rinit = function (x0, ...) {
      c(x = x0)
    },
    skeleton = map(
      function (r, t, x, omega, ...) {
        c(x=x*exp(r*sin(omega*t)))
      },
      delta.t=1
    ),
    params = c(r=1,x0=1,omega=4)
  ) -> po

  plot(po)


 # takes too long for R CMD check
  ## generate a bifurcation diagram for the Ricker map
  p <- parmat(coef(ricker()),nrep=500)
  p["r",] <- exp(seq(from=1.5,to=4,length=500))
  trajectory(
    ricker(),
    times=seq(from=1000,to=2000,by=1),
    params=p,
    format="array"
  ) -> x
  matplot(p["r",],x["N",,],pch='.',col='black',
    xlab=expression(log(r)),ylab="N",log='x')