eeulermultinom at R level

DESCRIPTION

@@ -1,8 +1,8 @@
 Package: pomp
 Type: Package
 Title: Statistical Inference for Partially Observed Markov Processes
-Version: 6.0.2.3
-Date: 2024-12-18
+Version: 6.0.3.0
+Date: 2024-12-19
 Authors@R: c(person(given=c("Aaron","A."),family="King",role=c("aut","cre"),email="[email protected]",comment=c(ORCID="0000-0001-6159-3207")),
 	  person(given=c("Edward","L."),family="Ionides",role="aut",comment=c(ORCID="0000-0002-4190-0174"))	,
 	  person(given="Carles",family="Bretó",role="aut",comment=c(ORCID="0000-0003-4695-4902")),

NAMESPACE

@@ -26,6 +26,7 @@ export(dbetabinom)
 export(deulermultinom)
 export(discrete_time)
 export(ebolaModel)
+export(eeulermultinom)
 export(euler)
 export(expit)
 export(freeze)

R/eulermultinom.R

@@ -46,7 +46,6 @@
 ##'
 ##' For all of the functions described here, access to the underlying C routines is available:
 ##' see below.
-##'
 ##' @name eulermultinom
 ##' @rdname eulermultinom
 ##' @family implementation information
@@ -59,45 +58,22 @@
 ##' @param x matrix or vector containing number of individuals that have
 ##' succumbed to each death process.
 ##' @param log logical; if TRUE, return logarithm(s) of probabilities.
-##' @return
-##' \item{reulermultinom}{
-##'    Returns a \code{length(rate)} by \code{n} matrix.
-##'    Each column is a different random draw.
-##'    Each row contains the numbers of individuals that have succumbed to the corresponding process.
-##' }
-##' \item{deulermultinom}{
-##'    Returns a vector (of length equal to the number of columns of \code{x}) containing
-##'    the probabilities of observing each column of \code{x} given the specified parameters (\code{size}, \code{rate}, \code{dt}).
-##' }
-##' \item{rgammawn}{
-##'    Returns a vector of length \code{n} containing random increments of the integrated Gamma white noise process with intensity \code{sigma}.
-##' }
-##'
 ##' @section C API:
 ##' An interface for C codes using these functions is provided by the package.
 ##' Visit the package homepage to view the \href{https://kingaa.github.io/pomp/C_API.html}{\pkg{pomp} C API document}.
-##'
 ##' @author Aaron A. King
-##'
 ##' @references
-##'
 ##' \Breto2011
-##'
 ##' \He2010
-##'
 ##' @keywords distribution
-##' @examples
-##'
-##' print(dn <- reulermultinom(5,size=100,rate=c(a=1,b=2,c=3),dt=0.1))
-##' deulermultinom(x=dn,size=100,rate=c(1,2,3),dt=0.1)
-##' ## an Euler-multinomial with overdispersed transitions:
-##' dt <- 0.1
-##' dW <- rgammawn(sigma=0.1,dt=dt)
-##' print(dn <- reulermultinom(5,size=100,rate=c(a=1,b=2,c=3),dt=dW))
-##'
+##' @example examples/eulermultinom.R
 NULL
 
 ##' @rdname eulermultinom
+##' @return
+##' \code{reulermultinom} returns a \code{length(rate)} by \code{n} matrix.
+##' Each column is a different random draw.
+##' Each row contains the numbers of individuals that have succumbed to the corresponding process.
 ##' @export
 reulermultinom <- function (n = 1, size, rate, dt) {
   tryCatch(
@@ -107,6 +83,9 @@ reulermultinom <- function (n = 1, size, rate, dt) {
 }
 
 ##' @rdname eulermultinom
+##' @return
+##' \code{deulermultinom} returns a vector (of length equal to the number of columns of \code{x}).
+##' This contains the probabilities of observing each column of \code{x} given the specified parameters (\code{size}, \code{rate}, \code{dt}).
 ##' @export
 deulermultinom <- function (x, size, rate, dt, log = FALSE) {
   tryCatch(
@@ -116,6 +95,20 @@ deulermultinom <- function (x, size, rate, dt, log = FALSE) {
 }
 
 ##' @rdname eulermultinom
+##' @return
+##' \code{eeulermultinom} returns a \code{length(rate)}-vector
+##' containing the expected number of individuals to have succumbed to the corresponding process.
+##' @export
+eeulermultinom <- function (size, rate, dt) {
+  tryCatch(
+    .Call(P_E_Euler_Multinom,size,rate,dt),
+    error = function (e) pStop(who="eeulermultinom",conditionMessage(e))
+  )
+}
+
+##' @rdname eulermultinom
+##' @return
+##' \code{rgammawn} returns a vector of length \code{n} containing random increments of the integrated Gamma white noise process with intensity \code{sigma}.
 ##' @export
 rgammawn <- function (n = 1, sigma, dt) {
   tryCatch(

TODO.md

@@ -5,8 +5,7 @@
 ## For pomp:
 
 - Rosenzweig-MacArthur example
-- Make certain error checks in inline codes depend on NDEBUG flag
-- **R** interface to `eeulermultinom`
+- Make certain error checks in inline codes depend on NDEBUG flag?
 - replace `melt` with `data.table::melt`?
 - YAML interface
 - "Getting Started" vignette:

examples/eulermultinom.R

@@ -0,0 +1,18 @@
+## Simulate 5 realizations of Euler-multinomial random variable:
+
+dn <- reulermultinom(5,size=100,rate=c(a=1,b=2,c=3),dt=0.1)
+dn
+
+## Compute the probability mass function at each of the 5 realizations:
+
+deulermultinom(x=dn,size=100,rate=c(1,2,3),dt=0.1)
+
+## Compute the expectation of an Euler-multinomial:
+
+eeulermultinom(size=100,rate=c(a=1,b=2,c=3),dt=0.1)
+
+## An Euler-multinomial with overdispersed transitions:
+
+dt <- 0.1
+dW <- rgammawn(sigma=0.1,dt=dt)
+reulermultinom(5,size=100,rate=c(a=1,b=2,c=3),dt=dW)

examples/spy.R

@@ -1,5 +1,9 @@
-ricker() |> spy()
+\donttest{
 
-sir() |> spy()
+  ricker() |> spy()
 
-sir2() |> spy()
+  sir() |> spy()
+
+  sir2() |> spy()
+
+}

src/decls.h

@@ -8,6 +8,7 @@ extern SEXP do_dinit(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP);
 /* src/distributions.c */
 extern SEXP R_Euler_Multinom(SEXP, SEXP, SEXP, SEXP);
 extern SEXP D_Euler_Multinom(SEXP, SEXP, SEXP, SEXP, SEXP);
+extern SEXP E_Euler_Multinom(SEXP, SEXP, SEXP);
 extern SEXP R_GammaWN(SEXP, SEXP, SEXP);
 extern SEXP R_BetaBinom(SEXP, SEXP, SEXP, SEXP);
 extern SEXP D_BetaBinom(SEXP, SEXP, SEXP, SEXP, SEXP);

src/distributions.c

@@ -68,6 +68,24 @@ SEXP D_Euler_Multinom (SEXP x, SEXP size, SEXP rate, SEXP deltat, SEXP log) {
   return f;
 }
 
+SEXP E_Euler_Multinom (SEXP size, SEXP rate, SEXP deltat) {
+  int ntrans = length(rate);
+  SEXP x, nm;
+  if (length(size)>1)
+    warn("in 'eeulermultinom': only the first element of 'size' is meaningful");
+  if (length(deltat)>1)
+    warn("in 'eeulermultinom': only the first element of 'dt' is meaningful");
+  PROTECT(size = AS_NUMERIC(size));
+  PROTECT(rate = AS_NUMERIC(rate));
+  PROTECT(deltat = AS_NUMERIC(deltat));
+  PROTECT(x = NEW_NUMERIC(ntrans));
+  PROTECT(nm = GET_NAMES(rate));
+  SET_NAMES(x,nm);
+  eeulermultinom(ntrans,*REAL(size),REAL(rate),*REAL(deltat),REAL(x));
+  UNPROTECT(5);
+  return x;
+}
+
 // This function draws a random increment of a gamma whitenoise process.
 // This will have expectation=dt and variance=(sigma^2*dt)
 // If dW = rgammawn(sigma,dt), then

src/init.c

@@ -12,6 +12,7 @@ static const R_CallMethodDef callMethods[] = {
   {"load_stack_decr", (DL_FUNC) &load_stack_decr, 1},
   {"R_Euler_Multinom", (DL_FUNC) &R_Euler_Multinom, 4},
   {"D_Euler_Multinom", (DL_FUNC) &D_Euler_Multinom, 5},
+  {"E_Euler_Multinom", (DL_FUNC) &E_Euler_Multinom, 3},
   {"R_GammaWN", (DL_FUNC) &R_GammaWN, 3},
   {"R_BetaBinom", (DL_FUNC) &R_BetaBinom, 4},
   {"D_BetaBinom", (DL_FUNC) &D_BetaBinom, 5},

tests/eulermultinom.R

@@ -8,6 +8,7 @@ reulermultinom(n=5,size=100,rate=c(1,-2,3),dt=0.1)
 reulermultinom(n=5,size=100,rate=c(1,NA,3),dt=0.1)
 reulermultinom(n=5,size=100.3,rate=c(1,2,3),dt=0.1)
 reulermultinom(n=0,size=100,rate=c(1,2,3),dt=0.1)
+reulermultinom(n=5,size=100,rate=c(1,2,3),dt=Inf)
 try(reulermultinom(n=-2,size=100,rate=c(1,2,3),dt=0.1))
 reulermultinom(n=1,size=100,rate=c(1,2e400,3),dt=0.1)
 reulermultinom(n=1,size=100,rate=c(1,2,3),dt=c(0.1,0.2,0.3,Inf))
@@ -31,6 +32,25 @@ deulermultinom(x=c(3,4,0),size=10,rate=c(1,1,-1),dt=0.1)
 deulermultinom(x=c(3,-4,0),size=10,rate=c(1,1,1),dt=0.1)
 deulermultinom(x=c(3,6,3),size=10,rate=c(1,1,0),dt=0.1,log=TRUE)
 
+eeulermultinom(size=100,rate=c(1,2,3),dt=0.1)
+eeulermultinom(size=-3,rate=c(1,2,3),dt=0.1)
+eeulermultinom(size=100,rate=c(1,-2,3),dt=0.1)
+eeulermultinom(size=100,rate=c(1,NA,3),dt=0.1)
+eeulermultinom(size=100.3,rate=c(1,2,3),dt=0.1)
+eeulermultinom(size=100,rate=c(1,2,3),dt=0.1)
+eeulermultinom(size=100,rate=c(1,2,3),dt=Inf)
+eeulermultinom(size=100,rate=c(1,2e400,3),dt=0.1)
+eeulermultinom(size=100,rate=c(1,2,3),dt=c(0.1,0.2,0.3,Inf))
+eeulermultinom(size=100,rate=c(1,2,3),dt=c(0.1,0.2,0.3,0.5))
+eeulermultinom(size=c(100,200,300),rate=c(1,2,3),dt=0.2)
+eeulermultinom(size=0,rate=c(1,2,3),dt=0.2)
+eeulermultinom(size=10,rate=c(1,Inf,1),0.1)
+eeulermultinom(size=Inf,rate=c(1,100,1),0.1)
+eeulermultinom(size=NA,rate=c(1,2,3),dt=1)
+eeulermultinom(size=100,rate=c(1,NA,3),dt=1)
+eeulermultinom(size=100,rate=c(1,2,3),dt=NA)
+eeulermultinom(size=100,rate=c(0,0,0,0),dt=1)
+
 rgammawn(n=5,sigma=2,dt=0.1)
 rgammawn(n=5,sigma=1:5,dt=0.1)
 rgammawn(n=5,sigma=1,dt=rep(1,5))