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)
rgammawn(n = 1, sigma, dt)
integer; number of random variates to generate.
scalar integer; number of individuals at risk.
numeric vector of hazard rates.
numeric scalar; duration of Euler step.
matrix or vector containing number of individuals that have succumbed to each death process.
logical; if TRUE, return logarithm(s) of probabilities.
numeric scalar; intensity of the Gamma white noise process.
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.
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).
Returns a vector of length n containing random increments of the integrated Gamma white noise process with intensity sigma.
An interface for C codes using these functions is provided by the package. Visit the package homepage to view the pomp C API document.
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 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
dn <- 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).
Breto & 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 {\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 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.
For all of the functions described here, access to the underlying C routines is available: see below.
C. and E. L. Ionides. Compound Markov counting processe and their applications to modeling infinitesimally over-dispersed systems. Stochastic Processes and their Applications 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. Journal of the Royal Society Interface 7, 271--283, 2010.
More on implementing POMP models:
Csnippet,
accumulators,
basic_components,
covariate_table(),
dmeasure_spec,
dprocess_spec,
parameter_trans(),
pomp-package,
prior_spec,
rinit_spec,
rmeasure_spec,
rprocess_spec,
skeleton_spec,
transformations,
userdata
# NOT RUN {
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))
# }
Run the code above in your browser using DataLab