userdata: Facilities for making additional information to basic... In pomp: Statistical Inference for Partially Observed Markov Processes

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

 `1` ``` 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:

1. as an R function,

2. as a C snippet, or

3. 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:

 ```1 2 3 4``` ``` 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:

 `1` ``` simulate(..., rmeasure = f, theta = 42, ...) ```

where the `...` represent the other `simulate` arguments we might want to supply. When we do so, a message will be generated, informing us that `theta` is available for use by the POMP basic components. This warning helps forestall accidental triggering of this facility due to typographical error.

When the basic component is specified via a C snippet

A C snippet implementation of the aforementioned measurement model is:

 ```1 2 3 4``` ``` f <- Csnippet(" 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 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

 `1` ``` simulate(..., rmeasure = f, 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 faster and more elegant, when using C snippets, to put the needed objects directly into the C snippet library. The `globals` argument does this. See the example below.

More on implementing POMP models: `Csnippet`, `accumulators`, `basic_components`, `covariate_table()`, `distributions`, `dmeasure_spec`, `dprocess_spec`, `parameter_trans()`, `pomp-package`, `prior_spec`, `rinit_spec`, `rmeasure_spec`, `rprocess_spec`, `skeleton_spec`, `transformations`
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59``` ```## The familiar Ricker example ## For some bizarre reason, we wish to pass 'phi' via the userdata facility. ## 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(" double e = (sigma > 0.0) ? rnorm(0,sigma) : 0.0; N = r*N*exp(-N+e);" ), delta.t=1 ), rmeasure=Csnippet(" 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(" 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 ```