rprocess_spec: rprocess specification

In pomp: Statistical Inference for Partially Observed Markov Processes

rprocess specification

Description

Specification of the latent state process simulator, rprocess.

Usage

onestep(step.fun)

discrete_time(step.fun, delta.t = 1)

euler(step.fun, delta.t)

gillespie(rate.fun, v, hmax = Inf)

gillespie_hl(..., .pre = "", .post = "", hmax = Inf)

Arguments

step.fun	a C snippet, an R function, or the name of a native routine in a shared-object library. This gives a procedure by which one simulates a single step of the latent state process.
delta.t	positive numerical value; for euler and discrete_time, the size of the step to take
rate.fun	a C snippet, an R function, or the name of a native routine in a shared-object library. This gives a procedure by which one computes the event-rate of the elementary events in the continuous-time latent Markov chain.
v	integer matrix; giving the stoichiometry of the continuous-time latent Markov process. It should have dimensions nvar x nevent, where nvar is the number of state variables and nevent is the number of elementary events. v describes the changes that occur in each elementary event: it will usually comprise the values 1, -1, and 0 according to whether a state variable is incremented, decremented, or unchanged in an elementary event. The rows of v may be unnamed or named. If the rows are unnamed, they are assumed to be in the same order as the vector of state variables returned by rinit. If the rows are named, the names of the state variables returned by rinit will be matched to the rows of v to ensure a correct mapping. If any of the row names of v cannot be found among the state variables or if any row names of v are duplicated, an error will occur.
hmax	maximum time step allowed (see below)
...	individual C snippets corresponding to elementary events
.pre, .post	C snippets (see Details)

Discrete-time processes

If the state process evolves in discrete time, specify rprocess using the discrete_time plug-in. Specifically, provide

    rprocess = discrete_time(step.fun = f, delta.t),

where f is a C snippet or R function that simulates one step of the state process. The former is the preferred option, due to its much greater computational efficiency. The goal of such a C snippet is to replace the state variables with their new random values at the end of the time interval. Accordingly, each state variable should be over-written with its new value. In addition to the states, parameters, covariates (if any), and observables, the variables t and dt, containing respectively the time at the beginning of the step and the step's duration, will be defined in the context in which the C snippet is executed. See Csnippet for general rules on writing C snippets. Examples are to be found in the tutorials on the package website.

If f is given as an R function, its arguments should come from the state variables, parameters, covariates, and time. It may also take the argument ‘delta.t’; when called, the latter will be the timestep. It must also have the argument ‘...’. It should return a named vector of length equal to the number of state variables, representing a draw from the distribution of the state process at time t+delta.t conditional on its value at time t.

Continuous-time processes

If the state process evolves in continuous time, but you can use an Euler approximation, implement rprocess using the euler plug-in. Specify

    rprocess = euler(step.fun = f, delta.t)

in this case. As before, f can be provided either as a C snippet or as an R function, the former resulting in much quicker computations. The form of f will be the same as above (in the discrete-time case).

If you have a procedure that allows you, given the value of the state process at any time, to simulate it at an arbitrary time in the future, use the onestep plug-in. To do so, specify

    rprocess = onestep(step.fun = f).

Again, f can be provided either as a C snippet or as an R function, the former resulting in much quicker computations. The form of f should be as above (in the discrete-time or Euler cases).

Size of time step

The simulator plug-ins discrete_time, euler, and onestep all work by taking discrete time steps. They differ as to how this is done. Specifically,

1. onestep takes a single step to go from any given time t1 to any later time t2 (t1 <= t2). Thus, this plug-in is designed for use in situations where a closed-form solution to the process exists.

2. To go from t1 to t2, euler takes n steps of equal size, where

       n = ceiling((t2-t1)/delta.t).
   

3. discrete_time assumes that the process evolves in discrete time, where the interval between successive times is delta.t. Thus, to go from t1 to t2, discrete_time takes n steps of size exactly delta.t, where

       n = floor((t2-t1)/delta.t).
   

Exact (event-driven) simulations

If you desire exact simulation of certain continuous-time Markov chains, an implementation of Gillespie's algorithm (Gillespie 1977) is available, via the gillespie and gillespie_hl plug-ins. The former allows for the rate function to be provided as an R function or a single C snippet, while the latter provides a means of specifying the elementary events via a list of C snippets.

A high-level interface to the simulator is provided by gillespie_hl. To use it, supply

    rprocess = gillespie_hl(..., .pre = "", .post = "", hmax = Inf)

to pomp. Each argument in ... corresponds to a single elementary event and should be a list containing two elements. The first should be a string or C snippet; the second should be a named integer vector. The variable rate will exist in the context of the C snippet, as will the parameter, state variables, covariates, and the time t. The C snippet should assign to the variable rate the corresponding elementary event rate.

The named integer vector specifies the changes to the state variables corresponding to the elementary event. There should be named value for each of the state variables returned by rinit. The arguments .pre and .post can be used to provide C code that will run respectively before and after the elementary-event snippets. These hooks can be useful for avoiding duplication of code that performs calculations needed to obtain several of the different event rates.

Here's how a simple birth-death model might be specified:

    gillespie_hl(
        birth=list("rate = b*N;",c(N=1)),
        death=list("rate = m*N;",c(N=-1))
    )

In the above, the state variable N represents the population size and parameters b, m are the birth and death rates, respectively.

To use the lower-level gillespie interface, furnish

    rprocess = gillespie(rate.fun = f, v, hmax = Inf)

to pomp, where f gives the rates of the elementary events. Here, f may be an R function of the form

    f(j, x, t, params, ...)

When f is called, the integer j will be the number of the elementary event (corresponding to the column the matrix v, see below), x will be a named numeric vector containing the value of the state process at time t and params is a named numeric vector containing parameters. f should return a single numerical value, representing the rate of that elementary event at that point in state space and time.

Here, the stoichiometric matrix v specifies the continuous-time Markov process in terms of its elementary events. It should have dimensions nvar x nevent, where nvar is the number of state variables and nevent is the number of elementary events. v describes the changes that occur in each elementary event: it will usually comprise the values 1, -1, and 0 according to whether a state variable is incremented, decremented, or unchanged in an elementary event. The rows of v should have names corresponding to the state variables. If any of the row names of v cannot be found among the state variables or if any row names of v are duplicated, an error will occur.

It is also possible to provide a C snippet via the rate.fun argument to gillespie. Such a snippet should assign the correct value to a rate variable depending on the value of j. The same variables will be available as for the C code provided to gillespie_hl. This lower-level interface may be preferable if it is easier to write code that calculates the correct rate based on j rather than to write a snippet for each possible value of j. For example, if the number of possible values of j is large and the rates vary according to a few simple rules, the lower-level interface may provide the easier way of specifying the model.

When the process is non-autonomous (i.e., the event rates depend explicitly on time), it can be useful to set hmax to the maximum step that will be taken. By default, the elementary event rates will be recomputed at least once per observation interval.

Default behavior

The default rprocess is undefined. It will yield missing values (NA) for all state variables.

Note for Windows users

Some Windows users report problems when using C snippets in parallel computations. These appear to arise when the temporary files created during the C snippet compilation process are not handled properly by the operating system. To circumvent this problem, use the cdir and cfile options to cause the C snippets to be written to a file of your choice, thus avoiding the use of temporary files altogether.

See Also

rprocess

More on implementing POMP models: Csnippet, accumvars, basic_components, betabinomial, covariates, dinit_spec, dmeasure_spec, dprocess_spec, emeasure_spec, eulermultinom, parameter_trans(), pomp-package, pomp_constructor, prior_spec, rinit_spec, rmeasure_spec, skeleton_spec, transformations, userdata, vmeasure_spec

pomp documentation built on Sept. 13, 2024, 1:08 a.m.