pomp {pomp} | R Documentation |
Create a new pomp
object to hold a partially-observed Markov process model together with a uni- or multi-variate time series.
pomp(data, times, t0, ..., rprocess, dprocess, rmeasure, dmeasure, measurement.model, skeleton.map, skeleton.vectorfield, initializer, covar, tcovar, statenames, paramnames, covarnames, PACKAGE)
data |
An array holding the data.
This array is of dimensions nobs x ntimes , where nobs is the number of observed variables and ntimes is the number of times at which observations were made.
It is also possible to specify data as a data-frame, in which case the argument times must be the name of the column of observation times.
Note that if data is a data-frame, it will be internally coerced to an array with storage-mode ‘double’.
|
times |
vector of times corresponding to the observations:
times must be a strictly increasing numeric vector.
If data is a data-frame, times should be the name of the column of observation times.
|
t0 |
The zero-time.
This must be no later than the time of the first observation, times[1] .
|
rprocess |
optional; a function of prototype rprocess(xstart,times,params,...) which simulates from the unobserved process.
See below for details.
|
dprocess |
optional; a function of prototype dprocess(x,times,params,log,...) which evaluates the likelihood of a sequence of consecutive state transitions.
See below for details.
|
rmeasure |
optional; the measurement model simulator.
This can be specified in one of three ways:
(1) as a function of prototype rmeasure(x,t,params,...) which makes a draw from the observation process given states x , time t , and parameters params .
(2) as the name of a native (compiled) routine with prototype “pomp_measure_model_simulator” as defined in the header file “pomp.h”.
In the above cases, if the measurement model depends on covariates, the optional argument covars will be filled with interpolated values at each call.
(3) using the formula-based measurement.model facility (see below).
|
dmeasure |
optional; the measurement model probability density function.
This can be specified in one of three ways:
(1) as a function of prototype dmeasure(y,x,t,params,log,...) which computes the p.d.f. of y given x , t , and params .
(2) as the name of a native (compiled) routine with prototype “pomp_measure_model_density” as defined in the header file “pomp.h”.
In the above cases, if the measurement model depends on covariates, the optional argument covars will be filled with interpolated values at each call.
(3) using the formula-based measurement.model facility (see below).
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 and used to generate rmeasure and dmeasure functions.
If measurement.model is given it overrides any specification of rmeasure or dmeasure .
See below for an example.
NB: it will typically be possible to acclerate measurement model computations by writing dmeasure and/or rmeasure functions directly.
|
skeleton.map |
optional.
If we are dealing with a discrete-time Markov process, its deterministic skeleton is a map and can be specified using skeleton.map in one of two ways:
(1) as a function of prototype skeleton(x,t,params,...) which evaluates the deterministic skeleton at state x and time t given the parameters params , or
(2) as the name of a native (compiled) routine with prototype “pomp_vectorfield_map” as defined in the header file “pomp.h”.
If the deterministic skeleton depends on covariates, the optional argument covars will be filled with interpolated values of the covariates at the time t .
|
skeleton.vectorfield |
optional.
If we are dealing with a continuous-time Markov process, its deterministic skeleton is a vectorfield and can be specified using skeleton.vectorfield in either of the two ways described above, under skeleton.map .
Note that it makes no sense to specify the deterministic skeleton both as a map and as a vectorfield: an attempt to do so will generate an error.
|
initializer |
optional; a function of prototype initializer(params,t0,...) which yields initial conditions for the state process when given a vector, params , of parameters.
By default, i.e., if it is unspecified when pomp is called, the initializer assumes any names in params that end in “.0” correspond to initial value parameters.
These are simply copied over as initial conditions when init.state is called (see init.state-pomp ).
The names of the state variables are the same as the corresponding initial value parameters, but with the “.0” dropped.
|
covar, tcovar |
An optional table of covariates: covar is the table (with one column per variable) and tcovar the corresponding times (one entry per row of covar ).
This can be in the form of 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, i.e., whenever rprocess , dprocess , rmeasure , dmeasure , or init.state is evaluated.
The resulting interpolated values are passed to the corresponding functions as a numeric vector named covars .
|
statenames, paramnames, covarnames |
Optional character vectors specifying the names of state variables, parameters, or covariates, respectively.
These are only used in the event that one or more of the basic functions (rprocess , dprocess , rmeasure , dmeasure , skeleton ) are defined using native routines.
In that case, these name vectors are matched against the corresponding names and the indices of the names are passed to the native routines.
|
PACKAGE |
An optional string giving the name of the dynamically loaded library in which the native routines are to be found. |
... |
Any additional arguments are passed as arguments to each of the functions rprocess , dprocess , rmeasure , dmeasure , and initializer whenever they are evaluated.
|
It is not typically necessary (or easy) to define all of the functions rprocess
, dprocess
, rmeasure
, dmeasure
, and initializer
in any given problem.
Each algorithm makes use of a different subset of these functions.
The specification of process-model codes rprocess
and/or dprocess
can be somewhat nontrivial:
for this reason, plug-ins have been developed to streamline this process for the user.
At present, if the process model evolves in discrete-time or one is willing to make such an approximation (e.g., via an Euler approximation), then the onestep.simulate
and euler.simulate
plug-ins for rprocess
and onestep.density
plug-in for dprocess
are available.
The following describes how each of these functions should be written.
Examples are provided in the help documentation and the vignettes.
rprocess
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 a rank-2 array (matrix) with rows corresponding to state variables and columns corresponding to independent realizations of the process.
params
will similarly be a rank-2 array with rows corresponding to parameters and columns corresponding to independent realizations.
The columns of params
correspond to those of xstart
;
in particular, they will agree in number.
Both xstart
and params
will have rownames, which are available for use by the user.
rprocess
must return a rank-3 array with rownames.
Suppose x
is the array returned.
Then dim(x)=c(nvars,nreps,ntimes)
, where nvars
is the number of state variables (=nrow(xstart)
), nreps
is the number of independent realizations simulated (=ncol(xstart)
), and 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
.
As mentioned above, some plug-ins are available for rprocess
.
dprocess
x
, times
, params
, log
, and ...
.
It may 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 dprocess
, x
may be assumed to be an nvars
x nreps
x ntimes
array, where these terms have the same meanings as above.
params
will be a rank-2 array 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
will have rownames, available for use by the user.
dprocess
must return a rank-2 array (matrix).
Suppose d
is the array returned.
Then dim(d)=c(nreps,ntimes-1)
.
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 is returned.
As mentioned above, some plug-ins are available for dprocess
.
In writing this function, you may assume that the transitions are consecutive. It should be quite clear that, but for this assumption, it would be quite difficult in general to write the transition probabilities. In fact, from one perspective, the algorithms in pomp are designed to overcome just this difficulty.
rmeasure
x
, t
, params
, and ...
.
It may take additional arguments, which will be filled with user-specified data as above.
x
may be assumed to be a named vector of length nvars
, (which has the same meanings as above).
t
is a scalar quantity, the time at which the measurement is made.
params
may be assumed to be a named vector of length npars
.
rmeasure
must return a named vector.
If y
is the returned vector, then length(y)=nobs
, where nobs
is the number of observable variables.
dmeasure
y
, x
, times
, params
, log
, and ...
.
y
may be assumed to be a named vector of length nobs
containing (actual or simulated) values of the observed variables;
x
will be a named vector of length nvar
containing state variables
params
, a named vector containing parameters;
and t
, a scalar, the corresponding observation time.
It may take additional arguments which will be filled with user-specified data as above.
dmeasure
must return a single numeric value, the pdf of y
given x
at time t
.
If log=TRUE
, then the log of the pdf is returned.
skeleton
skeleton
is an R function, it 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.
The optional argument covars
is a numeric vector containing the values of the covariates at the time t
.
covars
will have one value for each column of the covariate table specified when the pomp
object was created.
covars
is constructed from the covariate table (see covar
, below) by interpolation.
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
.
The return value is interpreted as the 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
.
If skeleton
is the name of a native routine, this routine must be of prototype “pomp_vectorfield_map” as defined in the header “pomp.h” (see the “include” directory in the installed package directory).
initializer
params
, t0
, and ...
.
params
is a named vector of parameters.
t0
will be the time at which initial conditions are desired.
initializer
must return a named vector of initial conditions.
An object of class pomp
.
Some error checking is done, but complete error checking is impossible.
If the user-specified functions do not conform to the above specifications (see Details), then the results may be invalid.
Each algorithm that uses a pomp
object uses some subset of the four basic components (rprocess
, dprocess
, rmeasure
, dmeasure
).
It is unlikely that any algorithm will use all of these.
Aaron A. King kingaa at umich dot edu
pomp-class,
pomp-methods,
rprocess
,
dprocess
,
rmeasure
,
dmeasure
,
skeleton
,
init.state
,
euler
## Simple example: a 2-D Brownian motion, observed with normal error ## first, set up the pomp object: rw2 <- pomp( rprocess = function (xstart, times, params, ...) { ## this function simulates two independent random walks with intensities s1, s2 nreps <- ncol(params) ntimes <- length(times) dt <- diff(times) x <- array(0,dim=c(2,nreps,ntimes)) rownames(x) <- rownames(xstart) noise.sds <- params[c('s1','s2'),] x[,,1] <- xstart for (j in 2:ntimes) { ## we are mimicking a continuous-time process, so the increments have SD ~ sqrt(dt) ## note that we do not have to assume that 'times' are equally spaced x[,,j] <- rnorm(n=2*nreps,mean=x[,,j-1],sd=noise.sds*dt[j-1]) } x }, dprocess = function (x, times, params, log, ...) { ## given a sequence of consecutive states in 'x', this function computes the p.d.f. nreps <- ncol(params) ntimes <- length(times) dt <- diff(times) d <- array(0,dim=c(2,nreps,ntimes-1)) noise.sds <- params[c('s1','s2'),] for (j in 2:ntimes) d[,,j-1] <- dnorm(x[,,j]-x[,,j-1],mean=0,sd=noise.sds*dt[j-1],log=TRUE) if (log) { apply(d,c(2,3),sum) } else { exp(apply(d,c(2,3),sum)) } }, measurement.model=list( y1 ~ norm(mean=x1,sd=tau), y2 ~ norm(mean=x2,sd=tau) ), skeleton = function (x, t, params, covars, ...) { ## since this is a continuous-time process, the skeleton is the vectorfield ## since the random walk is unbiased, the vectorfield is identically zero rep(0,length=length(x)) }, times=1:100, data=rbind( y1=rep(0,100), y2=rep(0,100) ), t0=0 ) ## specify some parameters p <- rbind(s1=c(2,2,3),s2=c(0.1,1,2),tau=c(1,5,0),x1.0=c(0,0,5),x2.0=c(0,0,0)) ## simulate examples <- simulate(rw2,params=p) plot(examples[[1]]) ## construct an (almost) equivalent pomp object using a plugin: rw2 <- pomp( rprocess = euler.simulate, step.fun = function (x, t, params, delta.t, ...) { y <- rnorm(n=2,mean=x[c("x1","x2")],sd=params[c("s1","s2")]*delta.t) names(y) <- c("x1","x2") y }, rmeasure = function (x, t, params, ...) { y <- rnorm(n=2,mean=x[c("x1","x2")],sd=params["tau"]) names(y) <- c("y1","y2") y }, dmeasure = function (y, x, t, params, log, ...) { d <- dnorm(x=y[c("y1","y2")],mean=x[c("x1","x2")],sd=params["tau"],log=TRUE) print(d) if (log) sum(d,na.rm=T) else exp(sum(d,na.rm=T)) }, delta.t=1, data=data.frame( time=1:100, y1=rep(c(1,2,3,4,NA),20), y2=0 ), times="time", t0=0 ) rw2 <- simulate(rw2,params=c(s1=2,s2=0.1,tau=1,x1.0=0,x2.0=0)) ## For more examples, see the vignettes.