pomp: Partially-observed Markov process object.

Description

Create a new pomp object to hold a partially-observed Markov process model together with a uni- or multi-variate time series.

Usage

pomp(data, times, t0, ..., rprocess, dprocess, rmeasure, dmeasure,
       measurement.model, skeleton.map, skeleton.vectorfield,
       initializer, covar, tcovar, statenames, paramnames, covarnames,
       PACKAGE)

Arguments

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

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.

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

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,

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 specif

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,...) whic

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.

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 c

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 eith

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,

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.

`Value`

An object of class pomp.

`Warning`

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.

`Details`

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.
  [object Object],[object Object],[object Object],[object Object],[object Object],[object Object]

`See Also`

pomp-class,
  pomp-methods,
  rprocess,
  dprocess,
  rmeasure,
  dmeasure,
  skeleton,
  init.state,
  euler

`Examples`

Run this code## 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.
Run the code above in your browser using DataLab