dprocess examples?

[sdwfrost]

Hi @kingaa! While I understand that pomp mostly deals with models where the likelihood of the process is unknown, I'd like to use some pseudolikelihood approaches where I specify (or misspecify) dprocess. In going through the old materials on the ou2 example, there is the following code:

ou2 <- function (
  alpha_1 = 0.8, alpha_2 = -0.5, alpha_3 = 0.3, alpha_4 = 0.9,
  sigma_1 = 3, sigma_2 = -0.5, sigma_3 = 2,
  tau = 1,
  x1_0 = -3, x2_0 = 4,
  times = 1:100, t0 = 0
)
{

  simulate(
    times=times, t0=t0,
    seed=787832394,
    rprocess=discrete_time("_ou2_step"),
    dprocess="_ou2_pdf",
    dmeasure = "_ou2_dmeasure",
    rmeasure = "_ou2_rmeasure",
    skeleton = map("_ou2_skel",delta.t=1),
    PACKAGE="pomp",
    paramnames = c(
      "alpha_1","alpha_2","alpha_3","alpha_4",
      "sigma_1","sigma_2","sigma_3",
      "tau"
    ),
    statenames = c("x1","x2"),
    obsnames = c("y1","y2"),
    params=c(
      alpha_1=alpha_1, alpha_2=alpha_2, alpha_3=alpha_3, alpha_4=alpha_4,
      sigma_1=sigma_1, sigma_2=sigma_2, sigma_3=sigma_3,
      tau=tau,
      x1_0=x1_0, x2_0=x2_0
    )
  )
}

This calls the following C code:

// onestep transition probability density for use in 'onestep.dens' plug-in
// transition from x to z as time goes from t1 to t2
void _ou2_pdf (double *f,
  double *x, double *z, double t1, double t2, const double *p,
  const int *stateindex, const int *parindex, const int *covindex,
  const double *covars)
{
  if (t2-t1 != 1)
    errorcall(R_NilValue,"ou2_pdf error: transitions must be consecutive");
  f[0] = dens_ou2(x[X1],x[X2],z[X1],z[X2],ALPHA1,ALPHA2,ALPHA3,ALPHA4,
    SIGMA1,SIGMA2,SIGMA3,1);
}

I have a few questions:

1. How do I rewrite the dprocess in vanilla R code?
2. How do I rewrite the dprocess in a Csnippet.
3. I see that onestep.dens is deprecated. Does this affect things?

I can rewrite rprocess, rmeasure, and dmeasure in either R or Csnippet, but I can't work out the syntax for dprocess for the ou2 example.

[kingaa]

Good to hear from you, @sdwfrost!

Specifying dprocess is somewhat more straightforward in versions 2+ than it was before. The user must write code that evaluates the likelihood of a single transition from arbitrary state x1 at arbitrary time t1 to arbitrary state x2 at arbitrary time t2, where of course t2 >= t1. It can be assumed that the transition is elementary if this is helpful.

You ask three questions; I'll take them in reverse order.

1. onestep.dens is gone as of pomp version 2 because it is no longer needed. Instead, one furnishes the dprocess code directly via the dprocess argument.
2. The preferred way to specify dprocess, as with all other basic model components, if via a C snippet. Access the full documentation by executing ?dprocess_spec. From those docs:

   The goal of a dprocess C snippet is to fill the variable loglik with the log probability density. In the context of such a C snippet, the parameters, and covariates will be defined, as will the times t_1 and t_2. The state variables at time t_1 will have their usual name (see statenames) with a "_1" appended. Likewise, the state variables at time t_2 will have a "_2" appended.

3. The specification of dprocess using an R function is sometimes useful for pedagogical or debugging purposes, though one ultimately wants to move to C snippets because of their much greater computational speed. As with all other basic model components in pomp versions >=2, one writes a function with arguments chosen from among the state variables, covariates, and parameters. The only wrinkle special to dprocess is that the state variables will have _1 or _2 appended to their names according to whether they correspond to times t_1 or t_2. Again, the ?dprocess_spec page has all this spelled out.

Because all of this may seem a bit abstract, here's an example for the simple case where we have a Brownian motion as our latent state process.

R function approach

pomp(
     ...,
     dprocess = function (t_1, t_2, x_1, x_2, sigma, ...) {
         standdev <- sigma * sqrt(t_2 - t_1)
     loglik <- dnorm(x=x_2, mean=x_1, sd=standdev, log=TRUE)
         loglik          
     },
     ...
    )

C snippet approach

pomp(
     ...,
     dprocess = Csnippet("
         loglik = dnorm(x_2,x_1,sigma*sqrt(t_2-t_1),1);
     "),
      ...,
      statenames=c("x"),
      paramnames=c("sigma"),
      ...
     )

I hope this helps. If it does not, or if more is needed, please don't hesitate to reopen the issue. I hope to hear more about your experiences!

[sdwfrost]

Thanks! I didn't know about dprocess_spec - the _1/_2 magic is precisely what I was looking for.