eulermultinom {pomp} | R Documentation |
Eulermultinomial and gamma-whitenoise distributions
Description
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.
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. |
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. |
Details
If individuals face constant hazards of death in
ways
at rates
,
then in an interval of duration
,
the number of individuals remaining alive and dying in each way is multinomially distributed:
where is the number of individuals remaining alive and
is the number of individuals dying in way
over the interval.
Here, the probability of remaining alive is
and the probability of dying in way is
In this case, we say that
where .
Draw
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 are the numbers of individuals who have died in each of the
ways over the interval
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 with
,
where
is the increment of an integrated Gamma white noise process with intensity
.
That is,
has mean
and variance
.
The resulting process is overdispersed and converges (as
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.
Value
reulermultinom |
Returns a |
deulermultinom |
Returns a vector (of length equal to the number of columns of |
rgammawn |
Returns a vector of length |
C API
An interface for C codes using these functions is provided by the package. Visit the package homepage to view the pomp C API document.
Author(s)
Aaron A. King
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. 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. 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
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))