trajectory: Trajectory of a deterministic model

Description

Compute trajectories of the deterministic skeleton of a Markov process.

Usage

# S4 method for missing
trajectory(
  ...,
  t0,
  times,
  params,
  skeleton,
  rinit,
  ode_control = list(),
  format = c("pomps", "array", "data.frame"),
  verbose = getOption("verbose", FALSE)
)
# S4 method for 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 pomp
trajectory(
  object,
  ...,
  params,
  skeleton,
  rinit,
  ode_control = list(),
  format = c("pomps", "array", "data.frame"),
  verbose = getOption("verbose", FALSE)
)
# S4 method for traj_match_objfun
trajectory(object, ..., verbose = getOption("verbose", FALSE))

Value

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

Arguments

...

additional arguments are passed to pomp.

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.

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.

Examples

Run this code

# \donttest{
  ## 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')