---
title: "Getting started with npfseir"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Getting started with npfseir}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

## Overview

`npfseir` implements the online Bayesian inference framework of Feng and Wang (2025)
for stochastic SEIR epidemic models with a time-varying transmission rate.
The transmission rate $\beta_t$ is modelled as a latent Ornstein--Uhlenbeck
(OU) process in log-space, and inference is performed via the nested particle
filter (NPF) of Crisan and Miguez (2018).

## 1. Simulate a synthetic epidemic

```{r 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")
```

## 2. Run the nested particle filter

The `npf_seir()` function takes only the observation vector and the fixed
parameter list. By default it uses $K = 100$ outer particles and $M = 200$
inner particles; we use smaller values here for speed.

```{r 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
)
```

## 3. Inspect results

```{r print}
print(fit)
```

```{r summary}
summary(fit, burn = 30)
```

## 4. Plot the R_t trajectory

```{r plot-Rt}
plot(fit, burn = 30)
abline(h = 1, lty = 2, col = "red")
```

## 5. Forecast

```{r 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))
```

## 6. Cori-style R_t comparison (illustrative only)

`cori_rt()` implements an in-house Cori-style renewal estimator for
**illustrative comparison of estimands** only. It is *not* the R `EpiEstim`
package.

```{r 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)
```

## References

Crisan, D. and Miguez, J. (2018). Nested particle filters for online
parameter estimation in discrete-time state-space Markov models.
*Bernoulli*, 24(4A):3039--3086.

Feng, W. and Wang, W. (2025). Bayesian sequential inference for a stochastic SEIR model
with latent transmission dynamics. *Preprint*.