## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE, comment = "#>",
  fig.width = 7, fig.height = 5
)
library(npfseir)
set.seed(42)

## ----simulate-----------------------------------------------------------------
fixed <- list(
  eps        = 1/5,        # incubation rate (5-day mean)
  gamma      = 1/10,       # recovery rate (10-day mean)
  delta      = 1/(70*365), # background mortality
  b          = 1/(70*365), # birth rate
  q          = 0.30,       # case detection probability
  N          = 1e6,        # reference population size
  sigma_comp = 0.01        # compartment diffusion coefficient
)

x0   <- c(S = 1e6 - 100, E = 50, I = 50, R = 0, Psi = log(0.25))
traj <- seir_simulate(n_steps = 300, kappa = 0.05, sigma_psi = 0.10,
                      mu_psi = log(0.25), x0 = x0, fixed = fixed)

plot(traj$day, traj$obs, type = "l", col = "gray40",
     xlab = "Day", ylab = "Daily confirmed cases",
     main = "Simulated epidemic (300 days)")
lines(traj$day, exp(traj$Psi) * 3000, col = "steelblue", lty = 2)
legend("topright", c("Observed counts", "True beta_t (scaled)"),
       col = c("gray40","steelblue"), lty = c(1,2), bty = "n")

## ----fit, cache=TRUE----------------------------------------------------------
obs <- traj$obs[-1]   # observations P_1, ..., P_T (exclude day-0 seed)

fit <- npf_seir(
  observations = obs,
  fixed        = fixed,
  K            = 50,    # outer (parameter) particles
  M            = 100,   # inner (state) particles
  seed         = 1
)

## ----print--------------------------------------------------------------------
print(fit)

## ----summary------------------------------------------------------------------
summary(fit, burn = 30)

## ----plot-Rt------------------------------------------------------------------
plot(fit, burn = 30)
abline(h = 1, lty = 2, col = "red")

## ----forecast-----------------------------------------------------------------
fc <- predict(fit, horizon = 14, n_draws = 500)
cat("14-day ahead forecast (mean ± 95% PI):\n")
print(round(data.frame(day = 1:14, mean = fc$mean,
                       lo = fc$lo, hi = fc$hi), 1))

## ----cori---------------------------------------------------------------------
ct <- cori_rt(obs, tau = 7, si_mean = 5.5, si_sd = 2.5)
plot(ct$t, ct$Rt_mean, type = "l", col = "darkred",
     ylim = c(0, 4), xlab = "Day", ylab = expression(R[t]),
     main = "Cori-style Rt (illustrative only)")
abline(h = 1, lty = 2)