Inference in camdl

How the particle filter, IF2, PGAS, and NUTS work — including maximum-likelihood and Bayesian fitting on the ODE backend — what the diagnostics mean, and how the inference pipeline fits together.

---

The inference problem

A compartmental model defines a stochastic process over a latent state — the compartment populations (S, E, I, R) evolving through time via stochastic transitions. You never observe this state directly. What you observe is a noisy, incomplete projection: weekly case reports, which are some fraction of recoveries plus measurement noise.

The goal of inference is to estimate model parameters (transmission rate, recovery rate, reporting probability, etc.) from the observed data. This requires evaluating the likelihood — the probability of the observed data given the parameters:

p(y_{1:T} \mid \theta) = \int p(y_{1:T} \mid x_{0:T}, \theta) \, p(x_{0:T} \mid \theta) \, dx_{0:T}

This integral is over all possible latent state trajectories x_{0:T}. For a compartmental model with thousands of individuals tracked over hundreds of weekly observations, this integral is intractable — you can’t evaluate it analytically.

But you can simulate trajectories from p(x_{0:T} \mid \theta). That’s what camdl does: given parameters, generate a stochastic realization of (S, E, I, R) over time. The particle filter exploits this simulation ability to estimate the intractable integral via Monte Carlo.

Why this is hard

Three things make compartmental model inference harder than standard statistical problems:

1. Intractable likelihood. The stochastic process (chain-binomial, Gillespie) doesn’t have an analytic transition density. You can simulate from it but you can’t evaluate p(x_t \mid x_{t-1}, \theta) in closed form. This rules out MCMC methods that require pointwise likelihood evaluation.

2. High-dimensional latent state. The state at each timepoint is the full vector of compartment populations. Over T observations, the latent trajectory is T-dimensional. The likelihood integral is over this entire path space.

3. Nonlinear dynamics. Small parameter changes can produce qualitatively different behavior — biennial vs annual epidemic cycles, fadeout vs persistence, early vs late peak timing. The likelihood surface has ridges, local optima, and flat regions where different parameter combinations produce similar dynamics.

How the algorithms relate

The three particle methods in camdl — IF2, PMMH and PGAS — are built on Sequential Monte Carlo (SMC), the particle filter. They differ in what they do with the particles and what they produce. (The ode-backend methods, mh and the nl-* optimizers, need no filter: they fit the deterministic ODE’s likelihood, which is computed directly.)

IF2 (Iterated Filtering) perturbs parameters inside the particle filter and cools toward the MLE. It’s a stochastic optimization algorithm, not a sampler — it finds the best-fit parameters but doesn’t characterize uncertainty. Fast, robust, good for finding the right basin.

PMMH (Particle Marginal Metropolis-Hastings) uses the particle filter’s log-likelihood estimate as the acceptance ratio in a Metropolis sampler. It marginalizes out trajectories — the PF integrates over all possible latent state paths, and PMMH only sees the marginal likelihood number \hat{p}(y|\theta). Any process model works (plug-and-play), but the PF likelihood estimate is noisy, which slows mixing.

PGAS (Particle Gibbs with Ancestor Sampling) conditions on a specific trajectory. It holds one complete latent trajectory X fixed and evaluates the exact complete-data likelihood p(y, X | \theta) — no estimation noise. Parameters are updated via NUTS or MH using this exact likelihood. The trajectory is then refreshed via CSMC-AS (a particle filter conditioned on the old trajectory). The Gibbs alternation (\theta | X then X | \theta) samples from the full joint posterior p(\theta, X | y).

Marginalizing vs conditioning on trajectories

This is the fundamental design choice:

	PMMH	PGAS
Trajectories	Marginalized out by PF	Conditioned on, explicitly sampled
Likelihood	Estimated (noisy)	Exact (no PF variance)
Process model	Any (plug-and-play)	Chain-binomial only (needs transition density)
Output	Posterior p(\theta \mid y) (trajectories available but low-quality due to path degeneracy)	Posterior p(\theta, X \mid y) jointly (high-quality trajectory samples via CSMC)
Bottleneck	PF variance → slow mixing	Trajectory convergence → slow on long series

PMMH is more general (works with any simulator) but pays for it with noisy likelihood estimates. PGAS is more efficient (exact likelihood) but requires the ability to evaluate transition densities — currently only the chain-binomial (Euler-multinomial) backend supports this.

Which one to use

PGAS is the default for the posterior of a chain-binomial model, and PMMH is a second sampler to check it against; IF2 gives a point estimate and, through camdl profile, profile likelihoods. The table, and the main reason to prefer a posterior (it shows parameters the data cannot separate, where a point estimate returns one point and hides them), are in camdl docs workflow, “Choose a method”.

Chains start from draws from the priors by default. If those starts are slow to reach the region where the likelihood is high, a short run can seed a longer one:

short.toml    [method] algorithm = "pgas", sweeps = 200                → @short
fit.toml      [method] algorithm = "pgas", starts = { from_posterior = "@short" }

from_posterior draws one start per chain from the short run’s posterior, so the chains still start apart and R̂ stays meaningful. from_mle = "@mle" starts every chain at one IF2 point instead, and R̂ is then reported as not assessed.

The particle filter and these three methods all target the stochastic-process likelihood p(y \mid \theta) on the chain-binomial backend. The ODE backend fits a different statistical object — the deterministic marginal likelihood — with its own methods: the NLopt optimizers for an MLE and gradient-free mh for a posterior; see Fitting on the ODE backend, below.

The particle filter

The particle filter (sequential Monte Carlo) estimates the likelihood by running many parallel simulations and letting the data select which ones survive.

The key insight: instead of integrating over all possible state trajectories at once, do it sequentially — one observation at a time. At each observation, use importance sampling to focus computational effort on trajectories that are consistent with the data seen so far.

Particles are state trajectories

Each of the N particles is an independent stochastic simulation of the full compartmental model. At any time t, particle i has its own state vector (S_i, E_i, I_i, R_i)_t — its own realization of the epidemic. The particles all share the same parameters \theta but differ in their random draws (which individuals get infected, when they recover, etc.).

The ensemble of N particles approximates the filtering distribution p(x_t \mid y_{1:t}, \theta) — the posterior over the latent state given all data up to time t.

Weights score particles against data

At each observation time t, each particle i gets a weight proportional to how well it predicts the observed data:

w_i^{(t)} = p(y_t \mid x_i^{(t)}, \theta)

This is the observation model likelihood — for example, the discretized Normal probability of seeing 500 reported cases given that particle i’s projected recoveries (scaled by reporting rate \rho) predicted 490.

If particle i predicted well, w_i is large. If it predicted poorly (e.g., projected 50 cases when 500 were observed), w_i is tiny.

Resampling focuses effort

After weighting, bootstrap resampling draws N new particles from the current N, with probability proportional to weights. Particles that predicted well get duplicated. Particles that predicted poorly are discarded.

After resampling, all particles have equal weight, but they cluster around state trajectories that are consistent with the data. The filter has used the observation to update its belief about the latent state — this is Bayesian updating via Monte Carlo.

The likelihood estimate

The marginal likelihood of each observation is the mean weight:

\hat{p}(y_t \mid y_{1:t-1}, \theta) = \frac{1}{N} \sum_{i=1}^{N} w_i^{(t)}

The total log-likelihood is the sum over all observations:

\hat{\ell}(\theta) = \sum_{t=1}^{T} \log \hat{p}(y_t \mid y_{1:t-1}, \theta)

This estimate is unbiased (in expectation, it equals the true log-likelihood). With more particles, the variance decreases. The estimate is always a lower bound on the true log-likelihood — more particles can only improve it.

Effective sample size (ESS)

After weighting but before resampling, the weights are unequal. The effective sample size measures how many particles are actually contributing useful information:

\text{ESS}_t = \frac{\left(\sum_i w_i^{(t)}\right)^2}{\sum_i \left(w_i^{(t)}\right)^2}

• \text{ESS} \approx N: all weights are similar — every particle is useful. The observation is unsurprising given the model.
• \text{ESS} \approx 1: one particle has almost all the weight — the filter has degenerated. Only one trajectory out of N is consistent with the data. The log-likelihood estimate is unreliable.

ESS is the primary diagnostic. It drops during epidemic peaks (where the data is most informative and small differences in predicted incidence produce large weight differences) and recovers during inter-epidemic troughs (where all particles predict similar low incidence).

One-step-ahead predictions

Before resampling, the weighted particle ensemble gives the one-step-ahead prediction: what the filter expected to see at time t before observing y_t. The weighted mean and quantiles of \rho \times \text{projected}_i across particles give prediction intervals.

If 90% of data falls within the 90% prediction interval, the model is well-calibrated — its uncertainty is neither too wide nor too narrow. Systematic prediction bias (always overshooting peaks, always undershooting troughs) indicates model misspecification.

The mechanics above are the internal filter view. The user-facing way to get the one-step-ahead band that carries full posterior parameter uncertainty is camdl fit predict --fit fit.toml --horizon one_step: it re-runs the bootstrap filter once per posterior draw, samples ỹ from the propagated (pre-reweight) particles at each observation time, and pools over (particles × draws) into the q05…q95 columns of predictive/<stream>.tsv — horizon=one_step, treatment=posterior. This is chain-binomial only (an ODE fit’s one-step reduces to its free-forward band). See camdl docs workflow §7.

What happens at each observation time

1. PROPAGATE: advance all N particles from t_{k-1} to t_k
   For each particle i, for each sub-step dt:
     - Evaluate propensities from particle i's state
     - Draw events (multinomial for chain-binomial)
     - Accumulate flows (infection counts, recovery counts, etc.)
   After 7 sub-steps (one week), each particle has its own state
   and its own incidence count since the last observation.

2. WEIGHT: score each particle against the data
   For each particle i:
     projected_i = cumulative recovery flow since last observation
     weight_i = P(observed_cases | rho × projected_i, observation_model)
   Particles that predicted close to the observed value get high weight.
   Particles that predicted far from it get near-zero weight.

3. AGGREGATE: compute the log-likelihood increment
   ll_k = log(mean(weights))
   This is the marginal probability of this observation given all the
   particles. Sum these over all observations to get the total loglik.

4. DIAGNOSE: ESS and prediction quantiles
   ESS = 1 / sum(normalized_weights²)
   When all particles agree, ESS ≈ N. When one particle dominates,
   ESS ≈ 1. Low ESS means the filter is degenerating — most particles
   are useless and the loglik estimate is unreliable.

   Prediction quantiles (q05, q50, q95) show what the filter
   expected BEFORE seeing the data. If the data consistently falls
   outside the 90% interval, the model is misspecified.

5. RESAMPLE: keep the good particles, kill the bad ones
   Systematic resampling: select particles proportional to their
   weights. A particle with 3× the average weight gets ~3 copies.
   A particle with near-zero weight gets killed.

   After resampling, all particles are equally weighted again.
   The diversity has decreased (some particles are copies) but the
   surviving particles are all consistent with the data so far.

6. RESET: clear flow accumulators for the next observation interval

Incidence observations and the model origin

An incidence observation (incidence(...), i.e. a cumulative_flow projection) is scored against the flow accumulated over the period its row is declared to cover — the stream’s covers = ... line, or its window_start/window_stop columns (camdl docs data, “What a row covers”). Two consequences follow for the model origin (simulate.from, internal t_start).

A period cannot open before the origin. Under closing_at(time, 7 'days) a row labelled at the origin covers the seven days before it, which were never simulated; camdl pfilter / camdl fit reject the file before the filter runs, naming the period and the origin. The remedies are the obvious two: drop the row, or move simulate.from back to the period’s start or earlier. A row at the origin under day(time) or starting_on(...) covers the period starting there and is fine.

A warm-up before the first row is simulated, not scored. When the origin sits far before the data — a covariate-informed burn-in that starts dynamics years before the case series so births and campaigns shape the susceptible pool — the first period still opens where its declaration says, so the first datum is scored against one declared period and the pre-data span contributes no likelihood term. Nothing has to be set for this: the declaration is the boundary. (The former fit.toml key condition_from, which hand-placed that boundary, is removed; a leftover key is a hard error saying so — see camdl docs language-changes.)

Prevalence observations (current_pop) read state at the instant and have no period; neither point applies to them. See camdl docs fit-toml (“Where scoring begins”) and §3.9 of camdl-inference-spec.md.

Multi-cadence: streams on different schedules

Surveillance streams often arrive on different schedules — polio AFP (acute flaccid paralysis) reported roughly monthly, environmental surveillance (ES, poliovirus in sewage) sampled roughly biweekly. camdl fit / pfilter / profile accept this directly: you supply one data file per stream (each on its own time axis), and the filter merges them onto a union observation axis — the sorted set of all streams’ observation times.

Two things make this correct rather than a fudge:

• Each stream is scored only at its own observation times. At a union time where a stream has no observation (a sibling’s reporting date), that stream contributes no likelihood term — it is simply not observed then.
• Each incidence stream’s bin closes on its own cadence. The flow accumulator for an incidence stream resets only at its observations, never at a sibling’s. So an ES sampling date does not truncate AFP’s monthly count; AFP’s bin runs the full month regardless of how often ES reports in between. (Prevalence streams read state at the instant and never accumulate, so the reset doesn’t apply to them.)

Each stream’s first period opens where its own declaration says, so where scoring begins is settled per stream with no extra setting. A homogeneous model (every stream on one cadence) is the special case where the union is the shared axis and every stream is observed at every time — and it fits byte-for-byte as before.

camdl simulate --obs-dir <dir> writes one file per stream at its own cadence, which is the natural input back into a multi-cadence fit.

The --trace output

time  ll_increment  ESS    pred_mean  pred_q05  pred_q50  pred_q95  observed
7     -7.84         17.4   42.3       5         31        112       82
14    -5.37         217.7  51.2       12        45        98        98

ll_increment: How surprising this observation was. More negative = more surprising. A value of -3 means “this observation is about as likely as seeing a specific card drawn from a deck of 20.” A value of -10 means “this observation is extremely unlikely given the model.”

ESS: Effective sample size. Healthy range: 20-80% of N. Below 10% means the filter is collapsing — increase N or check the model. Above 90% means the observation is uninformative (the model already knew what to expect).

pred_q05/q50/q95: The filter’s prediction before seeing the data. If observed falls between q05 and q95 about 90% of the time, the model is well-calibrated.

The --pf-health output — “do I have enough particles?”

High-dimensional models (many spatial patches, fine age structure) can need far more particles than low-dimensional ones, because a particle filter degenerates as the effective observation dimension grows. --pf-health measures whether your particle count is adequate and estimates how many you would need.

camdl pfilter model.camdl --data cases.tsv --particles 5000 --pf-health health.tsv

It writes a per-observation TSV (time, ESS, ESS_frac, tau2) and prints a summary:

pf-health summary (156 obs, N=5000):
  ESS fraction:           median 41.2%, worst 3.1%
  log-weight variance τ²: median 1.80, max 9.4
  implied N to avoid collapse exp(τ²/2): ~2 (median step), ~110 (worst step)

ESS_frac is the realized effective sample size as a fraction of N — how many particles are actually contributing. Persistently low fractions (a worst-step a few percent of N) mean the cloud is collapsing on the tightest observations.

tau2 (τ²) is the variance of the per-particle log-weights at that observation — the predictor of how many particles you need. The ensemble size required to avoid weight collapse scales as exp(τ²/2) (Snyder, Bengtsson, Bickel & Anderson 2008, Obstacles to High-Dimensional Particle Filtering, Monthly Weather Review 136:4629–4640). The driver is the number of independently observed dimensions, not the raw state dimension: an aggregate observable keeps τ² small (a handful of particles suffice), while observing many strata independently makes τ² — and the required N — grow steeply.

The implied N line is that estimate. Treat it as an order-of-magnitude floor measured at the current parameters and N — not a guarantee. (If the filter is already collapsed you under-see the tail, so the true requirement can be larger.) When the worst-step implied N exceeds what you can afford, more particles will not rescue a fine-resolution fit on their own; the problem is the observation dimension, which calls for a different method (e.g. a block/local particle filter) rather than brute-force N.

CLI

camdl pfilter model.camdl --params p.toml --data cases.tsv \
    --particles 5000 --dt 1 --seed 42 \
    --trace -

--trace -: Write per-observation diagnostics (ll increment, ESS, the filter’s predictive quantiles vs the observed value) as a TSV. Use - for stdout, or a path to write a file.

Observation model: not a CLI flag — it is declared in the model file’s observations {} block (e.g. likelihood = neg_binomial(...)). To change the observation model, edit the model, not the command line.

Likelihood floor: when a particle predicts ~0 cases but the data shows 80, both “predicted 0” and “predicted 5” are equally wrong — the filter floors the per-particle likelihood so they are treated the same. Without the floor, “predicted 0” gets a 650 log-unit worse penalty than “predicted 5”, collapsing ESS. The default matches pomp.

Filtering marginals vs smoothing paths

At each observation step t, the bootstrap filter holds N particles weighted by p(y_t | x_t, θ). Two different distributions come out of this setup, and conflating them produces quietly-wrong plots:

• Filtering marginals p(x_t | y_{1..t}, θ) — the per-step distribution of particle states at time t, weighted by their log-weights. camdl pfilter --save-filtering PATH dumps these as a long-format TSV. Joining particles across t by index is NOT a sample path. Resampling between steps shuffles the swarm; the particle indexed i at step t+1 is not a descendant of particle i at step t.

• Smoothing paths — samples from p(x_{1:T} | y_{1:T}, θ). Each path is a coherent latent trajectory consistent with all observations. Obtained via ancestor tracing: at the final step, sample a particle proportional to its weight, walk its ancestor chain backwards to collect the state at each earlier step. camdl pfilter --save-paths N PATH writes N such paths.

For “does this fit match the data?” plots, use --save-paths. Its quantile ribbon over N paths estimates the smoothing marginal at each t — what the model believes the latent trajectory was given all the data.

For PF diagnostics (particle degeneracy, ESS decay, obs-model sanity checks, filter-implementation debugging), use --save-filtering. The per-step log-weights are what you need to detect those pathologies; they’re not what you need to compare trajectories to data.

The diagnostic plot: unconditional vs smoothing

A fitted stochastic compartmental model gives you three distinct views of the data:

1. Unconditional posterior predictive. camdl fit predict --fit fit.toml --horizon free_forward. “What does the fitted model predict a priori?” — sampled y_rep per posterior draw, pooled into the predictive/<stream>.tsv ribbon (averaged over the cloud, not run at a single plug-in point).
2. Smoothing over latent. camdl pfilter --save-paths N at the MLE. “What does the model think the latent trajectory was, given the data?”
3. Raw observations. camdl fit predict writes these to observed/<stream>.tsv in the same tidy keys, ready to overlay on (1).

Plot (1) and (2) as ribbons, (3) as points, side by side:

• If both ribbons track the data: well-specified model, inference worked.
• If (2) tracks the data but (1) misses it: diagnostic of over- flexible process noise papering over structural mis-specification. The PF log-likelihood is high because the model is flexible enough to thread through any data via stochastic fluctuations — not because it predicts well.
• If both miss the data: the fit is wrong.

The second case is pedagogically important and easy to misread. A reader seeing

1. alone miss the data will conclude “the fit is bad”; a reader seeing (2) alone track the data will conclude “the fit is good.” Neither is right. The divergence between them is the diagnostic — teach it that way.

Background: docs/dev/proposals/archive/pre-alpha/2026-04-19-pf-latent-trajectories.md.

---

IF2: turning the particle filter into an optimizer

IF2 (Iterated Filtering, Ionides et al. 2015) finds the maximum likelihood estimate (MLE) — the parameter values that make the data most probable. It does this without gradients, using only the ability to simulate forward.

The key idea

In a regular particle filter, all particles share the same parameters. In IF2, each particle carries its own parameter vector. Particle 1 might have R₀=57.2, particle 2 might have R₀=55.8. Each simulates with its own R₀.

When the filter resamples, particles with good R₀ values survive and particles with bad R₀ values die. The parameter cloud contracts around values that explain the data. Add a cooling schedule that shrinks the perturbation over time, and the cloud converges to a point — the MLE.

What happens at each observation time (IF2 vs PF)

The structure is identical to the particle filter, with two additions:

1. PERTURB: jitter each particle's parameters (NEW in IF2)
   For each particle i, for each estimated parameter:
     θ_i += Normal(0, sd_i × cooling) on the transformed scale, where
     sd_i is the natural-scale rw_sd converted to the transformed scale
     at particle i's current value (the delta method; see below)
   Parameters declared `perturb_only_at_t0` (initial conditions) are
   perturbed at t=0 only.

2. PROPAGATE: same as PF, but each particle uses its OWN params
   particle_i simulates with particle_params[i], not shared θ

3. WEIGHT: same as PF — score against data, at the SAME θ that
   just drove the propagation
4. RESAMPLE: states AND parameters are copied together (NEW in IF2)
   Good (state, θ) pairs survive. Bad pairs die.
5. RESET: same as PF

Perturbing before the propagation is not cosmetic. It is what makes one perturbed θ drive both the simulation of the latent state and the measurement density that scores it — the coupling in Algorithm 1 of Ionides et al. (2015) and in pomp’s mif2. Scoring at a θ the state was not simulated from breaks that coupling for any parameter appearing in both the process and the observation model.

The cooling schedule

The perturbation shrinks over time. After cooling_target_iters (50) full iterations of the filter, the perturbation SD is cooling_fraction (0.95) of the initial value.

Per-step cooling factor:
  c = 0.95 ^ (1 / (50 × n_observations))

After m iterations × n_obs steps each:
  effective_sd = rw_sd × c^(m × n_obs)

With 780 weekly observations, the cooling is very gentle per step (c ≈ 0.99993) but compounds over many iterations. After 50 iterations: SD is 95% of initial. After 100 iterations: ~90%. After 200: ~81%.

Early iterations: wide exploration (parameter cloud spans a broad range). Late iterations: fine tuning (cloud contracts to a tight point).

Parameter transforms and bounds

Parameters live on different scales. R₀=56.8 and ρ=0.488 need different perturbation strategies.

Parameter type	Transform	Why
positive in [a, b]	Scaled logit	Bounds enforced by construction
rate (unbounded)	Log	Multiplicative perturbation
probability in [0, 1]	Logit	Stays in (0,1)

The transform is derived automatically from the DSL parameter declaration. A parameter declared R0 : positive in [1, 100] uses scaled logit — the perturbation happens on (-∞, ∞) and the inverse transform maps back to [1, 100]. R₀ can never leave its bounds.

rw_sd is on the natural scale. --rw-sd "R0=5" means “perturb R₀ by about 5 units per step.” Internally, this is converted to the transformed scale via the delta method: for log-transformed params, the effective SD on log scale ≈ rw_sd / current_value. For R₀=56.8 with rw_sd=5, the perturbation is ~9% per step on the natural scale.

Scale warnings: If rw_sd is >50% of the parameter value, the perturbation is dangerously large. If <0.1%, the parameter isn’t exploring. The CLI warns with suggested adjustments.

Multi-chain and chain-agreement Â

Run multiple independent IF2 chains from different random seeds to detect multimodality and assess convergence. An IF2 fit is a fit whose [method] is algorithm = "if2":

# fit.toml
[method]
algorithm = "if2"
backend = "chain_binomial"
chains = 4
particles = 1000
iterations = 50
cooling = 0.7

camdl fit run fit.toml --seed 42

Chain-agreement  measures across-chain agreement (Gelman–Rubin 1992 form, applied to IF2’s per-iteration parameter-mean trajectory across chains; this is not a posterior mixing statistic — IF2 is an MLE optimizer, not a sampler, so  here measures whether the optimizer’s chains agreed on a basin, not whether a posterior has mixed). Computed from the last half of iterations:

• Â < 1.1: converged (✓) — chains agree
• Â 1.1–1.5: uncertain (~) — might need more iterations
• Â > 1.5: not converged (✗) — surface may be multimodal

Note: Bayesian (PGAS, PMMH) outputs continue to use the name rhat for their own posterior-mixing diagnostics; only IF2 uses chain_agreement / Â.

Two IF2 files: scout, then refine

A common MLE workflow is a cheap exploratory IF2 file followed by a sharper one that starts from its result. Each is a fit.toml with the same problem half and its own [method]; the second names the first by handle. Scout and refine are a convention, and the knobs (chains, particles, iterations, cooling) are starting points to adapt, not defaults the tool enforces:

Scout — 8 chains, 500 particles, 30 iterations, cooling = 0.70 (mild). Exploration: chains stay hot enough to wander across basins rather than quenching onto the first local optimum. Over the 30 iterations the perturbation SD shrinks only from 1.0× to 0.49× initial. Run this first to find problems: Is the surface multimodal? Which parameters are identifiable? Is the observation model appropriate? The cross-chain  at the end, combined with the loglik-eval decibans-spread gate (see camdl-inference-spec §6.1.1), is the multi-modality diagnostic.

# scout.toml
[method]
algorithm = "if2"
backend = "chain_binomial"
chains = 8
particles = 500
iterations = 30
cooling = 0.70

Refine — 4 chains, 1000 particles, 50 iterations, cooling = 0.05 (aggressive), starts = { from_mle = "@scout" }. Starts every chain at the scout’s best-chain parameters and collapses them tightly onto the local MLE — final SD is 0.25% of initial, so particle clouds concentrate near the scout’s endpoint. Â is then a statement about IF2’s own perturbation re-spreading the swarm, not about independent starts; the summary’s seeded from line says so. A scout whose verdict is not converged is refused as a source unless --allow-nonconverged-source records the choice.

# refine.toml
[method]
algorithm = "if2"
backend = "chain_binomial"
chains = 4
particles = 1000
iterations = 50
cooling = 0.05
starts = { from_mle = "@scout" }

Cooling follows pomp’s cooling.fraction.50: the perturbation SD reaches cooling × rw_sd after cooling_target_iters iterations (default 50, pomp’s fixed value) and keeps shrinking geometrically after that. Formula and worked numbers: docs/methods/cooling.md.

Initial-state parameters (perturb_only_at_t0)

Initial conditions (S₀, E₀, I₀) set the starting state but don’t change during simulation, so they are perturbed only at t=0, not at every observation. Declare that with perturb_only_at_t0 = true on the [estimate] entry. They are model-declared parameters: list them under [estimate] in the fit.toml like any other parameter, and the fit estimates them. The filter jitters their initial values once when particles initialize, then holds them fixed as it runs forward — whereas a parameter like R₀ is perturbed at every observation time. The flag is IF2’s schedule, so only an if2 config accepts it; any other method refuses the config at load rather than run without it.

The flag changes how IF2 perturbs the parameter, not how its convergence is judged: its tail chain agreement  is gated like every other estimate’s, and fit summary reports the dynamics-only sub-verdict beside the overall one.

Under IF2 that single t=0 jitter is the only channel through which the data can inform an initial-condition parameter, so each particle then builds its own initial compartment counts from its own jittered value: particle j starts at x_0^{(j)} = \text{init}(\theta^{(j)}), not at the swarm mean’s initial state. That is what puts spread in the swarm at t = 0, and it is why the first reweight can discriminate between particles — the precondition ic_free = true relies on. A parameter that appears in initial_conditions and in no rate and in no observation model reaches the likelihood through this path and no other.

The other way to put spread in the swarm at t = 0 is to declare the initial condition as a law in the model (init { I ~ poisson(rate = I0) }), which every filter draws per particle. That one works under pfilter and pmmh too, where \theta is shared by the whole swarm and the jitter above does not exist. See “Initial-state parameters” below.

---

Profile likelihoods

Fix a focal parameter at a grid of values, run IF2 at each to maximize over the remaining parameters. The resulting curve shows how the MLE changes — revealing identifiability, confidence intervals, and parameter correlations.

1D profile

camdl profile model.camdl --starts from_params=p.toml --data cases.tsv \
    --sweep "R0=lin(10,80,8)" \
    --rw-sd "sigma=0.01,gamma=0.01" \
    --particles 500 --iterations 30 --n-starts 3 --parallel 8

Output: TSV with R₀, max loglik at each grid point, and the estimated values of all other parameters.

A sharp peak means R₀ is well-identified. A flat profile means R₀ is not identifiable from the data (the model fits equally well across a range of R₀ values).

2D profile

camdl profile model.camdl --starts from_params=p.toml --data cases.tsv \
    --sweep "alpha=0.85,0.90,0.95,0.99" \
    --sweep "gamma=0.06,0.08,0.10,0.12" \
    --rw-sd "R0=2,sigma=0.01" \
    --particles 500 --n-starts 2 --parallel 8

Shows ridges and correlations between parameters. An elongated contour along the alpha-gamma diagonal means those parameters trade off — you can’t identify both independently.

Priors and precedence

camdl fit run resolves priors for each estimated parameter via a three-tier precedence chain: fit-toml > model-IR ~ syntax > flat. camdl profile resolves the same chain but only records the result: every profile cell is a maximization that reads no prior. The behaviour at tier 3 differs between the two subcommands — see below.

Tier 1 — fit-toml priors (highest). The fit toml’s [estimate.<param>.prior] block wins over any other source. For fit run the fit toml is always loaded; for profile it’s the --fit <toml> flag.

[estimate]
beta = { bounds = [0.01, 5.0], prior = { log_normal = { mu = -0.3, sigma = 0.5 } } }

Tier 2 — model-IR ~ priors (fallback). When the fit toml doesn’t declare a prior for an estimated parameter, the resolver falls through to whatever the .camdl file declared via ~ syntax.

parameters {
  beta : rate in [0.001, 5.0] ~ log_normal(mu = -0.3, sigma = 0.5)
}

This is the recommended source of truth for stable priors: the model file is the single artifact reviewers read for the structural priors, and individual fit tomls only override when doing sensitivity analysis. Stripping prior duplicates out of N fit tomls and into one model file is what gh#75 unblocked.

Tier 3 — Prior::Flat (last resort). Improper uniform. The behaviour at this tier is asymmetric across subcommands:

• camdl profile: tier 3 is recorded as flat_fallback in the profile-base sidecar, with no warning. The cells (if2, nl-*) maximize the likelihood, so the prior changes nothing they compute; only --starts from_prior reads it, to draw starting points.

• camdl fit run: tier 3 is a hard error at config-load time (the fit refuses to start). The downstream interpretation of a fit-run chain — canonical posterior in fit_summary.json, consumed by tooling that treats those samples authoritatively — is too high-stakes for a silent demotion to “scaled likelihood”. Users who genuinely want flat priors declare it explicitly via

  [estimate.beta]
  prior = { flat = {} }

  This opt-in path is fully accountable: the toml records the intent, run.json’s resolved_priors records the source as flat_explicit, and no warning fires (the user said what they meant). Silent fallback to flat is unreachable.

When all three tiers fail (no fit-toml prior, no model ~, no explicit-flat opt-in), camdl fit run errors with a 2-column table naming every offending parameter plus all three remedies:

error: [method] (algorithm = "pmmh") has parameters with no resolved prior:

  beta        no prior in fit toml, no `~` in model file
  gamma       no prior in fit toml, no `~` in model file

To proceed, do one of:

  (i)   Declare `prior = { <dist> = { ... } }` in the fit toml's
        [estimate.<param>] for each listed parameter.
  (ii)  Declare a `~ <dist>(...)` prior in the model file for
        each listed parameter.
  (iii) Opt into flat priors explicitly via
        `prior = { flat = {} }` in the fit toml — only do this if you
        intentionally want the chain to target the unconditioned
        likelihood (scaled-likelihood posterior).

# Bayesian fit with priors in the model file (no duplication in N fits)
camdl fit run fits/synth.toml --seed 0

Precedence rules for parameter values (the unified chain shipped in the 2026-05-25 CLI UX revision; see docs/camdl-run-spec.md §1.3 for the authoritative tier list and docs/dev/proposals/archive/post-alpha/2026-05-25-cli-init-and-params-ux.md §“Precedence (last wins)” for the design rationale). Each tier overrides the tier above it:

1. Compiled-in value. The value already in the compiled IR. A .camdl never writes one (parameters {} declares names, kinds, bounds and priors), so this tier is empty unless camdlc --set NAME=V pinned the parameter at compile time.
2. fit.toml [fixed] block. When --fit is in scope, every key in [fixed] overrides the compiled-in value.
3. --fixed-file <toml> (repeatable). A flat parameter TOML; top-level keys are parameter names. Layered in declaration order — later files override earlier ones.
4. Scenario preset (--scenario NAME or composed --enable/--disable). The active scenario’s params.set / params.scale directives override everything above. Scenarios travel with the model and beat user-supplied params files by design — choosing a scenario is choosing the model author’s documented bundle.
5. --fixed NAME=VALUE (repeatable, highest). The per-invocation override; wins over every value the scenario does not set or scale. Use this for “I want to change one value for this one run.” A --fixed on a parameter the scenario sets or scales is refused (gh#571): it would replace the scenario’s value while the run is still labelled with the scenario. Put the value in the model’s scenarios {} block, or, for a scale, in a --fixed-file the scale applies over.

The legacy --params and --param flags on inference subcommands (profile, if2, fit run, survey) were removed in the same revision. Their replacements are:

• --fixed-file <toml> for the “load many values from a file” case.
• --fixed NAME=VALUE for the “change one value” case.
• --starts from_params=<toml>, on fit run and on profile alike, for the warm-start chain origin case — when the file is a starting point for inference, not a pin.

--fixed/--fixed-file on inference subcommands also removes the listed parameter from the [estimate] set if it was there — so --fixed gamma=0.1 --sweep tau=lin(-35,-1,30) is the canonical slice-while-holding-gamma pattern. Each such parameter gets one info: line on stderr, so the override is never silent: info: gamma: `--fixed gamma=0.1` pins it, so it is not estimated (the fit config lists it in [estimate]).

Override visibility. Whenever one tier replaces, with a different value, a value a lower tier supplied — --fixed over the active scenario, a later --fixed-file over an earlier one — the resolver prints one line to stderr naming the parameter, both values and both sources, and records every replacement in run.json’s parameters_provenance block:

"beta": {
  "value": 0.5,
  "source": "fixed_cli",
  "role": "fixed",
  "overrides": [
    {
      "replaced": { "source": "scenario", "value": 0.3, "scenario": "worst_case" },
      "by": { "source": "fixed_cli", "value": 0.5 }
    }
  ]
}

The line is info: for a replacement the command’s own flags asked for, and warning: when the replaced value is a fit’s estimate (a params file carrying a fit [provenance] block). CLI override of a scenario value is a legitimate quick-test workflow; the line and the provenance record ensure six-months-later auditing of which value actually ran isn’t blocked by archaeology. The full rule is docs/camdl-run-spec.md §1.3.

Profile focal parameters. The focal swept parameter(s) are always removed from the estimated set, even when they appear in the fit toml’s [estimate]. Listing a focal parameter in --fixed/--fixed-file (or in the fit toml’s [fixed]) is a hard error — a parameter cannot be simultaneously swept and pinned. A scenario that sets or scales a focal parameter is replaced by the grid at every point, and the profile says so once per run: info: scenario `baseline` sets beta = 0.4; beta is profiled, so the grid values replace it.

Provenance. Both subcommands write per-parameter resolved_priors into run.json:

"resolved_priors": [
  { "param": "beta",  "source": "fit_toml" },
  { "param": "gamma", "source": "model_ir" },
  { "param": "sigma", "source": "flat_explicit" }
]

The four wire values are "fit_toml", "model_ir", "flat_explicit" (gh#75 — fit run only), and "flat_fallback" (profile only — the silent-fallback case). Reviewers reading a fit dir’s run.json can audit at a glance whether the chain targeted a posterior with priors or the unconditioned likelihood.

The fit CAS identity includes the model IR (the model digest in FitDigest, via cas::fit_level_hash), so re-running against the same fit toml after editing a ~ prior in the model file produces a different cache dir. For profile, the CAS hash additionally keys on fit_toml_hash + resolved per-parameter prior sources so re-running against the same model with a different --fit flag produces a different umbrella.

Reparameterization for better sampler geometry

The parameterization a model is written in — not just its priors — determines how hard the posterior is to sample. camdl assigns the inference-scale transform from a parameter’s constraint (positive → log, probability → logit, with the Jacobian handled), exactly as Stan does; it does not second-guess your parameterization by role. So when a parameter’s natural scale is a poor sampling geometry, the fix is to reparameterize the model, not to expect a magic transform.

Overdispersion parameters are the classic case. For a neg_binomial(mean, r) or beta_binomial(mean, concentration) model, the dispersion parameter (r, concentration — call it φ) has a treacherous scale: the likelihood is flat as φ → ∞ (the Poisson/binomial no-overdispersion limit is a plateau), so a generic prior and a boundary-avoiding init put mass on, and strand chains in, the near-no-overdispersion region. This is a documented failure mode (Stan Prior-Choice Recommendations: “putting priors on the over-dispersion parameter of the negative binomial” fails because “most of the prior mass is on models with a large amount of over-dispersion”).

The fix (Stan’s recommendation): put the parameter on the reciprocal scale. The generic prior “works much much better on the parameter 1/φ. Even better, you can use 1/√φ” — because 1/φ → 0 is the no-overdispersion limit, so a half-normal / exponential prior there has mass exactly where the data usually lives and regularizes the plateau. In camdl, declare the reciprocal as the estimated parameter and derive φ from it:

parameters {
  phi_inv : positive   # the SAMPLED parameter, prior lives here
    ~ half_normal(sigma = 1)
}
# ... where the concentration is used:
#   ... concentration = 1.0 / phi_inv

Then phi_inv → 0 is the binomial limit, the init lands near it, and a gradient or adaptive-MH sampler moves off the plateau. This is a modeling choice you make — camdl (like Stan) gives you the constraint transform and the prior machinery; the reciprocal parameterization is yours to apply where the geometry calls for it.

Per-cell diagnostics

Profile output gains a fixed-schema block of per-cell convergence columns appended after the existing focal / loglik / parameter columns (gh#74 Option B). The columns are written into the umbrella summary.tsv (the file --output mirrors) and into each per-seed profile.tsv.

Read by column name, not column index. The diagnostic columns are the API; their order is stable across runs, but defensive consumers should look them up by name so future schema additions land cleanly.

The full list:

Column	Algorithm	Meaning
loglik_spread_starts	all	max − min of per-start final log-likelihoods. > ~5 nats means the starts disagree on the basin.
loglik_rhat_starts	IF2	Gelman–Rubin R̂ on the per-start log-likelihood traces (Gelman & Rubin 1992, Statist. Sci. 7(4); Brooks & Gelman 1998 corrected variant). > 1.05 is the conventional “chains haven’t agreed yet” threshold.
starts_n_completed	all	Count of starts that produced a finite final loglik (didn’t diverge).
iterations_used	IF2	Final cooling-step index reached (= --iterations on normal completion).
cooling_final	IF2	Mean across estimated parameters of the final iteration’s effective_rw_sd — the actual ending perturbation SD, not the target.

For algorithms that don’t supply a given column (e.g. iterations_used on an nl-* run), the cell renders as NaN (capital N — camdl’s TSV NaN convention).

The K<3 rule. loglik_rhat_starts is NaN when fewer than three of the K starts have a usable trace. Gelman–Rubin R̂ is undefined at K=1 and unstable at K=2; the rule prevents a spurious diagnostic from a single-chain spike. To get a finite R̂ supply --n-starts 3 or more.

The diagnostic R̂ is computed on log-likelihood traces, not on per-parameter posteriors, so it is a chain-level “are these starts walking the same basin” check rather than a full multivariate convergence diagnostic. The column name reflects that: it’s loglik_rhat_starts, not rhat.

Cross-seed aggregation. For multi-seed profile runs (--seeds 1,2,3), each cell’s diagnostic columns are averaged across seeds in summary.tsv (with starts_n_completed summing rather than averaging — it’s a total count). Per-seed values remain visible in each replicates/seed_<n>/profile.tsv.

---

PGAS (Particle Gibbs with Ancestor Sampling)

IF2 finds the MLE. PGAS characterizes the full posterior — credible intervals, parameter correlations, posterior trajectory samples.

How it works

PGAS is a Gibbs sampler alternating two steps per sweep:

Step 1: θ | X, y (parameter update). With the full latent trajectory X known, the complete-data log-likelihood is exact:

\log p(y, X \mid \theta) = \sum_s \log p(x_{s+1} \mid x_s, \theta) + \sum_t \log p(y_t \mid x_t, \theta)

No particle filter, no estimation noise. The transition density at each substep is a product of Binomial log-PMFs mirroring the Euler-multinomial decomposition in the simulation. Parameters are proposed via NUTS (gradient-based) or one-at-a-time MH.

Step 2: X | θ, y (trajectory update). CSMC-AS (Conditional SMC with Ancestor Sampling) produces a new trajectory sample from p(X \mid \theta, y). One particle slot is clamped to the reference trajectory; ancestor sampling at each substep reconnects the reference to the free-particle cloud via the transition density. Trajectory renewal (fraction of the traceback from non-reference particles) measures CSMC health — near 0% means path degeneracy, above 50% means healthy mixing.

NUTS gradient proposals

The complete-data log-likelihood is differentiable with respect to parameters: the Binomial log-PMF depends on rates via p = 1 - \exp(-\text{rate} \cdot dt), and the rates are differentiable expressions from the model.

The OCaml compiler performs source-to-source symbolic differentiation of rate expressions (autodiff.ml), emitting rate_grad fields in the IR JSON. The Rust backend evaluates these derivative expressions via the same eval_expr interpreter — no runtime autodiff, no finite differences.

NUTS (No-U-Turn Sampler, Hoffman & Gelman 2014) uses these gradients to propose all parameters jointly via Hamiltonian dynamics. A two-phase warmup adapts both the step size (dual averaging) and the diagonal mass matrix (empirical variance from burn-in). The mass matrix rescales each parameter by its posterior variance, so NUTS takes appropriately-sized steps in every direction.

Running PGAS

# From the priors (the default whenever every estimated parameter has one):
camdl fit run fit.toml --seed 42

# From a stored fit's posterior, one draw per chain — R̂ stays meaningful:
#   [method] starts = { from_posterior = "@short" }
camdl fit run fit.toml

# Every chain at a stored point estimate — R̂ is then not assessed:
camdl fit run fit.toml --starts from_mle=@mle

# Force MH-within-Gibbs instead of NUTS
camdl fit run fit.toml --no-nuts

Configuration in fit.toml:

[method]
algorithm = "pgas"
chains = 4
sweeps = 10000
particles = 100
burn_in = 2000
thin = 5
n_trajectories = 200 # posterior trajectory samples per chain

Output per chain: trace.tsv (parameters + log-likelihood per sweep) and trajectories.tsv (posterior latent state draws). The latter is a tidy/long file — all saved draws stacked, with leading chain draw time id columns, then the compartment columns, flow_<transition> columns, and inc_<stream> columns (the model’s incidence projection of the latent flows, so posterior incidence needs no finite-differencing). A sibling trajectories.json manifest records the method, time granularity, column list, and the conditioned: true flag (these are smoother paths X | θ, y, not the free-forward posterior-predictive a simulate --obs run produces). Load a posterior band in two lines:

import pandas as pd
df = pd.read_csv("chain_0/trajectories.tsv", sep="\t", comment="#")
band = df.groupby("time")[["S", "E", "I", "R"]].quantile([0.05, 0.5, 0.95])

Initial-state parameters

A parameter that only sets the initial state (an initial infected count I0, an initial susceptible fraction) reaches the complete-data log-likelihood through one channel and no other: the trajectory’s initial state is stored, not recomputed, so under PGAS the transition and observation terms are invariant to it. For the data to inform it, the initial state has to be a random variable.

That is a declaration in the model, not something the sampler infers:

init { I ~ poisson(rate = I0)
       S = N0 - I }

With it, PGAS’s CSMC draws each particle’s own x_0 and \log p(x_0 \mid \theta) enters the complete-data likelihood as its own term — reported per sweep in the initial_state_ll column of the chain trace — so the MH ratio constrains I0. IF2, pfilter and pmmh use the same declaration for the draw alone. Without it the term is exactly zero, on every chain and from every start.

Nothing auto-detects initial-value parameters. PGAS used to: it finite- differenced a rounded initial count from the chain’s own starting draw and attached a Binomial density to whatever compartment the parameter moved, so two chains of one fit could disagree about whether the term existed at all. Both the detector and that density were removed (gh#719). The perturb_only_at_t0 flag is a separate thing — an IF2 perturbation schedule, refused at config load under every other method.

Why a deterministic initial condition identifies its parameter poorly

This is the part that surprises people, and it is worth understanding before reading a posterior for an initial-state parameter.

An integer compartment’s initial value is rounded. Written I = I0, the compartment starts at round(I0), so I0 = 50.0 and I0 = 50.4 produce the same initial state, the same trajectory density, and the same likelihood. The likelihood is a step function of I0: its derivative is zero almost everywhere, and non-existent at the boundaries.

Two consequences follow, and both show up in a posterior:

• A gradient-based move learns nothing about it. NUTS sees a flat direction, so I0 moves only when a proposal happens to cross a rounding boundary. Its marginal ends up close to its prior — not because the data are uninformative about the seed, but because the model cannot express the difference.
• Its correlations with other parameters are hidden. Any real coupling — between the seed and the growth rate, say, or the reporting fraction — is smeared across each rounding bin, because within a bin the data cannot distinguish one value of I0 from another. Pair plots look less structured than the model actually is.

Writing I ~ poisson(rate = I0) removes both. The rate is continuous while the draw is integral, so log p(x_0 \mid \theta) is smooth in I0; the parameter has a real gradient everywhere and the sampler explores (x_0, θ) jointly. Correlations that were previously invisible become visible. They were always in the model — the rounding was hiding them.

Expect the drawn parameterisation to give a wider marginal for I0, not just a better-resolved one: the law puts a layer between I0 and the realised seed, so I0 is one step further from the data. A narrower posterior after the change is worth understanding before trusting it.

A second, smaller effect: writing the complement as a coefficient (S = N0 - 3.386 * I0) does not survive rounding, because round(a·I0) + round(b·I0) ≠ round((a+b)·I0). Such a model carries a small population imbalance that varies with I0. Writing S = N0 - I - E1 - E2 - E3 against the drawn values makes the budget exact.

The ~ form arrived in d4e89659 (2026-08-24), which also removed the auto-detector described above; before it, every initial condition was deterministic-and-rounded and behaved as this section describes.

Spatial models and seeding (iota)

Spatial models with inter-patch coupling need care to ensure inference works correctly. Two issues arise that don’t affect single-patch models:

Seeding terms. If the infection rate for patch i is \beta \cdot S_i \cdot I_i / N_i, it goes to exactly zero when I_i = 0. The stochastic simulator can still draw events from near-zero floating-point rates (importation coupling creates tiny nonzero values), but the density evaluator computes the rate as exactly zero and rejects the trajectory.

Fix: add a small seeding term to the infection rate: \beta \cdot S_i \cdot (I_i + \iota) / N_i where \iota \approx 10^{-6}. This ensures the infection rate is never exactly zero, allowing importation-driven infections to have finite (though very small) density. pomp spatial models use the same pattern.

Not all spatial models need iota. Models with constant importation via events {} blocks, or models where the rate expression already includes an additive term, are fine without it.

Discrete-event seeding and PGAS

A model with discrete seeding via events { ... add(E, n_seed) at [tau] } runs under PGAS (gh#80). The chain-binomial simulator and the density evaluator are consistent at the event substep: at substep s_{\text{event}}, counts_before is the pre-event state, all stochastic flows are drawn from pre-event rates (which are zero for the seeded compartment’s downstream transitions), and the event delta is applied atomically with transition deltas at end-of-substep. The transition log-density at that substep is exactly zero (finite). No Gaussian-pulse workaround is needed.

Use PMMH to fit the seed timing or size. PGAS scores the recorded latent path, and the channel through which a parameter enters an effect — its amount (n_seed), its time (tau) — or a plain init {} assignment (init { I = I0 }) appears in no term of that score, so the sampler cannot move the parameter through it (a parameter entering only that way keeps its prior). PGAS refuses such a parameter by name, and IF2 refuses an estimated parameter in an effect’s schedule (it resolves the schedule once per run). PMMH re-simulates the path at each proposed θ and fits all three. For an initial size, declaring a law (I ~ poisson(rate = I0)) makes it estimable under PGAS. (A replay θ-move that would score these under PGAS is gh#986.)

With a fixed tau outside the simulation window, the seed never fires, predicted incidence stays zero, and the observation density goes to -\infty against real data with nonzero cases. camdl prints a split-by-component diagnostic at startup so the failure mode is identifiable (transition_ll vs observation_ll vs ivp_ll).

During CSMC ancestor sampling on event-using models, the density evaluator returns -\infty for a free particle whose pre-step state has zero rate for a transition that fired in the reference’s flow record. This is correct: the free particle cannot be the ancestor and is excluded from the categorical.

Time step size. The Euler-multinomial approximation assumes exit probabilities are small per substep. In spatial models with high R_0 and dt = 1, p_{\text{total}} can approach 1, causing overdrafts where total withdrawals from a compartment exceed its population (resolved by clamping). Use a smaller dt (e.g., 0.25) to keep p_{\text{total}} < 0.3 and avoid approximation breakdown.

MCMC initialization strategy

PGAS chains should be initialized at or near a known high-likelihood region, not from random or diffuse starting points. The recommended workflow:

1. IF2 scout: Run 8–16 chains with random starts to map the likelihood basins. More chains are needed for spatial models where the surface is multimodal (R0–sigma–amplitude ridges create multiple basins).
2. Profile likelihood: Run a 1D profile over R0 (the parameter most prone to basin structure) to confirm which basin has the highest likelihood.
3. Initialize PGAS: Start all chains at the best IF2 MLE ± small jitter (e.g., ±5% per parameter). This avoids wasting burn-in searching for a basin that IF2 already found.

Starting chains near the mode is standard MCMC practice (Gelman et al., BDA3; Stan’s default workflow optimizes first, then samples). MCMC convergence guarantees are asymptotic — initialization affects only burn-in length, not the target distribution. Starting from a good point reduces wasted computation; it does not bias the posterior.

When initialization matters most: Spatial models with seasonal forcing. The R0–sigma trade-off creates basins separated by 50+ log-likelihood units. IF2 with only 4 chains can land in the wrong basin (e.g., R0≈28 instead of the true R0≈20), and PGAS initialized there may never cross the barrier. More IF2 scout chains is the fix.

Where priors live, and what starts = "from_prior" reads. Parameter priors belong in the model, declared with ~: param : <kind> in [<lo>, <hi>] ~ <dist>(...); a fit.toml [estimate].<param>.prior overrides the model’s declaration for one fit. from_prior draws each chain’s start from the prior the sampler scores against — the same precedence, fit.toml over model — so the start and the target cannot disagree. A parameter whose resolved prior is flat (or hierarchical, with no hyperparameter values at start time) has no distribution to draw from and falls back to a bounds-uniform start, with a single startup warning: line naming it; declare a proper prior in either place to draw it from one. The same fallback applies to starts = { from_posterior = … } for parameters absent from the posterior source.

---

Fitting on the ODE backend

Everything above targets the stochastic-process likelihood p(y_{1:T} \mid \theta) — the integral over latent trajectories that the particle filter estimates (PF, IF2, PMMH) or that PGAS conditions on. The ODE backend offers a different object. Integrate the mean-field system (see camdl docs backends, the ODE section) and there is exactly one trajectory per \theta; the latent-path integral collapses, and the likelihood is read off directly:

p(y_{1:T} \mid \theta, \text{ODE}) = \prod_t p\bigl(y_t \mid \pi_t(x(\theta)), \theta\bigr)

where \pi_t(x(\theta)) is the model’s projection (incidence or prevalence) of the integrated state at observation t. There is no Monte Carlo and no filter variance — the same observation model (neg_binomial, discretized-normal, …) is scored against the ODE’s projected value instead of a particle cloud’s. This deterministic marginal likelihood is a genuinely different statistical object from p(y \mid \theta): the two coincide only in the low-noise / large-population limit, where the mean field is a faithful stand-in for the stochastic process. Choosing backend = "ode" for a fit is choosing to target this object — appropriate for equilibrium or large-population models where demographic stochasticity is negligible and running a particle filter would be structurally redundant.

Two things can be asked of this likelihood: where its maximum is, and what the posterior around it looks like. Three NLopt optimizers answer the first and the mh sampler answers the second; all four are [beta].

Maximum likelihood: nl-sbplx, nl-bobyqa, nl-lbfgs. Each runs a local NLopt algorithm on p(y | θ, ODE) inside the [estimate] box, from as many independent starts as chains says, and reports the best. They differ only in how they search. nl-sbplx (a Nelder–Mead variant) and nl-bobyqa (a quadratic trust region) are derivative-free: they learn the local shape by evaluating the likelihood at neighbouring points, so one objective evaluation is one ODE solve, and a search in d dimensions needs many of them before it can take a confident step. nl-lbfgs is gradient-based: each evaluation integrates the augmented (x, ∂x/∂θ) system once — forward sensitivities, carrying \partial x / \partial \theta alongside the state through fixed-step RK4, with the rate and observation derivatives differentiated source-to-source by the compiler (the same machinery behind PGAS’s rate_grad; see “NUTS gradient proposals” above) — and returns the log-likelihood together with its exact gradient, so the direction to move in is read off rather than sampled for.

The trade is one of cost per step against number of steps. A gradient evaluation costs more than a value evaluation, and more as the number of estimated parameters grows, because the sensitivity block carries one column per parameter; but it replaces the dozens of probing evaluations a derivative-free method spends discovering the same direction. On a five-parameter SEIR fitted to weekly incidence, nl-sbplx reached its own convergence in 3609 objective evaluations and 2.9 s where nl-lbfgs took 94 and 0.22 s — the same optimum, the same log-likelihood to six decimals. On a two-parameter SIR the counts were 216 and 14. Start with nl-sbplx, which fits anything the ODE backend can integrate; reach for nl-lbfgs when the search itself is the bottleneck and the model is differentiable.

Bayesian posteriors: mh. Adaptive Metropolis-Hastings on the deterministic marginal. It needs no gradient, so it runs on any ODE model the backend can integrate, and it is the robust default for identifiable or low-dimensional posteriors. It will diagnose a difficult posterior — a weakly-identified ridge shows up as low effective sample size — but random-walk MH mixes poorly through a curved ridge: adaptive covariance corrects a rotated Gaussian, not a banana. camdl has no gradient-based sampler on the ODE backend; on a stochastic backend the marginal likelihood is an intractable integral with no closed-form gradient, and the gradient sampler is the NUTS move inside pgas (use_nuts, on by default), conditioned on a sampled trajectory. Requesting an invalid pair — say pgas on ode — fails at config load with a message naming the right alternative, and algorithm = "nuts" is refused on every backend with a message pointing at pgas and mh.

The differentiability requirement

nl-lbfgs needs a closed-form gradient of the whole likelihood, so it accepts only a differentiable model. One capability check answers for it, and it refuses, with a message naming the reason, a model that has:

• a rate or observation term whose gradient the compiler cannot emit (a live-but-undifferentiated coefficient — see the gradient-refusal list in camdl docs features);
• an adaptive integrator (rk45) — its step sequence is data-dependent, so the sensitivities have no fixed step to propagate through; use fixed-step rk4;
• a scheduled effect — an interventions {} or events {} entry, a discrete state jump the smooth sensitivity flow cannot cross;
• an initial condition the ∂init/∂θ seed cannot express — one drawn from a distribution, or a parameterized one whose expression camdl emits no initial-condition gradient for.

Any of these is a hard error that points at the fix, raised at config validation — before a leaf is claimed, so a refusal costs a message rather than a half-written leaf. A model that genuinely needs one of them — a mid-run campaign, an estimated S_0 — is fit either with the gradient-free ODE methods (mh for a posterior, nl-sbplx for an MLE), or with the stochastic-process methods, which handle scheduled effects and estimated initial states.

The gate answers about the model the run will integrate, not the model as written. interventions { } is toggleable and defaults off, so a fit that declares a campaign and never enables it is not refused for a scheduled effect its trajectory does not contain.

A refusal that survives to the search is possible in one more way: the gradient may be finite everywhere the gate can check and still not be computable at some θ the optimizer reaches (a solve that diverges there). nl-lbfgs fails that chain with the reason and the θ rather than scoring the point badly and carrying on, because a gradient method given no derivative has no direction to search in — unlike the derivative-free methods, which treat such a point as merely very bad and steer away from it.

Running an ODE fit

These are methods like any other: pick the algorithm and set backend = "ode". The Bayesian ones additionally need a prior on every estimated parameter (see “Priors and precedence”); the MLE ones do not read priors at all.

The MLE methods take chains (independent multi-start optimizations), tolerance (NLopt’s xtol_rel) and max_evals (per-chain evaluation budget), and no prior:

[method]
algorithm = "nl-lbfgs" # gradient MLE; needs a differentiable model
backend = "ode"
chains = 4
starts = "uniform_unconstrained" # a spread rule; `single` collapses to one chain

Gradient-free mh takes iterations (total MCMC steps) with an optional burn_in:

[method]
algorithm = "mh"
backend = "ode"
chains = 4
iterations = 4000
burn_in = 1000

camdl fit run fit.toml

camdl fit methods prints the live registry — the one-liner, the “use for”, and the beta caveat for every (algorithm, backend) pair. The two backends compute different objects (chain_binomial → p(y|θ); ode → p(y|θ, ODE skeleton)), so read that guidance before choosing which likelihood a fit should target.

---

Diagnostic interpretation guide

Healthy pfilter trace

time  ll_increment  ESS    pred_mean  pred_q05  pred_q95  observed
7     -4.2          2800   45         12        95        52      ← data in interval
14    -3.8          3100   120        48        220       135     ← data in interval

ESS stays above 50% of N. Data falls within prediction interval. Log-likelihood increments are moderate (not extreme).

Degenerating filter

time  ll_increment  ESS    pred_mean  pred_q05  pred_q95  observed
7     -4.2          2800   45         12        95        52
14    -12.8         23     120        105       140       350     ← data far outside

ESS crashes to <1% of N. The data is very surprising given the model’s predictions. Causes: wrong parameters, wrong observation model, missing model features (e.g., no seasonal forcing when the data has seasonal epidemics).

IF2 convergence trace

iteration  loglik   R0      gamma
0          -6200    42.3    0.15     ← exploring
5          -4100    51.2    0.09     ← approaching
15         -3850    55.8    0.084    ← converging
30         -3810    56.5    0.083    ← stabilizing
50         -3805    56.8    0.083    ← converged

Log-likelihood should improve monotonically (with noise). Parameters should approach stable values. If loglik oscillates without improving, rw_sd is too large. If parameters haven’t moved after 20 iterations, rw_sd is too small.

IF2 chain-agreement diagnostics

 (across 4 chains, last 25 iterations):
  R0           Â=1.02 ✓ range=[55.2, 58.1]
  sigma        Â=1.01 ✓ range=[0.078, 0.080]
  gamma        Â=3.20 ✗ range=[0.065, 0.120]

R₀ and sigma have converged (Â < 1.1, tight range). Gamma has not (Â=3.2, wide range). This means gamma is either poorly identified or the surface is multimodal along the gamma axis. Run a profile likelihood for gamma to distinguish.

---

The fit workflow

The low-level commands (camdl pfilter, camdl profile) are building blocks. IF2 is not a standalone command — an IF2 fit is a fit whose [method] is algorithm = "if2". For all model fitting, camdl fit provides a structured workflow driven by a fit.toml configuration file, one problem and one method per file:

Two files that share a problem half share the fit-level hash (<h8>); each is its own segment, named by its stem, with its own --label. A method reads a stored fit only through a starts handle, and every multi-chain sampler writes chain_starts.tsv, the record of where each chain actually began. An IF2 leaf’s mle_params.toml is a standard params file with provenance hashing that feeds directly into camdl simulate.

from_mle reads the scout’s point estimate from its fit_state.toml and starts every chain there, so the posterior fit’s R̂ is not assessed: R̂ compares chains that began apart, and chains that all begin at one point can agree without having explored the posterior. An IF2 leaf writes no draws.tsv, so from_posterior = "@scout" is refused. To keep R̂, seed from a short posterior run instead (starts = { from_posterior = "@short" }, above).

camdl fit run    scout.toml --seed 1 --label scout
camdl fit run    posterior.toml --seed 1
camdl fit summary @scout
camdl fit summary results/fits/<dir>/

When start = is omitted in [estimate]

Each [estimate.X] entry’s start = field is optional. When you omit it (and the model file doesn’t already declare a value for the parameter via parameters { X : rate = 0.3 } or a scenario), the runner draws a single Transform-aware uniform value within the parameter’s bounds and uses that as the base point. From there the selected init mode perturbs per-chain as usual.

The draw is:

• Log-uniform for Log-typed parameters with strictly positive bounds (rates, positive quantities) — equivalent to drawing uniformly in [ln(lo), ln(hi)] and exponentiating, so a parameter with bounds [1e-6, 1.0] doesn’t collapse to a value near 0.5.
• Linear-uniform for Logit/None parameters and for any parameter whose bounds aren’t strictly positive.

The draw is deterministic per (seed, parameter_name): re-running with the same --seed gives the same fallback values, and two parameters with identical bounds at the same seed get different values (their names hash differently). Different seeds give different fallback points within bounds — useful when you want a seed sweep to also sweep over starting positions for the unspecified parameters.

This replaces an earlier bounds-midpoint heuristic that gave the same point at every seed and ignored the parameter’s transform.

Where the chains begin: starts

How chain (or per-cell) starting points are drawn. Set under [method] with the starts key, or overridden on the CLI with fit run --starts <spec>; camdl profile --starts <spec> takes the same rules for its per-cell starts, and --n-starts <N> says how many chains are drawn from them. Honoured by every method: IF2, PGAS, PMMH, mh, NLopt (nl-sbplx, nl-bobyqa, nl-lbfgs), and profile.

[method]
algorithm = "if2"
backend = "chain_binomial"
chains = 16
starts = "uniform_unconstrained" # every parameter has a prior? then the default is from_prior

The rules fall into two kinds. A spread rule gives each chain its own start:

starts	Behaviour	When to use
from_prior	One draw per chain from each parameter’s resolved prior (fit.toml over model). A parameter with a flat prior falls back to bounds-uniform, with a warning.	The default whenever every estimated parameter declares a prior. Carries the scale a wide bound does not: at province scale a bounds-uniform draw is routinely a start the filter cannot score.
uniform_unconstrained	Stan-style: each chain draws z ~ Uniform(-2, 2) i.i.d. per parameter on the unconstrained scale, squashes it to u = σ(z), and maps that into bounds the transform-aware way lhs does (log-uniform for Log rates, linear otherwise). σ(±2) ≈ (0.119, 0.881) is a fixed interior band, so starts never sit on a bound and the radius is bounds-independent.	The default when some parameter has no prior. Boundary-avoiding (no -inf/zero-gradient starts) and scale-invariant; Stan’s well-tested default. Over-dispersed i.i.d. starts are the textbook basis for R-hat.
lhs	Latin-hypercube stratified sampling, scale-aware via the parameter’s transform: Log-typed rates are sampled in log space and exponentiated, so a single LHS pass spans orders of magnitude. Logit/None-typed parameters are sampled linearly in [lo, hi].	When you want stratified full-bounds coverage rather than i.i.d. draws — most useful for low-chain-count IF2 basin-finding.
uniform	Per-chain uniform random within natural-scale bounds. Chain 0 keeps the seeded start.	Legacy mode. Linear in natural space, so it clumps for Log-typed parameters at low chain count. Kept for reproducibility of pre-lhs results.
{ from_posterior = "<handle>" }	One row of a stored fit’s posterior per chain (or of a draws TSV named directly).	A warm start that keeps R̂ meaningful: the chains begin apart, where the earlier fit found mass.

A point rule puts every chain at one point, and R̂ is then reported as not assessed — chains that began together agree by construction:

starts	Behaviour	When to use
single	Every chain at the seeded [estimate].start (or its Transform-aware uniform fallback when start is omitted). Chains differ only by per-chain RNG.	See “When single is the right choice” below.
{ from_mle = "<handle>" }	Every chain at a stored fit’s point estimate.	Refining an optimizer’s estimate with more particles, where the chain-agreement  is about IF2’s own re-spreading of the swarm and the summary says so.
{ from_params = "<toml>" }	Every chain at the values in a flat params TOML.	A hand-picked start.

A handle is @label, a fit-id prefix, a method leaf’s run_id prefix, the leaf directory, or a fit.toml. A source whose stored verdict is not converged is refused; fit run --allow-nonconverged-source records the choice to start from it anyway.

Why uniform_unconstrained is the fallback. Stan initializes by drawing Uniform(-2, 2) on the unconstrained (transformed) scale and mapping back into the constrained space (mc-stan.org Reference Manual, “Initialization”). It’s boundary-avoiding — the sigmoid squash keeps every draw in the interior of [lo, hi], so no chain begins where the likelihood is -inf or the gradient degenerates — and scale-invariant: the same radius works whether a parameter is O(1) or O(1e6). camdl keeps the same transform-aware mapping lhs uses (Log rates in log space), so it retains the log-scale awareness that drove lhs past the legacy linear uniform, and simply replaces LHS’s stratified draw with i.i.d. over-dispersed draws — the standard recipe for diagnosing MCMC convergence via R-hat. The one thing it gives up versus lhs is guaranteed stratified coverage of each dimension’s range; a low-chain-count fit that needs that (e.g. the typhoid stratified scout, where 30 LHS chains beat 8 uniform-random chains by ~80,000 nats) can set starts = "lhs".

When single is the right choice. Three legitimate cases:

1. Single-chain runs (chains = 1) — there’s no per-chain spread to draw, so the three bounds-based rules collapse to the same draw.
2. Reproducibility-critical tests — single gives byte-identical chain starts across runs at the same seed; LHS/uniform draws shift if the RNG order changes upstream.
3. Deterministic NLopt with no spread desired — nl-sbplx, nl-bobyqa and nl-lbfgs are all deterministic, so single + chains > 1 gives N identical optimisations and the chain-agreement gate is uninformative; the run collapses to one chain and says so.

camdl profile dispatches the same way at each grid cell: --starts lhs --n-starts N draws N stratified per-cell starts across the non-focal estimated parameters; the focal parameters stay pinned to the grid point. --starts single reproduces the historical “every start at the same point, IF2 RNG provides the spread” behaviour.

Out-of-sample validation

The fit.toml holdout keys are applied at fit load (gh#585). Declaring holdout_after = <time> under [data] (a model-time number, a date under a calendar-anchored model, or last_obs - 6 weeks) truncates training to observations at t ≤ <time>; a [data.holdout] section binds explicit held-out files, which must lie strictly after each stream’s last training time (tail-only; an interior holdout is refused). The fit banner prints the training window, and the applied window is recorded in the fit’s fit.meta.json — the machine-checked proof camdl compare requires before labelling any score hold_out_tail.

[data]
holdout_after = "2020-03-01"

[data.observations]
weekly_cases = "data/cases.tsv"

camdl compare then scores such fits held-out by default: it runs the filter over the full series at the sealed θ̂, scores only the observations past the training boundary (each held-out week one-step-ahead, with earlier observations assimilated), verifies no leakage into θ̂, and stamps the trace hold_out_tail. Pass --in-sample to force the old full-series in-sample scoring. Comparing a holdout fit against a non-holdout fit is refused — those scores answer different questions.

The equivalent manual construction (also the fallback for a fit produced before the holdout was applied):

camdl data split data/cases.tsv --at-time 5474
camdl fit run fit_train.toml
camdl fit summary @train --params-only
camdl pfilter model.camdl --params theta.toml --data weekly_cases=data/cases.tsv --save-prequential preq --score-from 5474

camdl data split writes the train and holdout TSVs; point the fit’s [data.observations] at the training file, and write the fit’s θ̂ (from fit summary --params-only) to theta.toml. --score-from windows the trace to the held-out rows while still assimilating the training window.

Score the full series, not the holdout file alone. Binding only the held-out rows starts the filter from the model’s prior at the first held-out observation, with nothing assimilated from the training window: that scores p(y_holdout | θ̂) under unanchored state uncertainty, which penalises a stochastic model for exactly the uncertainty the training data would have resolved — and it usually trips the wide-first-window warning as well. Over the full series, each row of preq.tsv with t past the split time is a one-step-ahead score whose filter state assimilated the training window while θ̂ never saw the held-out observations. That subset is the held-out quantity.

Sum the log_score column over those rows with stream = "joint" for the held-out elpd. The TSV is tidy/long — a joint row per step followed by one row per scheduled stream — so summing every row double-counts each step.

Two numbers not to mis-read:

• the loglik pfilter prints covers the whole series, training window included. It is not a held-out score.
• the emitted trace is stamped conditioning = "in_sample" whatever window you score; the stamp is hardcoded in the trace builder today, and gh#585’s fix is what adds an honest one. The row subset is the honest number, not the trace’s own label.

Prediction quantiles

The pfilter trace includes both observation-space and state-space prediction quantiles:

• obs_mean, obs_q05, obs_q50, obs_q95 — full predictive distribution (process + observation noise). Data should fall inside the 5-95 ribbon ~90% of the time.
• state_mean, state_q05, state_q50, state_q95 — latent state quantiles mapped through the observation model mean. Process uncertainty only.

Both are on the observation scale (reported cases, not latent recoveries). The gap between the obs and state ribbons shows the observation model’s contribution to uncertainty.

Pfilter replicates

camdl pfilter model.camdl --params mle.toml --data d.tsv \
    --replicates 100 --output logliks.tsv

Runs N independent particle filters at different seeds. Reports loglik = -3804.9 ± 5.2 (100 replicates, N=5000).

See docs/camdl-inference-spec.md for the full specification.

Forecasting from the filtered state

camdl pfilter --save-final-state writes the particle ensemble at the last observation time — an unweighted sample from p(x_T | y_{1:T}) at the parameter vector the filter ran at. camdl simulate --init-state reads it back and runs forward from those states instead of the model’s init {} block, so a forecast does not re-run the filter:

# Filter the observed period, saving the state at the last observation.
camdl pfilter model.camdl --data train.tsv --params mle.toml \
    --particles 1000 --save-final-state final_particles.tsv

# Forecast eight weeks past the last observation, one replicate per particle.
camdl simulate model.camdl --params mle.toml \
    --init-state final_particles.tsv --replicates 1000 \
    --to "last_obs + 8 weeks" --fit fit.toml

The file is a TSV with a # camdl-final-state v1 t=<time> header, then one row per particle and one column per integer compartment — nothing else. The header’s t is the forecast origin and becomes the run’s t_start; --to sets the horizon (camdl simulate --help documents its anchored, absolute and calendar forms).

What a restart does and does not carry. Compartment counts are restored, by name. Interval accumulators start fresh at the origin: the forecast opens a new accumulation window, so the origin row carries zero flow, exactly as t_start does on an ordinary run.

That has one consequence worth expecting in the output. For an incidence stream, the observation at the forecast origin is a structural zero — its interval has zero length, so there is nothing yet to count. It is not a prediction that cases stop; read the forecast from the first row after the origin, or drop the origin row before plotting.

Everything that cannot be restored faithfully is refused rather than quietly defaulted:

Situation	Why it is refused
a real-valued compartment	the filter’s state is integer counts only, so the reservoir’s value at the origin is recorded nowhere
a reactive intervention	the policy’s observation history, once/cooldown gating, pending-effect queue and surveillance RNG stream are not in the file
--backend gillespie / ode	neither has a start-from-state seam, and the states came from the chain-binomial filter
--draws	the states belong to one θ; pairing them with unrelated posterior draws forecasts a state and a parameter vector that never went together
--replicates ≠ the row count	replicate i restores row i, and a prefix of a post-resampling swarm is not an exchangeable subsample

The origin must coincide with an output-emit time — a time between snapshots is an error, never snapped to a neighbour.

An effect scheduled at the origin time is not applied again: the saved state is post-effect, because the run that saved it applied every effect due at that time before recording it. Effects after the origin fire on the forecast’s own grid, origin + k·dt (language spec §13.2).

Parameter uncertainty is not propagated by this form: every replicate runs at the filter’s θ, so the spread is process noise plus filtering uncertainty in x_T, not posterior uncertainty in θ.

What that costs: the intervals are too narrow. Conditioning on a single θ̂ discards a variance component the forecast genuinely has, so a 90% band from this workflow will cover the truth less than 90% of the time, and increasingly so the further out you forecast (parameter error compounds over the horizon while process noise averages). Treat these bands as a lower bound on uncertainty, and do not report them as calibrated predictive intervals. Pairing each latent state with its own posterior draw is what fixes it — that is --init-state fit, below.

Forecasting from the paired (θ, X) posterior

--init-state <file> conditions on where the epidemic is but runs at one θ; --draws posterior propagates parameter uncertainty but starts every draw from init {} at t = 0, which is a replay rather than a forecast. --init-state fit is both at once: draw i restores the terminal state of its own saved latent path and runs forward under its own θ_i.

camdl simulate model.camdl --fit results/fits/<stem>-<hash>/ \
    --draws posterior --init-state fit \
    --to "last_obs + 8 weeks" --scenario control_25

The origin is the last snapshot of each posterior draw’s saved latent path, which a PGAS run writes to <method>/chain_*/trajectories.tsv. At the terminal observation time the smoothing distribution equals the filtering distribution — no future data remains to condition on — so that row is a draw from p(x_T | y_{1:T}) carrying its own parameters. Interior origins are deliberately not offered: iterating forward from a smoothing distribution is not one of the three objects in Särkkä’s taxonomy (Särkkä 2013, Bayesian Filtering and Smoothing, CUP, doi:10.1017/CBO9781139344203, §1.3). This is the published idiom for the workflow — Funk, Camacho, Kucharski, Lowe, Eggo & Edmunds (2019), PLoS Comput Biol 15(2):e1006785, doi:10.1371/journal.pcbi.1006785, use “the final values of estimated state trajectories as initial values for the forecasts”.

Only draws with a saved path can be forecast, and the run says how many. n_trajectories and thin are independent knobs, so a run that retains 606 posterior draws may have saved only 300 latent paths. The run reports the subset it used against the full posterior:

simulate: --init-state fit → 300/606 posterior draws have a saved latent path;
forecasting that paired (θ_i, X_i) subset from t = 196 (the other 306 draws have
no state to fork)
  method 'pgas' (results/fits/…/draws.tsv)

Raise n_trajectories on the fit’s [method] if you want the whole posterior forecastable. A PGAS fit that saved no paths — and every PMMH or particle-filter fit, which never writes any — is refused by name rather than falling back to init {}, which would look like a forecast and be a free-forward replay. An ODE fit is refused separately and for a different reason: it stores no paths because X is recomputed from θ, and the post-origin re-integration seam is not wired (gh#325).

The refusals for a reactive policy, a non-chain-binomial --backend and an off-grid origin apply here exactly as they do to the file source. --scenario overlays and anchored --to horizons compose with both.