pomp (version 1.10)

eulermultinom: The Euler-multinomial distributions and Gamma white-noise processes

Description

This page documents both the Euler-multinomial family of distributions and the package's simulator of Gamma white-noise processes.

Usage

reulermultinom(n = 1, size, rate, dt)
deulermultinom(x, size, rate, dt, log = FALSE)
rgammawn(n = 1, sigma, dt)

Arguments

n

integer; number of random variates to generate.

size

scalar integer; number of individuals at risk.

rate

numeric vector of hazard rates.

sigma

numeric scalar; intensity of the Gamma white noise process.

dt

numeric scalar; duration of Euler step.

x

matrix or vector containing number of individuals that have succumbed to each death process.

log

logical; if TRUE, return logarithm(s) of probabilities.

Value

C API

An interface for C codes using these functions is provided by the package. At an R prompt, execute

file.show(system.file("include/pomp.h",package="pomp"))

to view the ‘pomp.h’ header file that defines and explains the API.

Details

If $N$ individuals face constant hazards of death in $k$ ways at rates $r1,r2,\dots,rk$, then in an interval of duration $dt$, the number of individuals remaining alive and dying in each way is multinomially distributed: $$(N-\sum_{i=1}^k \Delta n_i, \Delta n_1, \dots, \Delta n_k) \sim \mathrm{multinomial}(N;p_0,p_1,\dots,p_k),$$ where $dni$ is the number of individuals dying in way $i$ over the interval, the probability of remaining alive is $p0=exp(-\sum(ri dt))$, and the probability of dying in way $j$ is $$p_j=\frac{r_j}{\sum_i r_i} (1-\exp(-\sum_i r_i \Delta t)).$$ In this case, we say that $$(\Delta n_1, \dots, \Delta n_k) \sim \mathrm{eulermultinom}(N,r,\Delta t),$$ where $r=(r1,\dots,rk)$. Draw $m$ random samples from this distribution by doing

dn <- reulermultinom(n=m,size=N,rate=r,dt=dt),

where r is the vector of rates. Evaluate the probability that $x=(x1,\dots,xk)$ are the numbers of individuals who have died in each of the $k$ ways over the interval $$dt, by doing

deulermultinom(x=x,size=N,rate=r,dt=dt).

Breto & Ionides (2011) discuss how an infinitesimally overdispersed death process can be constructed by compounding a binomial process with a Gamma white noise process. The Euler approximation of the resulting process can be obtained as follows. Let the increments of the equidispersed process be given by

reulermultinom(size=N,rate=r,dt=dt).

In this expression, replace the rate $r$ with $r {\Delta W}/{\Delta t}$, where $\Delta\!W \sim \mathrm{Gamma}(\Delta\!t/\sigma^2,\sigma^2)$ is the increment of an integrated Gamma white noise process with intensity $\sigma$. That is, $\Delta\!W$ has mean $\Delta\!t$ and variance $\sigma^2 \Delta\!t$. The resulting process is overdispersed and converges (as $\Delta t$ goes to zero) to a well-defined process. The following lines of R code accomplish this:

dW <- rgammawn(sigma=sigma,dt=dt)

dn <- reulermultinom(size=N,rate=r,dt=dW)

or

dn <- reulermultinom(size=N,rate=r*dW/dt,dt=dt).

He et al. use such overdispersed death processes in modeling measles. For all of the functions described here, access to the underlying C routines is available: see below.

References

C. Breto & E. L. Ionides, Compound Markov counting processes and their applications to modeling infinitesimally over-dispersed systems. Stoch. Proc. Appl., 121:2571--2591, 2011.

D. He, E. L. Ionides, & A. A. King, Plug-and-play inference for disease dynamics: measles in large and small populations as a case study. J. R. Soc. Interface, 7:271--283, 2010.

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))