How to use ddgamma(x,shape,scale) in rprocess of Csnippets?

[SwethaLal]

I was trying to use the discrete gamma distribution in an rprocess step and it showed me the error:

**In function ‘__pomp_stepfn’:
/tmp/Rtmp0oaYOd/114435/pomp_de2c59a6b36029764925da605840af96.c:73:6: warning: implicit declaration of function ‘ddgamma’; did you mean ‘dgamma’? [-Wimpl**

It looks like gamma distribution is defined and can be used here, whereas discrete gamma distribution is not to be found.

The R code I wrote contained the following. mu is the mean and var is the variance of the discrete gamma distribution and k and r gives the shape and rate respectively for the corresponding mu and var values. mu and var are parameters in the code. I've also used library(extraDistr) for enabling the use of the discrete gamma distribution commands.

Bstep <- Csnippet("
double l=0.0,v, k, r;
k = ((mu*mu)/var); 
r = (mu/var);
double *x = &X1;
int i, n = (int) N;
for (i = n-1; i > 0; i--) {
v = ddgamma(i, shape=k, rate=r);
l = l+gv;
x[i] = x[i-1];
}
v = ddgamma(0, shape=k, rate=r);
l = l+v;
x[0] = rpois(l);
")

Binit <- Csnippet("
double *x = &X1;
int i, n = (int) N;
for (i=1; i < n; i++) x[i] = 0.0;
x[0] = 1;
")

dmeas <- Csnippet("
lik = dbinom(X,X1,rho,give_log);
")

rmeas <- Csnippet("
X = rbinom(X1,rho);
")

paramets <- c(N=7,R0=4,mu=5.5,var=5.5,rho=1)

dat %>%
  pomp(
    times="Day",t0=0,
    rprocess=euler(Bstep,delta.t=1),
    rinit=Binit,
    rmeasure=rmeas,
    dmeasure=dmeas,
    partrans=parameter_trans(
      log=c("R0"),
      ),
    statenames=sprintf("X%d",1:7),
    paramnames=c("N","R0","mu","var","rho"),
    params=paramets
  ) -> Bsim

Bsim %>%
  simulate(params=paramets,
    nsim=10,
    format="data.frame",
    include.data=TRUE

  ) -> sims

This gave the error mentioned above. Can we somehow add or call or install the discrete gamma distribution to pomp such that it can be used inside the Csnippets? It works fine otherwise in R, but I was wondering if it could somehow be included for usage inside Csnippets.

Thank you!:)

My R-studio version and system version is as below.

-----------------------------------
 sessionInfo 
-----------------------------------
R version 4.0.3 (2020-10-10)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.2 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0

locale:
 [1] LC_CTYPE=en_GB.UTF-8       LC_NUMERIC=C               LC_TIME=en_GB.UTF-8        LC_COLLATE=en_GB.UTF-8     LC_MONETARY=en_GB.UTF-8    LC_MESSAGES=en_GB.UTF-8   
 [7] LC_PAPER=en_GB.UTF-8       LC_NAME=C                  LC_ADDRESS=C               LC_TELEPHONE=C             LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] parallel  stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] doRNG_1.8.2       rngtools_1.5      doParallel_1.0.16 iterators_1.0.13  foreach_1.5.1     extraDistr_1.9.1  pomp_3.2          forcats_0.5.1     stringr_1.4.0    
[10] dplyr_1.0.5       purrr_0.3.4       readr_1.4.0       tidyr_1.1.3       tibble_3.1.0      ggplot2_3.3.3     tidyverse_1.3.0  

loaded via a namespace (and not attached):
 [1] deSolve_1.28      tidyselect_1.1.0  reshape2_1.4.4    haven_2.3.1       lattice_0.20-41   colorspace_2.0-0  vctrs_0.3.6       generics_0.1.0    viridisLite_0.3.0
[10] rlang_0.4.10      pillar_1.5.1      glue_1.4.2        withr_2.4.1       DBI_1.1.1         dbplyr_2.1.0      modelr_0.1.8      readxl_1.3.1      lifecycle_1.0.0  
[19] plyr_1.8.6        munsell_0.5.0     gtable_0.3.0      cellranger_1.1.0  rvest_1.0.0       mvtnorm_1.1-1     codetools_0.2-18  coda_0.19-4       labeling_0.4.2   
[28] broom_0.7.5       Rcpp_1.0.6        scales_1.1.1      backports_1.2.1   jsonlite_1.7.2    farver_2.1.0      fs_1.5.0          hms_1.0.0         digest_0.6.27    
[37] stringi_1.5.3     grid_4.0.3        cli_2.3.1         tools_4.0.3       magrittr_2.0.1    crayon_1.4.1      pkgconfig_2.0.3   ellipsis_0.3.1    xml2_1.3.2       
[46] reprex_1.0.0      lubridate_1.7.10  assertthat_0.2.1  httr_1.4.2        rstudioapi_0.13   R6_2.5.0          compiler_4.0.3 

[kingaa]

With a C snippet, you can access any C function, provided you can link to it. By default, pomp gives you access to the R C API, which includes many distribution functions, but not all. In particular, the discrete Gamma distribution is not provided. I do not know whether any add-on packages provide a C interface for this function.

I was not familiar with this distribution, but I read (in the documentation of the ddgamma function in extraDistr, that this distribution has p.d.f. f(x) = F(x+1)-F(x), where F is the c.d.f. of the Gamma distribution. Is this right? If so, you could simply do

v = pgamma(I+1,k,1/r) - pgamma(I,k,1/r);

This is assuming that k is the shape parameter and 1/r is the scale parameter. Note that this uses the C API for R which provides pgamma (the c.d.f. of the Gamma distribution) using the shape/scale (not shape/rate) parameterization.

Note also that C does not support named arguments in function calls.