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

userdata

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, 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:

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

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

pomp documentation built on Aug. 8, 2023, 1:08 a.m.