accumvars: accumulator variables

In pomp: Statistical Inference for Partially Observed Markov Processes

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, 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, 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()

pomp documentation built on Sept. 13, 2024, 1:08 a.m.