eulermultinom function

Eulermultinomial and Gamma-whitenoise distributions

pomp provides a number of probability distributions that have proved useful in modeling partially observed Markov processes. These include the Euler-multinomial family of distributions and the the Gamma white-noise processes.

reulermultinom(n = 1, size, rate, dt)

deulermultinom(x, size, rate, dt, log = FALSE)

eeulermultinom(size, rate, dt)

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.
• 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.
• sigma: numeric scalar; intensity of the Gamma white noise process.

Returns

reulermultinom returns a length(rate) by n matrix. Each column is a different random draw. Each row contains the numbers of individuals that have succumbed to the corresponding process.

deulermultinom returns a vector (of length equal to the number of columns of x). This contains the probabilities of observing each column of x given the specified parameters (size, rate, dt).

eeulermultinom returns a length(rate)-vector containing the expected number of individuals to have succumbed to the corresponding process. With the notation above, if

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

then the $k$-th element of the vector dn will be $pk N$.

rgammawn returns a vector of length n containing random increments of the integrated Gamma white noise process with intensity sigma.

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:

$$(\Delta{n_0}, \Delta{n_1}, \dots, \Delta{n_K}) \sim \mathrm{Multinomial}(N;p_0,p_1,\dots,p_K),(dn0, dn1, \dots, dnK) ~ multinomial(N; p0, p1, \dots, pK),$$

where $dn0 = N - sum(dnk,k=1:K)$ is the number of individuals remaining alive and $dnk$ is the number of individuals dying in way $k$ over the interval. Here, the probability of remaining alive is

$$p_0=\exp(-\sum_k r_k \Delta{t})p0 = exp(-sum(rk, k=1:K)*dt)$$

and the probability of dying in way $k$ is

$$p_k=\frac{r_k}{\sum_j r_j} (1-p_0).pk = (1-p0)*rk/sum(rj, j=1:K).$$

In this case, we say that

$$(\Delta{n_1},\dots,\Delta{n_K})\sim\mathrm{Eulermultinom}(N,r,\Delta t),(dn1,\dots,dnK) ~ eulermultinom(N,r,dt),$$

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=$dt, by doing

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

Bretó & Ionides (2011) discuss how an infinitesimally overdispersed death process can be constructed by compounding a multinomial 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*dW/dt$, where $dW ~ Gamma(dt/sigma^2,sigma^2)$

is the increment of an integrated Gamma white noise process with intensity $sigma$. That is, $dW$ has mean $dt$ and variance $sigma^2*dt$. The resulting process is overdispersed and converges (as $dt$ goes to zero) to a well-defined process. The following lines of 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. (2010) use such overdispersed death processes in modeling measles and the "Simulation-based Inference" course discusses the value of allowing for overdispersion more generally.

For all of the functions described here, access to the underlying C routines is available: see below.

C API

An interface for C codes using these functions is provided by the package. Visit the package homepage to view the c("list("pomp")", " C API document").

Examples

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

References

C. Bretó and E. L. Ionides. Compound Markov counting processes and their applications to modeling infinitesimally over-dispersed systems. Stochastic Processes and their Applications 121 , 2571--2591, 2011. tools:::Rd_expr_doi("10.1016/j.spa.2011.07.005") . D. He, E.L. Ionides, and A.A. King. Plug-and-play inference for disease dynamics: measles in large and small populations as a case study. Journal of the Royal Society Interface 7 , 271--283, 2010. tools:::Rd_expr_doi("10.1098/rsif.2009.0151") .

See Also

More on implementing POMP models: Csnippet, accumvars, basic_components, betabinomial, covariates, dinit_spec, dmeasure_spec, dprocess_spec, emeasure_spec, parameter_trans(), pomp-package, pomp_constructor, prior_spec, rinit_spec, rmeasure_spec, rprocess_spec, skeleton_spec, transformations, userdata, vmeasure_spec

Author(s)

Aaron A. King