Constructor of the basic POMP object

Description

This function constructs a pomp object, encoding a partially-observed Markov process model together with a uni- or multivariate time series. One implements the model by specifying its components, each of which can be written as R functions or, for much greater computational efficiency, using C code. The preferred way to specify most components (as detailed below) is through the use of Csnippets, snippets of C that are compiled and linked into a running R session.

Usage

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
## S4 method for signature 'data.frame'
pomp(data, times, t0, ..., rprocess, dprocess, rmeasure, dmeasure,
       measurement.model,
       skeleton, skeleton.type = c("map","vectorfield"), skelmap.delta.t = 1,
       initializer, rprior, dprior, params, covar, tcovar,
       obsnames, statenames, paramnames, covarnames, zeronames,
       PACKAGE, fromEstimationScale, toEstimationScale, globals)
## S4 method for signature 'numeric'
pomp(data, times, t0, ..., rprocess, dprocess, rmeasure, dmeasure,
       measurement.model,
       skeleton, skeleton.type = c("map","vectorfield"), skelmap.delta.t = 1,
       initializer, rprior, dprior, params, covar, tcovar,
       obsnames, statenames, paramnames, covarnames, zeronames,
       PACKAGE, fromEstimationScale, toEstimationScale, globals)
## S4 method for signature 'matrix'
pomp(data, times, t0, ..., rprocess, dprocess, rmeasure, dmeasure,
       measurement.model,
       skeleton, skeleton.type = c("map","vectorfield"), skelmap.delta.t = 1,
       initializer, rprior, dprior, params, covar, tcovar,
       obsnames, statenames, paramnames, covarnames, zeronames,
       PACKAGE, fromEstimationScale, toEstimationScale, globals)
## S4 method for signature 'pomp'
pomp(data, times, t0, ..., rprocess, dprocess, rmeasure, dmeasure,
       measurement.model, skeleton, skeleton.type, skelmap.delta.t,
       initializer, rprior, dprior, params, covar, tcovar,
       obsnames, statenames, paramnames, covarnames, zeronames,
       PACKAGE, fromEstimationScale, toEstimationScale, globals)

Arguments

data, times

The time series data and times at which observations are made. data can be specified as a vector, a matrix, a data-frame. Alternatively, a pomp object can be supplied in the data argument.

If data is a numeric vector, times must be a numeric vector of the same length.

If data is a matrix, it should have dimensions nobs x ntimes, where nobs is the number of observed variables and ntimes is the number of times at which observations were made (i.e., each column is a distinct observation of the nobs variables). In this case, times must be given as a numeric vector (of length ntimes).

If data is a data-frame, times must name the column of observation times. Note that, in this case, data will be internally coerced to an array with storage-mode ‘double’.

times must be numeric and strictly increasing.

t0

The zero-time, at which the stochastic dynamical system is to be initialized. This must be no later than the time of the first observation, i.e., t0 <= times[1].

rprocess

optional function; a function of prototype

rprocess(xstart,times,params,\dots)

that simulates from the unobserved process. The form of this function is given below. pomp provides a number of plugins that construct appropriate rprocess arguments corresponding to several common stochastic simulation algorithms. See below for more details.

dprocess

optional function; a function of prototype

dprocess(x,times,params,log,\dots)

that evaluates the likelihood of a sequence of consecutive state transitions. The form of this function is given below. It is not typically necessary (or even feasible) to define dprocess. This functionality is provided to support future algorithm development.

rmeasure

optional; the measurement model simulator. This can be specified in one of four ways:

  1. as a function of prototype

    rmeasure(x,t,params,\dots)

    that makes a draw from the observation process given states x, time t, and parameters params.

  2. as the name of a native (compiled) routine of type

    pomp_measure_model_simulator

    as defined in the header file ‘pomp.h’. (To view the header file, execute

    file.show(system.file("include/pomp.h",package="pomp"))

    in an R session.)

  3. using the formula-based measurement.model facility (see below).

  4. as a snippet of C code (via Csnippet) that draws from the observation process as above. The last is typically the preferred option, as it results in much faster code execution.

dmeasure

optional; the measurement model probability density function. This can be specified in one of four ways:

  1. as a function of prototype

    dmeasure(y,x,t,params,log,\dots)

    that computes the p.d.f. of y given x, t, and params.

  2. as the name of a native (compiled) routine of type

    pomp_measure_model_density

    as defined in the header file ‘pomp.h’. (To view the header file, execute

    file.show(system.file("include/pomp.h",package="pomp"))

    in an R session.)

  3. using the formula-based measurement.model facility (see below).

  4. as a snippet of C code (via Csnippet) that computes the p.d.f. as above.

The last is typically the preferred option, as it results in much faster code execution. As might be expected, if log=TRUE, this function should return the log likelihood.

measurement.model

optional; a formula or list of formulae, specifying the measurement model. These formulae are parsed internally to generate rmeasure and dmeasure functions. If measurement.model is given it overrides any specification of rmeasure or dmeasure. NB: This is a convenience function, primarily designed to facilitate exploration; it will typically be possible to acclerate measurement model computations by writing dmeasure and/or rmeasure using Csnippets.

skeleton, skeleton.type, skelmap.delta.t

The function skeleton specifies the deterministic skeleton of the unobserved Markov process. If we are dealing with a discrete-time Markov process, its deterministic skeleton is a map: indicate this by specifying skeleton.type="map". In this case, the default assumption is that time advances 1 unit per iteration of the map; to change this, set skelmap.delta.t to the appropriate time-step. If we are dealing with a continuous-time Markov process, its deterministic skeleton is a vectorfield: indicate this by specifying skeleton.type="vectorfield".

The skeleton function can be specified in one of three ways:

  1. as an R function of prototype

    skeleton(x,t,params,\dots)

    that evaluates the deterministic skeleton at state x and time t given the parameters params,

  2. as the name of a native (compiled) routine of type

    pomp_skeleton

    as defined in the header file ‘pomp.h’. (To view the header file, execute

    file.show(system.file("include/pomp.h",package="pomp"))

    in an R session.)

  3. as a snippet of C code (via Csnippet) that performs this evaluation. The latter is typically the preferred option, for reasons of computational efficiency.

initializer

The initializer gives the parameterization of the initial state of the unobserved Markov process. Specifically, given a vector of parameters, params and an initial time, t0, the initializer determines the state vector at time t0.

By default, any parameters in params, the names of which end in “.0”, are assumed to be initial values of states. To initialize the unobserved state process, these are simply copied over as initial conditions. The names of the resulting state variables are obtained by dropping the “.0” suffix.

A custom initializer can be supplied here in one of two formats:

  1. as an R function of prototype

    initializer(params,t0,\dots)

    that yields initial conditions for the state process when given a vector, params, of parameters.

  2. as a snippet of C code (via Csnippet). See the Csnippet help for rules on writing a custom initializer.

rprior

optional; function drawing a sample from a prior distribution on parameters. This can be specified in one of three ways:

  1. as an R function of prototype

    rprior(params,\dots)

    that makes a draw from the prior distribution given params,

  2. as the name of a native (compiled) routine of type

    pomp_rprior

    as defined in the header file ‘pomp.h’, or (To view the header file, execute

    file.show(system.file("include/pomp.h",package="pomp"))

    in an R session.)

  3. as a snippet of C code (via Csnippet).

As above, the latter is typically preferable.

dprior

optional; function evaluating the prior distribution. This can be specified in one of three ways:

  1. as an R function of prototype

    dprior(params,log=FALSE,\dots)

    that evaluates the prior probability density,

  2. as the name of a native (compiled) routine of type

    pomp_dprior

    as defined in the header file ‘pomp.h’, or (To view the header file, execute

    file.show(system.file("include/pomp.h",package="pomp"))

    in an R session.)

  3. as a snippet of C code (via Csnippet).

As above, the latter is typically preferable.

params

optional named numeric vector of parameters. This will be coerced internally to storage mode double.

covar, tcovar

An optional matrix or data frame of covariates: covar is the table of covariates (one column per variable); tcovar the corresponding times (one entry per row of covar).

covar can be specified as either a matrix or a data frame. In either case the columns are taken to be distinct covariates. If covar is a data frame, tcovar can be either the name or the index of the time variable.

If a covariate table is supplied, then the value of each of the covariates is interpolated as needed. The resulting interpolated values are passed to the corresponding functions as a numeric vector named covars; see below.

obsnames, statenames, paramnames, covarnames

Optional character vectors specifying the names of observables, state variables, parameters, and covariates, respectively. These are only used in the event that one or more of the basic functions (rprocess, dprocess, rmeasure, dmeasure, skeleton, rprior, dprior) are defined using Csnippet or native routines.

zeronames

optional character vector specifying the names of accumulator variables (see below).

PACKAGE

An optional string giving the name of the dynamically loaded library in which any native routines are to be found.

fromEstimationScale, toEstimationScale

Optional functions specifying parameter transformations. Many algorithms for parameter estimation search an unconstrained space of parameters. When working with such an algorithm and a model for which the parameters are constrained, it can be useful to transform parameters. toEstimationScale and fromEstimationScale are transformations from the model scale to the estimation scale, and vice versa, respectively. These functions must have arguments params and .... toEstimationScale should transform parameters from the scale that rprocess, dprocess, rmeasure, dmeasure, and skeleton use internally to the scale used in estimation. fromEstimationScale should be the inverse of toEstimationScale. The parameter transformations can be defined (as above) using either R functions, native routines, or Csnippets.

Note that it is the user's responsibility to make sure that these transformations are mutually inverse. If obj is the constructed pomp object, and coef(obj) is non-empty, a simple check of this property is

      x <- coef(obj,transform=TRUE)
      obj1 <- obj
      coef(obj1,transform=TRUE) <- x
      identical(coef(obj),coef(obj1))
      identical(coef(obj1,transform=TRUE),x).
    

By default, both functions are the identity transformation. See the demos,

demo(package="pomp"),

pompExample, and the tutorials on the package website for examples.

globals

optional character; C code that will be included in the source for (and therefore hard-coded into) the shared-object library created when the call to pomp uses Csnippets. If no Csnippets are used, globals has no effect.

...

Any additional arguments given to pomp will be stored in the pomp object and passed as arguments to each of the basic functions whenever they are evaluated.

Value

pomp returns an object of class pomp. If data is an object of class pomp, then by default the returned pomp object is identical to data. If additional arguments are given, these override the defaults.

Important note

It is not typically necessary (or even feasible) to define all of the components rprocess, dprocess, rmeasure, dmeasure, and skeleton in any given problem. Each algorithm makes use of only a subset of these components. Any algorithm requiring a component that has not been defined will return an informative error.

The state process model

Specification of process-model codes rprocess and/or dprocess in most cases is facilitated by pomp's so-called plugins, which have been developed to handle common use-cases. Currently, if one's process model evolves in discrete time or one is willing to make such an approximation (e.g., via an Euler approximation), then the euler.sim, discrete.time.sim, and onestep.sim plugins for rprocess and onestep.dens plugin for dprocess are available. In addition, for exact simulation of certain continuous-time Markov chains, an implementation of Gillespie's algorithm is available (see gillespie.sim). To learn more about the use of plugins, consult the help documentation (plugins) and the tutorials on the package website. Several of the demos and examples make use of these as well.

In specific cases, it may be possible to obtain increased computational efficiency by writing custom versions of rprocess and/or dprocess instead of using the plugins. If such custom versions are desired, the following describes how these functions should be written.

rprocess

If the plugins are not used rprocess must be an R function with at least the following arguments: xstart, times, params, and .... It can also take additional arguments. It is guaranteed that these will be filled with the corresponding elements the user has included as additional arguments in the construction of the pomp object.

In calls to rprocess, xstart can be assumed to be an nvar x nrep matrix; its rows correspond to components of the state vector and columns correspond to independent realizations of the process. params will similarly be an npar x nrep matrix with rows corresponding to parameters and columns corresponding to independent realizations. Note that the columns of params correspond to those of xstart; in particular, they will agree in number. Both xstart and params are guaranteed to have rownames.

rprocess must return a rank-3 array with rownames. Suppose x is the array returned. Then dim(x)=c(nvars,nrep,ntimes), where ntimes is the length of the vector times. x[,j,k] is the value of the state process in the j-th realization at time times[k]. In particular, x[,,1] must be identical to xstart. The rownames of x must correspond to those of xstart.

dprocess

If the plugins are not used, dprocess must have at least the following arguments: x, times, params, log, and .... It may take additional arguments: again, these will be filled with the corresponding elements the user defines when the pomp object is constructed.

In calls to dprocess, x may be assumed to be an nvars x nrep x ntimes array, where these terms have the same meanings as above. params will be a matrix with rows corresponding to individual parameters and columns corresponding to independent realizations. The columns of params correspond to those of x; in particular, they will agree in number. Both x and params are guaranteed to have rownames.

dprocess must return a nrep x ntimes-1 matrix. Suppose d is the array returned. d[j,k] is the probability density of the transition from state x[,j,k-1] at time times[k-1] to state x[,j,k] at time times[k]. If log=TRUE, then the log of the pdf must be returned.

In writing this function, you may assume that the transitions are consecutive. It should be clear that, but for this assumption, it will in general be impossible to write the transition probabilities explicitly. In such cases, algorithms that make no use of dprocess, which are said to have the “plug and play” property, are useful. Most of the algorithms in pomp have this property. In particular, at present, no methods in pomp make use of dprocess.

The observation process model

The following is a guide to writing the measurement model components as R functions. For a description on how to write these components using Csnippets, see the tutorials on the package website.

rmeasure

if provided, must take at least the arguments x, t, params, and .... It may take additional arguments, which will be filled with user-specified data as above. x will be a named numeric vector of length nvars (which has the same meaning as above). t will be a scalar quantity, the time at which the measurement is made. params will be a named numeric vector of length npars. The rmeasure function may take additional arguments which will be filled with user-specified data as above.

rmeasure must return a named numeric vector of length nobs, the number of observable variables.

dmeasure

if provided, must take at least the arguments y, x, t, params, log, and .... y will be a named numeric vector of length nobs containing (actual or simulated) values of the observed variables; x will be a named numeric vector of length nvar containing state variables; params will be a named numeric vector containing parameters; and t will be a scalar, the corresponding observation time. The dmeasure function may take additional arguments which will be filled with user-specified data as above.

dmeasure must return a single numeric value, the probability density of y given x at time t. If log=TRUE, then dmeasure should return the log of the probability density.

The deterministic skeleton

The following describes how to specify the deterministic skeleton as an R function. For a description on how to write this component using Csnippets, see the tutorials on the package website and the Csnippet help.

If skeleton if provided, must have at least the arguments x, t, params, and .... x is a numeric vector containing the coordinates of a point in state space at which evaluation of the skeleton is desired. t is a numeric value giving the time at which evaluation of the skeleton is desired. Of course, these will be irrelevant in the case of an autonomous skeleton. params is a numeric vector holding the parameters. skeleton may take additional arguments, which will be filled, as above, with user-specified data.

skeleton must return a numeric vector of the same length as x, which contains the value vectorfield (if the dynamical system is continuous) or the value of the map (if the dynamical system is discrete), at the point x at time t.

The state-process initializer

if provided, must have at least the arguments params, t0, and .... params will be a named numeric vector of parameters. t0 will be the time at which initial conditions are desired. initializer must return a named numeric vector of initial states.

Covariates

If the pomp object contains covariates (via the covar argument; see above), then whenever any of the R functions described above are called, they will each be supplied with an additional argument covars. This will be a named numeric vector containing the (interpolated) values of the covariates at the time t. In particular, covars will have one value for each column of the covariate table.

Accumulator variables

In formulating models, one often wishes to define a state variable that will accumulate some quantity over the interval between successive observations. pomp provides a facility to make such features more convenient. Specifically, variables named in the pomp's zeronames argument will be set to zero immediately following each observation. See euler.sir and the tutorials on the package website for examples.

Warning

Some error checking is done by pomp, but complete error checking is impossible. If the user-specified functions do not conform to the above specifications, then the results may be invalid. In particular, if both rmeasure and dmeasure are specified, the user should verify that these two functions correspond to the same probability distribution. If skeleton is specified, the user is responsible for verifying that it corresponds to a deterministic skeleton of the model. Each pomp-package algorithm uses some subset of the five basic functions (rprocess, dprocess, rmeasure, dmeasure, skeleton). If an algorithm requires a component that has not been specified, an informative error will be generated.

Author(s)

Aaron A. King kingaa at umich dot edu

See Also

pomp methods, pomp low-level interface, process model plugins

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
## Not run: 
pompExample()
pomp.home <- system.file("examples",package="pomp")
pomp.examples <- list.files(pomp.home)
file.show(
          file.path(pomp.home,pomp.examples),
          header=paste("======",pomp.examples,"=======")
         )

## End(Not run)

Want to suggest features or report bugs for rdrr.io? Use the GitHub issue tracker.