pomp (version 0.39-3)

eulermultinom: Euler-multinomial death process

Description

Density and random-deviate generation for the Euler-multinomial death process with parameters size, rate, and dt.

Usage

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

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.

Value

• reulermultinom
  Returns a length(rate) by n matrix. Each column is a different random draw. Each row contains the numbers of individuals succumbed to the corresponding process.
• deulermultinom
  Returns a vector (of length equal to the number of columns of x) containing the probabilities of observing each column of x given the specified parameters (size, rate, dt).

Details

If $N$ individuals face constant hazards of death in $k$ ways at rates $r_1, r_2, \dots, r_k$, then in an interval of duration $\Delta t$, 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 $\Delta n_i$ is the number of individuals dying in way $i$ over the interval, the probability of remaining alive is $p_0=\exp(-\sum_i r_i \Delta t)$, 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 can say that $$(\Delta n_1, \dots, \Delta n_k) \sim \mathrm{eulermultinom}(N,r,\Delta t),$$ where $r=(r_1,\dots,r_k)$. Draw $m$ random samples from this distribution by doing reulermultinom(n=m,size=N,rate=r,dt=dt),

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

deulermultinom(x=x,size=N,rate=r,dt=dt). Direct access to the underlying C routines is available: see the header file pomp.h, included with the package.

Examples

print(x <- reulermultinom(5,size=100,rate=c(a=1,b=2,c=3),dt=0.1))
deulermultinom(x,size=100,rate=c(1,2,3),dt=0.1)