pomp (version 3.3)

accumulators: accumulators

Description

Accumulator variables

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, variables named in the pomp's accumvars argument will be set to zero immediately following each observation. See sir and the tutorials on the package website for examples.

See Also

sir

More on implementing POMP models: Csnippet, basic_components, covariate_table(), distributions, dmeasure_spec, dprocess_spec, parameter_trans(), pomp-package, prior_spec, rinit_spec, rmeasure_spec, rprocess_spec, skeleton_spec, transformations, userdata

Examples

# NOT RUN {
## 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()
# }