rprocess_spec {pomp} | R Documentation |
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 |
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 |
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,
-
onestep
takes a single step to go from any given timet1
to any later timet2
(t1 <= t2
). Thus, this plug-in is designed for use in situations where a closed-form solution to the process exists. To go from
t1
tot2
,euler
takesn
steps of equal size, wheren = ceiling((t2-t1)/delta.t).
-
discrete_time
assumes that the process evolves in discrete time, where the interval between successive times isdelta.t
. Thus, to go fromt1
tot2
,discrete_time
takesn
steps of size exactlydelta.t
, wheren = 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
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