userdata {pomp} | R Documentation |
Facilities for making additional information to basic components
Description
When POMP basic components need information they can't get from parameters or covariates.
Details
It can happen that one desires to pass information to one of the POMP model basic components (see here for a definition of this term) outside of the standard routes (i.e., via model parameters or covariates). pomp provides facilities for this purpose. We refer to the objects one wishes to pass in this way as user data.
The following will apply to every basic model component.
For the sake of definiteness, however, we'll use the rmeasure
component as an example.
To be even more specific, the measurement model we wish to implement is
y1 ~ Poisson(x1+theta), y2 ~ Poisson(x2+theta),
where theta
is a parameter.
Although it would be very easy (and indeed far preferable) to include theta
among the ordinary parameters (by including it in params
), we will assume here that we have some reason for not wanting to do so.
Now, we have the choice of providing rmeasure
in one of three ways:
as an R function,
as a C snippet, or
as a procedure in an external, dynamically loaded library.
We'll deal with these three cases in turn.
When the basic component is specified as an R function
We can implement a simulator for the aforementioned measurement model so:
f <- function (t, x, params, theta, ...) { y <- rpois(n=2,x[c("x1","x2")]+theta) setNames(y,c("y1","y2")) }
So far, so good, but how do we get theta
to this function?
We simply provide an additional argument to whichever pomp algorithm we are employing (e.g., simulate
, pfilter
, mif2
, abc
, etc.).
For example:
simulate(..., rmeasure = f, userdata = list(theta = 42), ...)
where the ...
represent other arguments.
When the basic component is specified via a C snippet
A C snippet implementation of the aforementioned measurement model is:
f <- Csnippet(r"{ double theta = *get_userdata_double("theta"); y1 = rpois(x1+theta); y2 = rpois(x2+theta); }")
Here, the call to get_userdata_double
retrieves a pointer to the stored value of theta
.
Note that, by using R string literals (r"{}"
) we avoid the need to escape the quotes in the C snippet text.
It is possible to store and retrieve integer objects also, using get_userdata_int
.
One must take care that one stores the user data with the appropriate storage type.
For example, it is wise to wrap floating point scalars and vectors with as.double
and integers with as.integer
.
In the present example, our call to simulate might look like
simulate(..., rmeasure = f, userdata = list(theta = as.double(42)), ...)
Since the two functions get_userdata_double
and get_userdata_int
return pointers, it is trivial to pass vectors of double-precision and integers.
A simpler and more elegant approach is afforded by the globals
argument (see below).
When the basic component is specified via an external library
The rules are essentially the same as for C snippets.
typedef
declarations for the get_userdata_double
and get_userdata_int
are given in the ‘pomp.h’ header file and these two routines are registered so that they can be retrieved via a call to R_GetCCallable
.
See the Writing R extensions manual for more information.
Setting globals
The use of the userdata facilities incurs a run-time cost.
It is often more efficient, when using C snippets, to put the needed objects directly into the C snippet library.
The globals
argument does this.
See the example below.
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
,
rprocess_spec
,
skeleton_spec
,
transformations
,
vmeasure_spec
Examples
## The familiar Ricker example.
## Suppose that for some reason we wish to pass 'phi'
## via the userdata facility instead of as a parameter.
## C snippet approach:
simulate(times=1:100,t0=0,
phi=as.double(100),
params=c(r=3.8,sigma=0.3,N.0=7),
rprocess=discrete_time(
step.fun=Csnippet(r"{
double e = (sigma > 0.0) ? rnorm(0,sigma) : 0.0;
N = r*N*exp(-N+e);}"
),
delta.t=1
),
rmeasure=Csnippet(r"{
double phi = *get_userdata_double("phi");
y = rpois(phi*N);}"
),
paramnames=c("r","sigma"),
statenames="N",
obsnames="y"
) -> rick1
## The same problem solved using 'globals':
simulate(times=1:100,t0=0,
globals=Csnippet("static double phi = 100;"),
params=c(r=3.8,sigma=0.3,N.0=7),
rprocess=discrete_time(
step.fun=Csnippet(r"{
double e = (sigma > 0.0) ? rnorm(0,sigma) : 0.0;
N = r*N*exp(-N+e);}"
),
delta.t=1
),
rmeasure=Csnippet("
y = rpois(phi*N);"
),
paramnames=c("r","sigma"),
statenames="N",
obsnames="y"
) -> rick2
## Finally, the R function approach:
simulate(times=1:100,t0=0,
phi=100,
params=c(r=3.8,sigma=0.3,N_0=7),
rprocess=discrete_time(
step.fun=function (r, N, sigma, ...) {
e <- rnorm(n=1,mean=0,sd=sigma)
c(N=r*N*exp(-N+e))
},
delta.t=1
),
rmeasure=function (phi, N, ...) {
c(y=rpois(n=1,lambda=phi*N))
}
) -> rick3