--- title: "Measles: Building the He et al. Model" draft: true execute: echo: false --- ```{python} #| output: false import subprocess, io, os, sys, re import polars as pl import numpy as np import matplotlib.pyplot as plt import matplotlib as mpl from IPython.display import HTML sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname("__file__"), ".."))) from _lib.cli import run_cli GREY = "#888888" mpl.rcParams.update({ "figure.facecolor": "none", "axes.facecolor": "none", "savefig.facecolor": "none", "figure.dpi": 150, "font.size": 11, "text.color": GREY, "axes.labelcolor": GREY, "xtick.color": GREY, "ytick.color": GREY, "axes.edgecolor": GREY, "legend.edgecolor": GREY, }) VDIR = "he2010" def run_sim(model, params=None, seed=42, extra_args=None, cwd=VDIR): cmd = ["camdl", "simulate", model, "--seed", str(seed)] if params: if isinstance(params, dict): for k, v in params.items(): cmd.extend(["--param", f"{k}={v}"]) else: cmd.extend(["--params", params]) if extra_args: cmd.extend(extra_args) result = subprocess.run(cmd, capture_output=True, text=True, check=True, cwd=cwd) lines = result.stdout.split("\n") tsv = "\n".join(l for l in lines if not l.startswith("#")) return pl.read_csv(io.StringIO(tsv), separator="\t") ``` The [previous chapters](../guide/fitting-data/report.html) fit an SIR model to 14 days of a boarding-school influenza outbreak — a closed population, two free parameters, a single epidemic wave. That was the simplest interesting inference problem. Measles in London is a different beast. The 1978 boarding school had 763 boys and lasted two weeks. London in the 1950s had 3.5 million people and measles came back every year or two for decades. To model that, we need machinery the boarding-school SIR doesn't have: births that refill the susceptible pool, a latent period, seasonal forcing from the school calendar, and an observation model that accounts for incomplete reporting. This chapter builds the He et al. (2010) London measles SEIR step by step. Each section adds one mechanism and shows what it does to the dynamics. By the end, every line of the full model is recognizable. ## From outbreak to endemic: why demography matters The boarding-school SIR is a **closed** population — no one enters or leaves. The epidemic burns through the susceptible pool and dies out. That's fine for a two-week flu. But measles in a city persists for decades because new children are born susceptible. Adding births and deaths to the SIR is the single biggest conceptual step in this chapter: it turns an outbreak model into an **endemic** model. ```{python} #| output: asis print("```camdl") with open(os.path.join(VDIR, "models/sir_endemic.camdl")) as f: print(f.read().rstrip()) print("```") ``` Three lines changed from the boarding-school SIR: - `birth : --> S @ mu * N0` — new susceptibles enter at rate μ × N - `death_S : S --> @ mu * S` — individuals leave each compartment at the same rate - The population stays roughly constant: births ≈ deaths when both happen at rate μ. The death rate μ = 0.02/yr (≈ 0.0000548/day) is tiny on any given day, but over 50 years it completely replaces the population. For measles with R₀ ≈ 15–20, the susceptible pool drains in weeks but refills over years — and that interplay is what generates recurrent epidemics. ```{python} #| label: fig-sir-endemic #| fig-cap: "SIR with demography at two R₀ values over 50 years (N = 500,000). Susceptibles (blue) drain during each epidemic, then slowly refill via births. Higher R₀ (right) produces shorter inter-epidemic periods because the susceptible threshold for a new outbreak is reached sooner." fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 4), sharey=True) for ax, R0, title in [(ax1, 8.0, "R₀ = 8"), (ax2, 16.0, "R₀ = 16")]: N0 = 500000 gamma = 1/12.0 # 12-day infectious period mu = 0.02/365.25 beta = R0 * (gamma + mu) traj = run_sim("models/sir_endemic.camdl", params={"beta": beta, "gamma": gamma, "mu": mu, "N0": N0}, seed=7) t_yr = traj["t"].to_numpy() / 365.25 ax.plot(t_yr, traj["S"].to_numpy(), color="#3366cc", linewidth=0.5, alpha=0.8, label="S") ax.plot(t_yr, traj["I"].to_numpy(), color="#cc3333", linewidth=0.5, alpha=0.8, label="I") ax.set_xlabel("Years") ax.set_title(title, fontsize=11, color=GREY) ax.spines[["top", "right"]].set_visible(False) ax.legend(frameon=False, fontsize=9) ax1.set_ylabel("Count") plt.tight_layout() plt.show() ``` Compare this to the boarding-school SIR: one epidemic, done. Here, susceptibles accumulate via births until the epidemic threshold is crossed, the epidemic fires, and the cycle repeats. The inter-epidemic period depends on R₀ — higher R₀ means fewer susceptibles needed to trigger an outbreak, so epidemics recur faster. This is already a useful model of endemic disease. But it's missing features that matter for measles specifically. ## Adding a latent period: from SIR to SEIR Measles has an incubation period of 8–13 days — you're infected but not yet infectious. The SIR model doesn't capture this: infection and infectiousness happen simultaneously. The **SEIR** model adds an **Exposed** (E) compartment where individuals incubate before moving to I. ```{python} #| output: asis print("```camdl") with open(os.path.join(VDIR, "models/seir_endemic.camdl")) as f: for line in f: if line.strip().startswith("#") and "SEIR" in line: print(line.rstrip()) elif "compartments" in line or "latency" in line or "sigma" in line: print(line.rstrip()) print("```") ``` The new pieces: `compartments { S, E, I, R }`, a `sigma : rate` parameter (1/σ = mean incubation period), and the transition `latency : E --> I @ sigma * E`. Everything else is identical to the SIR. ```{python} #| label: fig-seir-vs-sir #| fig-cap: "SEIR (left) vs SIR (right) with demography. Same R₀, same population. The latent period delays and broadens each epidemic wave — exposed individuals incubate for ~12 days before becoming infectious, which smooths the dynamics." fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 4), sharey=True) N0 = 500000 gamma = 1/12.0 sigma = 1/12.0 # 12-day incubation mu = 0.02/365.25 R0 = 12.0 beta = R0 * (gamma + mu) # SEIR traj_seir = run_sim("models/seir_endemic.camdl", params={"beta": beta, "gamma": gamma, "sigma": sigma, "mu": mu, "N0": N0}, seed=7) t_yr = traj_seir["t"].to_numpy() / 365.25 ax1.plot(t_yr, traj_seir["I"].to_numpy(), color="#cc3333", linewidth=0.5, alpha=0.8) ax1.plot(t_yr, traj_seir["E"].to_numpy(), color="#e67e22", linewidth=0.5, alpha=0.4, label="E") ax1.set_title("SEIR (σ = 1/12 day⁻¹)", fontsize=11, color=GREY) ax1.set_xlabel("Years") ax1.set_ylabel("Count") ax1.spines[["top", "right"]].set_visible(False) # SIR traj_sir = run_sim("models/sir_endemic.camdl", params={"beta": beta, "gamma": gamma, "mu": mu, "N0": N0}, seed=7) ax2.plot(traj_sir["t"].to_numpy() / 365.25, traj_sir["I"].to_numpy(), color="#cc3333", linewidth=0.5, alpha=0.8) ax2.set_title("SIR (no latent period)", fontsize=11, color=GREY) ax2.set_xlabel("Years") ax2.spines[["top", "right"]].set_visible(False) plt.tight_layout() plt.show() ``` The latent period acts as a low-pass filter on the dynamics: each epidemic wave is broader and lower-peaked than the SIR equivalent, because the exposed-to-infectious delay smears the surge over time. For inference, this means σ and γ interact — a longer latent period paired with a shorter infectious period can produce similar epidemic curves, creating a compensation ridge in the likelihood surface (a theme we'll return to in the fitting chapter). ## School-term forcing: where biennial epidemics come from The endemic SEIR produces recurrent epidemics, but their timing is irregular — determined by when the susceptible pool happens to cross the epidemic threshold. Real measles in pre-vaccination London showed a strikingly regular **biennial** pattern: large epidemics every other year, with smaller waves in between. The mechanism is **school-term forcing** (Fine & Clarkson 1982, Schenzle 1984). Children mix intensely during school terms and barely at all during holidays. The transmission rate β(t) oscillates between a high value (term-time) and a low value (holidays) with the school calendar. ```{python} #| label: fig-school-forcing #| fig-cap: "UK school calendar as a forcing function. The transmission rate scales by 1 + a during term (grey bands) and 1 - a × (277/88) during holidays, where a is the amplitude parameter. With a = 0.554 (the He et al. MLE), term-time transmission is ~50% higher than average." days = np.arange(0, 366) # UK school terms: days 7–100, 115–199, 252–300, 308–356 school = np.zeros_like(days, dtype=float) for lo, hi in [(7,100), (115,199), (252,300), (308,356)]: school[(days >= lo) & (days <= hi)] = 1.0 amplitude = 0.554 seas = 1.0 - amplitude + amplitude * (1.0 + 0.2411/0.7589) * school fig, ax = plt.subplots(figsize=(10, 3)) ax.fill_between(days, 0, 1, where=school > 0, alpha=0.08, color="#3366cc", transform=ax.get_xaxis_transform()) ax.plot(days, seas, color="#cc3333", linewidth=2) ax.axhline(1.0, color=GREY, linewidth=0.5, linestyle=":") ax.set_xlabel("Day of year") ax.set_ylabel("β(t) / β̄") ax.set_xlim(0, 365) ax.spines[["top", "right"]].set_visible(False) # Label holidays for x, label in [(53, "Spring\nterm"), (157, "Summer\nterm"), (276, "Autumn\nterm"), (332, "Winter\nterm")]: ax.text(x, 0.3, label, ha="center", fontsize=8, color="#3366cc", alpha=0.6) for x, label in [(107, "Easter"), (225, "Summer"), (304, "Half\nterm")]: ax.text(x, 1.55, label, ha="center", fontsize=8, color=GREY, alpha=0.7) plt.tight_layout() plt.show() ``` In camdl, school-term forcing uses a `periodic` block in the `forcing {}` section: ```camdl forcing { school : periodic { period = 365.25 'days step = 1 'days on = [7:100, 115:199, 252:300, 308:356] } } ``` The `school(t)` function returns 1 during term days and 0 during holidays. The transmission rate then modulates as `β_base × seas` where `seas = 1 - a + a × (1 + 0.2411/0.7589) × school(t)`, ensuring the time-averaged forcing equals 1. ```{python} #| label: fig-seir-seasonal #| fig-cap: "SEIR with school-term forcing (amplitude = 0.554) over 30 years, N = 3,000,000. The interaction between seasonal forcing and susceptible depletion produces the characteristic biennial pattern: a large epidemic depletes susceptibles below the winter threshold, so the following year's outbreak is smaller, allowing susceptibles to accumulate for a large epidemic the year after." N0 = 3000000 traj = run_sim("models/seir_seasonal.camdl", params={"R0": 17.0, "sigma": 1/12.0, "gamma": 1/12.0, "mu": 0.02/365.25, "amplitude": 0.554, "N0": N0}, seed=3) fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 6), sharex=True, gridspec_kw={"height_ratios": [1, 2]}) t_yr = traj["t"].to_numpy() / 365.25 mask = t_yr < 30 ax1.plot(t_yr[mask], traj["S"].to_numpy()[mask] / N0 * 100, color="#3366cc", linewidth=0.8) ax1.set_ylabel("Susceptible (%)") ax1.spines[["top", "right"]].set_visible(False) ax2.plot(t_yr[mask], traj["I"].to_numpy()[mask], color="#cc3333", linewidth=0.8) ax2.set_xlabel("Years") ax2.set_ylabel("Infectious (I)") ax2.spines[["top", "right"]].set_visible(False) plt.tight_layout() plt.show() ``` The biennial pattern emerges from the **interaction** between two clocks: the annual school calendar and the multi-year susceptible replenishment cycle. A large epidemic in year 1 depletes susceptibles so far below the threshold that the next winter's forcing can't trigger a full outbreak — only a small wave. That leaves more susceptibles for year 3, producing a large epidemic again. The period depends on R₀, the amplitude of forcing, and the birth rate — and it's not always biennial. At different parameter values the model produces annual, biennial, or chaotic dynamics, which is part of what makes inference hard. ## The London measles data Now we can see what we're actually trying to fit. The data are weekly notified measles cases in London, 1944–1965 — 1096 observations over 21 years. ```{python} #| label: fig-london-data #| fig-cap: "Weekly measles notifications in London, 1944–1965. The biennial pattern is clear: large epidemics in odd years with smaller inter-epidemic waves. Peak weeks reach 4,000+ cases; troughs drop below 100. This is the dataset He et al. (2010) fitted their SEIR model to." cases = pl.read_csv(os.path.join(VDIR, "data/he2010_london_cases.tsv"), separator="\t") fig, ax = plt.subplots(figsize=(12, 4)) t_yr = cases["time"].to_numpy() / 365.25 + 1944.0 # t=0 is 1944-01-07 ax.fill_between(t_yr, 0, cases["cases"].to_numpy(), color="#cc3333", alpha=0.3, linewidth=0) ax.plot(t_yr, cases["cases"].to_numpy(), color="#cc3333", linewidth=0.5) ax.set_xlabel("Year") ax.set_ylabel("Weekly cases") ax.spines[["top", "right"]].set_visible(False) plt.tight_layout() plt.show() ``` ## Can the simple SEIR explain this data? Before adding more machinery, it's worth asking: how far does the seasonal SEIR already get? We can simulate the seasonal model at roughly measles-like parameters (R₀ ≈ 17, 12-day latent and infectious periods, amplitude ≈ 0.55) and overlay on the observed data. The model doesn't have reporting (ρ), so we compare the shape and timing rather than the absolute scale. ```{python} #| label: fig-seir-vs-data #| fig-cap: "Seasonal SEIR simulation (blue, right axis) overlaid on London measles data (red, left axis). The model produces biennial epidemics with roughly the right timing, but there are clear gaps: the simulated peaks are too regular, the amplitude doesn't vary across decades, and the absolute scale requires a reporting fraction to match. These gaps motivate the full He et al. model." N0 = 3500000 traj = run_sim("models/seir_seasonal.camdl", params={"R0": 17.0, "sigma": 1/12.0, "gamma": 1/12.0, "mu": 0.02/365.25, "amplitude": 0.554, "N0": N0}, seed=5) # Aggregate model I to weekly for comparison t_days = traj["t"].to_numpy() I_vals = traj["I"].to_numpy() fig, ax1 = plt.subplots(figsize=(12, 5)) ax2 = ax1.twinx() # Data t_data_yr = cases["time"].to_numpy() / 365.25 + 1944.0 ax1.fill_between(t_data_yr, 0, cases["cases"].to_numpy(), color="#cc3333", alpha=0.2, linewidth=0) ax1.plot(t_data_yr, cases["cases"].to_numpy(), color="#cc3333", linewidth=0.5, alpha=0.7, label="Observed cases") ax1.set_ylabel("Weekly reported cases", color="#cc3333") ax1.tick_params(axis="y", colors="#cc3333") # Model — shift to align with data years (model starts at t=0, data at 1944) t_model_yr = t_days / 365.25 + 1944.0 ax2.plot(t_model_yr, I_vals, color="#3366cc", linewidth=0.5, alpha=0.6, label="Simulated I(t)") ax2.set_ylabel("Simulated infectious (I)", color="#3366cc") ax2.tick_params(axis="y", colors="#3366cc") ax1.set_xlabel("Year") ax1.set_xlim(1944, 1965) ax1.spines["top"].set_visible(False) ax2.spines["top"].set_visible(False) # Combined legend lines1, labels1 = ax1.get_legend_handles_labels() lines2, labels2 = ax2.get_legend_handles_labels() ax1.legend(lines1 + lines2, labels1 + labels2, frameon=False, fontsize=9, loc="upper right") plt.tight_layout() plt.show() ``` The seasonal SEIR gets the basic biennial pattern right — epidemics recur roughly every two years, peaking in winter. But the comparison reveals exactly what's missing: - **Scale.** The model simulates true infections; the data are *reported* cases. He et al. estimate ρ ≈ 0.49, meaning roughly half of infections go unreported. Without a reporting fraction, we can't compare counts directly. - **Amplitude variation.** The model's epidemic peaks are roughly constant; the data shows a declining trend (the 1940s epidemics are much larger than the 1960s). This comes from London's changing population and birth rate — time-varying covariates the simple model doesn't have. - **Stochastic variation.** Even within the biennial cycle, real epidemics vary in size and timing in ways that a single stochastic trajectory can't match consistently. The `overdispersed()` primitive adds day-to-day noise in the force of infection. - **Irregular inter-epidemic troughs.** The data shows non-zero cases even between epidemics. The importation rate ι (a few cases/year from outside London) sustains low-level transmission during troughs. The next sections add these features to reach the full He et al. model. ## Fitting the seasonal SEIR: prior predictive check Before fitting, we should check what the model predicts *before seeing data*. A prior predictive check draws parameters from the prior and simulates — if the prior envelope doesn't overlap the data at all, the fit will fail before it starts. We use measles-informed ranges: R₀ ∈ [10, 60], amplitude ∈ [0.1, 0.9], ρ ∈ [0.1, 0.9], s₀ ∈ [0.01, 0.12]. These aren't arbitrary — they reflect the literature on measles dynamics in pre-vaccination populations. ```{python} #| label: fig-prior-predictive #| fig-cap: "Prior predictive check: 40 simulations at uniformly drawn parameters from the prior ranges, with weekly incidence scaled by a random ρ. Grey lines are individual trajectories; the red shaded area is the observed data. The prior envelope covers the data — the model can, in principle, generate epidemics of the right magnitude and timing." fig, ax = plt.subplots(figsize=(12, 5)) # Plot data first t_data_yr = cases["time"].to_numpy() / 365.25 + 1944.0 ax.fill_between(t_data_yr, 0, cases["cases"].to_numpy(), color="#cc3333", alpha=0.2, linewidth=0) ax.plot(t_data_yr, cases["cases"].to_numpy(), color="#cc3333", linewidth=0.8, alpha=0.6, label="Observed cases", zorder=3) rng = np.random.default_rng(42) for i in range(40): R0_draw = rng.uniform(10, 60) amp_draw = rng.uniform(0.1, 0.9) rho_draw = rng.uniform(0.1, 0.9) s0_draw = rng.uniform(0.01, 0.12) try: traj = run_sim("models/seir_seasonal.camdl", params={"R0": R0_draw, "sigma": 0.0791, "gamma": 0.0832, "mu": 0.02/365.25, "amplitude": amp_draw, "N0": 3500000}, seed=i+1) t_yr = traj["t"].to_numpy() / 365.25 + 1944.0 # Weekly incidence proxy: difference in R scaled by rho R_vals = traj["R"].to_numpy() weekly_inc = np.diff(R_vals[::7]) * rho_draw t_weekly = t_yr[::7][1:] if len(t_weekly) > 0: ax.plot(t_weekly, weekly_inc, color="#555555", linewidth=0.5, alpha=0.45, zorder=1) except Exception: pass ax.set_xlabel("Year") ax.set_ylabel("Weekly cases") ax.set_xlim(1944, 1965) ax.set_ylim(0, 15000) ax.legend(frameon=False, fontsize=9) ax.spines[["top", "right"]].set_visible(False) plt.tight_layout() plt.show() ``` ## Fitting with IF2 We add a NegBin observation model to the seasonal SEIR (weekly reported recoveries, reporting fraction ρ, overdispersion k) and fit to the London data with IF2. Five parameters are estimated — R₀, amplitude, ρ, k, and the initial susceptible fraction s₀ — with σ, γ, μ, and N₀ fixed at literature values. ```{python} #| echo: false #| output: asis with open(os.path.join(VDIR, "fit_seir_seasonal.toml")) as f: print(f"```toml\n{f.read().strip()}\n```") ``` ```{python} #| echo: false display(run_cli( f"cd {VDIR} && camdl fit run fit_seir_seasonal.toml --seed 42", echo="camdl fit run fit_seir_seasonal.toml --seed 42", cwd=VDIR, max_lines=25 )) ``` ### Scout convergence ```{python} #| label: fig-scout-traces #| fig-cap: "IF2 scout traces for the seasonal SEIR on London measles. 16 chains × 1000 particles × 150 iterations. R₀ and k are multimodal ( > 10) — chains occupy different basins. Amplitude converges; ρ is marginal ( ≈ 1.3). The simple seasonal SEIR without importation or process noise has severe identifiability problems on 21 years of data." import glob as _glob scout_dir = sorted(_glob.glob(os.path.join(VDIR, "results/fits/fit_seir_seasonal-*/real/fit_42/scout")))[-1] traces = [] for cdir in sorted(_glob.glob(os.path.join(scout_dir, "chain_*"))): if not os.path.isdir(cdir): continue # Skip comment lines before reading trace_path = os.path.join(cdir, "parameter_traces.tsv") with open(trace_path) as _f: tsv_lines = [l for l in _f if not l.startswith("#")] df = pl.read_csv(io.StringIO("".join(tsv_lines)), separator="\t", infer_schema_length=10000) cid = int(os.path.basename(cdir).split("_")[1]) df = df.with_columns(pl.lit(cid).alias("chain")) traces.append(df) traces_df = pl.concat(traces) params = ["R0", "amplitude", "rho", "k", "s0"] fig, axes = plt.subplots(1, 5, figsize=(16, 3.5)) for ax, param in zip(axes, params): for chain in traces_df["chain"].unique().sort().to_list(): sub = traces_df.filter(pl.col("chain") == chain) ax.plot(sub["iteration"], sub[param], alpha=0.35, linewidth=0.4, color="#3366cc") ax.set_ylabel(param) ax.set_xlabel("Iteration") ax.spines[["top", "right"]].set_visible(False) plt.suptitle("Scout traces — seasonal SEIR on London measles", fontsize=11, color=GREY, y=1.02) plt.tight_layout() plt.show() ``` ```{python} #| echo: false #| output: asis # Read  from fit_state.toml import tomllib state_path = os.path.join(scout_dir, "fit_state.toml") if os.path.exists(state_path): with open(state_path, "rb") as f: state = tomllib.load(f) agreements = state.get("tail_chain_agreement", {}) best_ll = state.get("best_loglik", None) if best_ll is not None: print(f"**Best scout log-likelihood:** {best_ll:.1f} (chain {state.get('best_chain', '?')})") print() print("| parameter |  | status |") print("|---|---|---|") for param in ["R0", "amplitude", "rho", "k", "s0"]: agreement = agreements.get(param, 0) if agreement == 0: continue color = "#27ae60" if agreement < 1.05 else ("#e67e22" if agreement < 1.10 else "#c0392b") status = "✓" if agreement < 1.05 else ("~" if agreement < 1.10 else "✗") print(f"| `{param}` | " f"{agreement:.3f} | {status} |") ``` ### What the diagnostics tell us The scout found a best log-likelihood of about −11,300 — the model can explain some of the data, but not well. The convergence diagnostics are worse: R₀ and k have  values above 10, meaning chains are exploring completely different basins. This is the seasonal SEIR telling us it's not the right model for this data. The specific failures point to specific gaps: - **R₀ multimodality** — without importation (ι), inter-epidemic troughs go to exactly zero infections, and the model can match the data equally well at very different R₀ values by adjusting the timing of extinction and re-seeding. - **k wandering** — the NegBin overdispersion parameter absorbs variance that should come from process noise (`overdispersed()`) and time-varying covariates. Without those, k has to do too much work. - **s₀ at lower bound** — the initial susceptible fraction is pinned at 1%, suggesting the model needs a different initialization strategy (warm-up period or covariate-informed initial conditions). These are not fitting failures — they're the model telling us what it needs. Each gap maps to a specific feature in the full He et al. model. ## What's still needed: the full He et al. model The seasonal SEIR gets the biennial pattern right but can't fit the data quantitatively. The remaining mechanisms — time-varying covariates, importation, process noise, cohort births, and the He et al. observation model — address the specific failures we just diagnosed: **Time-varying covariates.** `interpolated` forcing for pop(t) and birthrate(t) replaces the fixed N₀ and μ: ```camdl forcing { pop : interpolated { data = "data/he2010_london_covariates.tsv" time_col = t; value_col = pop; method = "linear" } birthrate : interpolated { data = "data/he2010_london_covariates.tsv" time_col = t; value_col = birthrate; method = "linear" } } ``` **Importation** (ι ≈ 3 cases/year from outside London) sustains low-level transmission during inter-epidemic troughs — without it, the model goes to zero and can only restart via stochastic fluctuation. **Process noise** via `overdispersed()` adds Gamma-distributed noise to the force of infection. This is the He et al. "environmental stochasticity" that we saw outperform NegBin observation noise in the [model comparison chapter](../guide/fitting-data/model_comparison.html). **Cohort births** via `events {}` — children enter school once per year, creating an annual pulse of new susceptibles. **The observation model** — weekly reported recoveries with reporting fraction ρ and overdispersion ψ (He et al. eq. 6). ```{python} #| output: asis print("The complete `he2010_london.camdl`:") print() print("```camdl") with open(os.path.join(VDIR, "models/he2010_london.camdl")) as f: text = f.read() lines = text.split("\n") start = next(i for i, l in enumerate(lines) if l.strip() and not l.strip().startswith("#")) print("\n".join(lines[start:]).rstrip()) print("```") ``` Twelve parameters, time-varying forcing from three sources, and a non-trivial observation model. The fitting chapter will show what happens when IF2 meets this full model — which parameters converge, which form ridges, and what the diagnostics look like on a real-world inference problem.