pomp (version 0.28-2)

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)

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

## 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.