renewal model time series example for User's Guide

Bắt đầu bằng cách xem xét cấu trúc của User's Guide cùng với mô hình Stan và trình mô phỏng CmdStanPy được đưa vào issue. Thêm ví dụ về mô hình tái tạo cùng mã mô phỏng và các tài liệu tham khảo Fraser và Steyn được trích dẫn, sau đó xác minh rằng ví dụ được đặt trong phần phù hợp của hướng dẫn và các liên kết cũng như mã có thể sử dụng được.

We should include a renewal model---they're popular for tracking infectious diseases. Here's the basic model with Stan code and CmdStanPy simulation.

functions {
  real random_walk_lpdf(vector y, real sigma, real mu0, real sigma0) {
    int N = rows(y);
    return normal_lpdf(y[1] | mu0, sigma0)
           + normal_lpdf(y[2:N] | y[1:N - 1], sigma);
  }
}
data {
  int<lower=0> T, U;
  simplex[U] omega;
  array[T]  int<lower=0> C;
  real mu0;
  real<lower=0> sigma0;
}
transformed data {
  row_vector<lower=0>[T] Cv = to_row_vector(C);
}
parameters {
  vector<lower=0>[T] R;
  real<lower=0> sigma;
}
model {
  // C[1] not modeled
  for (t in 2:U) {
    // shorten history with renormalized simplex
    C[t] ~ poisson(R[t] * (Cv[1:t-1] * omega[1:t-1]) / sum(omega[1:t-1]));
  }
  for (t in U + 1:T) {
    C[t] ~ poisson(R[t] * (Cv[t-U:t-1] * omega));
  }
  sigma ~ lognormal(0, 0.25);
  R ~ random_walk(sigma, mu0, sigma0);
}

And here's the simulator.

import numpy as np
import cmdstanpy as csp


def simulate_renewal_data(
    T: int = 150,
    U: int = 5,
    seed: int = 2020,
):
    rng = np.random.default_rng(seed)

    alpha = np.arange(U, 0, -1, dtype=float)
    omega = alpha / sum(alpha)

    mu0 = 1.0
    sigma0 = 0.05
    sigma = 0.025

    R = np.empty(T)
    R[0] = abs(rng.normal(mu0, sigma0))
    for t in range(1, T):
        R[t] = abs(rng.normal(R[t - 1], sigma))

    C = np.zeros(T, dtype=int)
    C[0] = 1000
    for t in range(1, U):
        lam = R[t] * np.dot(C[0:t], omega[0:t] / sum(omega[0:t]))
        C[t] = rng.poisson(lam)
    for t in range(U, T):
        lam = R[t] * np.dot(C[t-U:t], omega)
        C[t] = rng.poisson(lam)

    return {
        "T": T,
        "U": U,
        "omega": omega,
        "C": C,
        "R": R,
        "sigma": sigma,
        "mu0": mu0,
        "sigma0": sigma0,
        "alpha": alpha,
    }

sim = simulate_renewal_data()

print("alpha: ")
print(np.round(sim["alpha"], 3))
print("\nomega: ")
print(np.round(sim["omega"], 3))
print("\nsigma: ", sim["sigma"])
print("\nR: ")
print(np.round(sim["R"], 3))
print("\nC: ")
print(sim["C"])

model = csp.CmdStanModel(stan_file="poisson.stan")
fit = model.sample(data=sim)
print(fit.summary())

The model is taken from Steyn et al. (2025), Model 1 (section 6.1), which also introduces overdispersed alternatives (negative binomial instead of Poisson), latent infectiousness states, day-of-week effects, etc. They cite Fraser (2007) as the source of the basic model.

@article{steyn2025primer,
  title={A primer on inference and prediction with epidemic renewal models and sequential {M}onte {C}arlo},
  author={Steyn, Nicholas and Parag, Kris V and Thompson, Robin N and Donnelly, Christl A},
  journal={Statistics in Medicine},
  volume={44},
  number={18-19},
  pages={e70204},
  year={2025}
}

@article{fraser2007estimating,
  title={Estimating individual and household reproduction numbers in an emerging epidemic},
  author={Fraser, Christophe},
  journal={PloS one},
  volume={2},
  number={8},
  pages={e758},
  year={2007}
}

Here are links to the papers:

• https://onlinelibrary.wiley.com/doi/pdf/10.1002/sim.70204
• https://journals.plos.org/plosone/article/file?id=10.1371/journal.pone.0000758&type=printable