pomp (version 0.21-3)

pomp: Partially-observed Markov process object.

Description

Create a new pomp object.

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

The 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

Function of prototype rprocess(xstart,times,params,...) which simulates from the unobserved process.

dprocess

Function of prototype dprocess(x,times,params,log,...) which evaluates the likelihood of a sequence of consecutive state transitions.

rmeasure

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

dmeasure

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

measurement.model

An optional 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

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 th

skeleton.vectorfield

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

initializer

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 i

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 stored in a slot userdata and are passed as arguments to each of the functions rprocess, dprocess, rmeasure, dmeasure, and initializer whenever they

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.

Details

It is the user's responsibility to ensure that the rprocess, dprocess, rmeasure, dmeasure, and initializer functions satisfy the following conditions: [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

## 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)
rw2 <- examples[[1]]  ## by default, simulate produces a list of pomp objects
plot(rw2)

t <- seq(0,20)
X <- simulate(rw2,params=p[,1],nsim=10,states=TRUE,obs=TRUE,times=t)

## compute the process model likelihoods
f <- dprocess(
              rw2,
              x=X$states,
              times=t,
              params=matrix(
                p[,1],
                nrow=nrow(p),
                ncol=10,
                dimnames=list(rownames(p),NULL)
                ),
              log=TRUE
              )
apply(f,1,sum)

## compute the measurement likelihoods
g <- dmeasure(
              rw2,
              y=X$obs[,7,],
              x=X$states,
              times=t,
              params=matrix(
                p[,1],
                nrow=nrow(p),
                ncol=10,
                dimnames=list(rownames(p),NULL)
                ),
              log=TRUE
              )
apply(g,1,sum)

## For more examples, see the vignettes.