accumvars {pomp} | R Documentation |
accumulator variables
Description
Latent state variables that accumulate quantities through time.
Details
In formulating models, one sometimes wishes to define a state variable that will accumulate some quantity over the interval between successive observations.
pomp provides a facility to make such features more convenient.
Specifically, if is a state-variable named in the
pomp
's accumvars
argument, then for each interval ,
,
will be set to zero at prior to any
rprocess
computation over that interval.
Deterministic trajectory computation is handled slightly differently:
see flow
.
For examples, see sir
and the tutorials on the package website.
See Also
More on implementing POMP models:
Csnippet
,
basic_components
,
betabinomial
,
covariates
,
dinit_spec
,
dmeasure_spec
,
dprocess_spec
,
emeasure_spec
,
eulermultinom
,
parameter_trans()
,
pomp-package
,
pomp_constructor
,
prior_spec
,
rinit_spec
,
rmeasure_spec
,
rprocess_spec
,
skeleton_spec
,
transformations
,
userdata
,
vmeasure_spec
Examples
## A simple SIR model.
ewmeas |>
subset(time < 1952) |>
pomp(
times="time",t0=1948,
rprocess=euler(
Csnippet("
int nrate = 6;
double rate[nrate]; // transition rates
double trans[nrate]; // transition numbers
double dW;
// gamma noise, mean=dt, variance=(sigma^2 dt)
dW = rgammawn(sigma,dt);
// compute the transition rates
rate[0] = mu*pop; // birth into susceptible class
rate[1] = (iota+Beta*I*dW/dt)/pop; // force of infection
rate[2] = mu; // death from susceptible class
rate[3] = gamma; // recovery
rate[4] = mu; // death from infectious class
rate[5] = mu; // death from recovered class
// compute the transition numbers
trans[0] = rpois(rate[0]*dt); // births are Poisson
reulermultinom(2,S,&rate[1],dt,&trans[1]);
reulermultinom(2,I,&rate[3],dt,&trans[3]);
reulermultinom(1,R,&rate[5],dt,&trans[5]);
// balance the equations
S += trans[0]-trans[1]-trans[2];
I += trans[1]-trans[3]-trans[4];
R += trans[3]-trans[5];
"),
delta.t=1/52/20
),
rinit=Csnippet("
double m = pop/(S_0+I_0+R_0);
S = nearbyint(m*S_0);
I = nearbyint(m*I_0);
R = nearbyint(m*R_0);
"),
paramnames=c("mu","pop","iota","gamma","Beta","sigma",
"S_0","I_0","R_0"),
statenames=c("S","I","R"),
params=c(mu=1/50,iota=10,pop=50e6,gamma=26,Beta=400,sigma=0.1,
S_0=0.07,I_0=0.001,R_0=0.93)
) -> ew1
ew1 |>
simulate() |>
plot(variables=c("S","I","R"))
## A simple SIR model that tracks cumulative incidence.
ew1 |>
pomp(
rprocess=euler(
Csnippet("
int nrate = 6;
double rate[nrate]; // transition rates
double trans[nrate]; // transition numbers
double dW;
// gamma noise, mean=dt, variance=(sigma^2 dt)
dW = rgammawn(sigma,dt);
// compute the transition rates
rate[0] = mu*pop; // birth into susceptible class
rate[1] = (iota+Beta*I*dW/dt)/pop; // force of infection
rate[2] = mu; // death from susceptible class
rate[3] = gamma; // recovery
rate[4] = mu; // death from infectious class
rate[5] = mu; // death from recovered class
// compute the transition numbers
trans[0] = rpois(rate[0]*dt); // births are Poisson
reulermultinom(2,S,&rate[1],dt,&trans[1]);
reulermultinom(2,I,&rate[3],dt,&trans[3]);
reulermultinom(1,R,&rate[5],dt,&trans[5]);
// balance the equations
S += trans[0]-trans[1]-trans[2];
I += trans[1]-trans[3]-trans[4];
R += trans[3]-trans[5];
H += trans[3]; // cumulative incidence
"),
delta.t=1/52/20
),
rmeasure=Csnippet("
double mean = H*rho;
double size = 1/tau;
reports = rnbinom_mu(size,mean);
"),
rinit=Csnippet("
double m = pop/(S_0+I_0+R_0);
S = nearbyint(m*S_0);
I = nearbyint(m*I_0);
R = nearbyint(m*R_0);
H = 0;
"),
paramnames=c("mu","pop","iota","gamma","Beta","sigma","tau","rho",
"S_0","I_0","R_0"),
statenames=c("S","I","R","H"),
params=c(mu=1/50,iota=10,pop=50e6,gamma=26,
Beta=400,sigma=0.1,tau=0.001,rho=0.6,
S_0=0.07,I_0=0.001,R_0=0.93)
) -> ew2
ew2 |>
simulate() |>
plot()
## A simple SIR model that tracks weekly incidence.
ew2 |>
pomp(accumvars="H") -> ew3
ew3 |>
simulate() |>
plot()