kingaa / pomp

R package for statistical inference using partially observed Markov processes
https://kingaa.github.io/pomp
GNU General Public License v3.0
110 stars 26 forks source link

Can I access the files related to Csnippet and directly edit? #186

Closed Fuhan-Yang closed 1 year ago

Fuhan-Yang commented 1 year ago

Thank you for developing such a fantastic package. I'm new to pomp, and my question is how to locate all the files related to Csnippet (source code, header files, etc)? Besides, can I directly edit the source code or is it recommended to do so in Csnippet in R?

My specific question is that, is there an easier way to embed a long external data stream in Csnippet of rprocess? If I want to add a function to link the data to the source code, can this be done in Csnippet in R, or does it need to be done in the source code in C?

kingaa commented 1 year ago

Use spy to view the C file you've created with your C snippets. By default, this is stored in the temporary directory. You can change this using the cdir and cfile options.

Most people find it easiest to work with C snippets, but you can certainly link to a library of your own. See FAQ 3.7. The C API for pomp is described here. In particular, if you furnish your own basic model components, they must conform to the prototypes for such functions.

I am not sure what you have in mind by "a long external data stream". Are you envisioning that rprocess should depend on covariates? If so, does the existing facility for covariates solve the problem?

More generally, you can do almost anything you want in the C snippets, as long as they are valid C. That doesn't mean that the resulting calculation will be correct, or even meaningful, of course!

Fuhan-Yang commented 1 year ago

Thanks for the answers. Does that mean the declarations void pomp_* should always stay the same?

The covariates facility is what I need. But I'm having issues running simulations. Below is my test pomp model:

#### daily integrated pomp model ####

library(pomp)
library(dplyr)
library(ggplot2)

dat <- data.frame(
  daynum = c(1:3652),
  dat = rbinom(3652,size = 10, prob = 0.3) + c(1:3652)* 0.1 - 0.3
)

covar <- data.frame(
  daynum = dat$daynum,
  step = 1.1
)

### initial condition ###
step_rinit <- Csnippet("
  x = 0;
  y = 0;
")

### process model ###

## stochastic version ##
step_rproc <- Csnippet("

    double dx = 0.1;
    double dy = -0.1;

    dx = dx + step * s;
    dy = dy + 0.5 * step * s;

    y = y + dx + dy;
    x = x + dx + dy;

")

### measure model ###
step_dmeas <- Csnippet("
  lik = dbinom(dat, y, rho, give_log);
")

step_rmeas <- Csnippet("
  dat = rbinom(y, rho);
")

##### plug in pomp ####
test_pomp <- 
  pomp(data = dat,
       times = 'daynum',
       t0 = 1,
       rinit = step_rinit,
       rprocess = euler(step_rproc,delta.t = 1),
       rmeasure = step_rmeas,
       dmeasure = step_dmeas,
       covar = covariate_table(covar$step,times = covar$daynum),
       covarnames = 'step',
       obsnames = 'dat',
       statenames = c('x','y'),
       paramnames = c('s','rho'),
       cdir = 'c')

test_pomp |>
  simulate(nsim = 20, 
           params = c(s = 0.2,rho = 1),
           covar = covariate_table(covar$step,times = covar$daynum),
           covarnames = 'step',
           format = 'data.frame') |>
  ggplot() + 
  geom_line(aes(x = daynum, y = y,color = .id))

The model just tests if the process model can be affected by the covariates. There are no issues building the pomp model, but error appeared in simulate as "Error: in ‘simulate’: variable 'step' not found among the covariates." . Can you please tell me what happened?

kingaa commented 1 year ago

Cf. the following, in particular, the call to covariate_table. Note there is no need to reinsert covar in the call to simulate.

##### plug in pomp ####
dat |>
  pomp(
    times = "daynum",
    t0 = 1,
    rinit = step_rinit,
    rprocess = euler(step_rproc,delta.t = 1),
    rmeasure = step_rmeas,
    dmeasure = step_dmeas,
    covar = covariate_table(covar,times="daynum"),
    statenames = c("x","y"),
    paramnames = c("s","rho"),
    cdir = "c"
  ) -> test_pomp

test_pomp |>
  simulate(
    nsim = 20, 
    params = c(s = 0.2,rho = 1),
    format = "data.frame"
  ) |>
  ggplot(aes(x = daynum, y = y, color = .id)) + 
  geom_line()

While I think the above addresses your immediate question, and I take it that you've cobbled this up for purely illustrative purposes, I do point out that it makes little sense to include step as a covariate, since it is constant and could more properly be included as a parameter.

Fuhan-Yang commented 1 year ago

Thank you so much for the prompt answer! Yes, this model is oversimplified to test the code implementation and to visually check the effect of the simple covariate. Thanks again for all the guidance. My issue is resolved and I will close it now.