Learning in the Loop

Where machine learning belongs in a probabilistic robot — learned models inside an unchanged Bayes filter, calibration as the admission test, and the filter itself as a differentiable program you can train without ground truth.

The problem: an accurate model that lies

An engineer wires up Rusty's new LiDAR on a Tuesday afternoon. The datasheet says ranging accuracy $\pm 5$ cm, so the beam model of Chapter 10 gets $\sigma_{hit} = 0.05$ m, a small spike at $z_{max}$ for dropouts, and a thin uniform floor for whatever else happens. No unexpected-obstacle component: the room is empty, the model has no reason to expect anything in front of the wall.

The model is accurate. Its peak sits on the correct range. Ask it for the most likely reading and it is right more often than the model the robot shipped with. Then Monte Carlo localization is switched on and the robot loses itself inside twenty seconds — not by drifting, but by committing, hard, to a pose that is 30 cm wrong and refusing every subsequent scan that disagrees.

Watch one full sweep of the slider before reading on. Three things are happening, and each one has a piece of mathematics later in this chapter.

The reliability curve sags below the diagonal. The diagonal is the promise: of all the observations for which the model claimed "at most this far out", that fraction really should land at most that far out. A curve below the diagonal means every interval the model quotes catches fewer observations than advertised. At $\sigma_{hit} = 5$ cm the expected calibration error is $0.107$ — the model's 50% interval is closer to a 30% interval — and the $\chi^2$ statistic of its PIT histogram is over 1400 on 14 degrees of freedom. That is not a borderline call; that is a model that is lying in a way you could detect from a hundred readings.

The PIT histogram is a U. Readings keep landing in the tails the model declared empty. The shape of the miscalibration tells you its cause: a U means too narrow, a hump means too wide, a tilt means biased. This is a diagnostic instrument, not a score, and it is the single most useful plot in this chapter.

The filter downstream cares, and it does not care about what you expected. The two right-hand tiles come from a real importance-weighting step over twelve particle clouds in the Apartment. As the claimed width shrinks, the effective sample size collapses — at $\sigma_{hit} = 3$ cm a 150-particle cloud is doing the work of one particle — and the posterior error rises even as the model looks sharper and more confident.

Building intuition

Learning, 1999 style — and what it cannot fix

Learning is not something modern practice grafted onto Thrun, Burgard and Fox. The 1999–2000 draft already fits the beam model's intrinsic parameters by expectation-maximization (learn_intrinsic_parameters, Table 6.2, derived in Chapter 10); it already trains a neural network to invert the measurement model by sampling triplets from the forward model (§9.3.2); it learns maps with EM, and it learns POMDP value functions by nearest-neighbor function approximation. The epigraph to this chapter is that book noticing its own frontier.

So run it. Take 4000 range readings drawn from the true sensor in the Apartment and hand them to learn_intrinsic_parameters, initialized at the engineer's hand-tuned guess. EM does its job: the mean log-likelihood climbs from $-0.923$ to $-0.224$ nats per beam in forty iterations and $\sigma_{hit}$ moves from $0.050$ m to $0.112$ m, essentially the truth. And yet the fitted model is still measurably miscalibrated, with an ECE of $0.029$ against the true model's $0.021$ — because EM inherited a defect from its initialization that no amount of iterating can repair.

The engineer's model set $z_{short} = 0$: no unexpected-obstacle component at all. In EM, a component's responsibility is proportional to its weight, so a weight of exactly zero produces zero responsibility, which produces a new weight of exactly zero. A component initialized at zero can never come back. EM did the only thing left to it and inflated $z_{rand}$ from $0.05$ to $0.142$, smearing a uniform floor across the whole range to absorb readings that came from an exponential concentrated in front of the wall.

That is the first honest lesson about learning in a robot. Maximum likelihood inside a family only searches that family. When the family is wrong — and it always is, a little — the fit will spend its freedom compensating in whatever direction the parameterization allows, and the result can be an excellent likelihood attached to a badly shaped uncertainty.

The question is not "how accurate", it is "how honest"

Point accuracy asks: is the model's best guess close to the truth? Calibration asks something the filter cares about far more: when the model says 70%, does it happen 70% of the time?

Gneiting, Balabdaoui and Raftery gave the field the sentence that settles the priority. The goal of probabilistic forecasting, they argue, is to maximize the sharpness of the predictive distributions subject to calibration — sharpness being how concentrated the forecast is, a property of the forecast alone, and calibration being the statistical consistency between forecasts and what actually happens, a joint property of the forecast and the world. Sharpness is what you want. Calibration is what you are allowed. Take them in the wrong order and you get exactly the widget above: a model that wins on confidence and loses on reality.

Overconfidence is a filter-killer, not an aesthetic complaint

Here is the mechanism, before the algebra. A particle filter weights each particle by $w^{[i]} \propto p(z_t \mid x_t^{[i]}, m)$. Sharpen the likelihood and you sharpen the weights: one particle takes almost all the mass, the effective sample size collapses, and resampling clones that single particle four hundred times. The cloud now has one hypothesis. If the truth is not underneath it, there is nothing left to re-weight and no recovery is possible — the particle deprivation of Chapter 12, arriving not from bad luck but from an overconfident model.

The knob in that widget, $\kappa$, is the honest way to talk about this. It does not change what the model says, only how loudly: $p(z \mid x, m)^\kappa$. And it is not an artificial device — Chapter 10 showed that treating $K$ correlated beams as independent multiplies every log-odds by roughly $K / K_{\text{eff}}$, which is exactly a tempering exponent greater than one that nobody chose and nobody wrote down.

The mathematics

Notation for this chapter

Proper scoring rules

You cannot train a model without a loss, and not every loss is safe. A scoring rule that a forecaster can improve by misreporting its beliefs is a machine for producing liars.

Strict propriety is what makes the negative log-likelihood the default training loss for anything that will live inside a Bayes filter: it is the loss whose unique optimum is the truth, so a model that minimizes it has no incentive left to shade its uncertainty in either direction.

Calibration, PIT, and ECE

The PIT is the whole diagnostic in one line, and it is worth seeing why it works. If $z$ really is drawn from the density the model claims, then $F(z)$ is uniform — the classic inverse-transform argument run backwards. So any structure in the histogram of PIT values is evidence of misspecification, and its shape names the defect:

One technical repair matters for range sensors. The beam model has an atom: a lump of probability mass $z_{max}$ sitting exactly on the maximum range. For a discontinuous CDF the plain PIT is not uniform even for a perfect model — every dropout beam would pile up at $v = 1$. The standard fix is the randomized PIT, which spreads each atom uniformly across the jump it makes:

For a continuous forecast $F(z^-) = F(z)$ and this reduces to the plain PIT, so it costs nothing to always use it. Skip it and you will spend an afternoon "fixing" the calibration of a model whose only sin was reporting dropouts honestly.

The two numbers this returns are not redundant, and it is important to hold both: the mean NLL measures accuracy and sharpness together, and the ECE measures honesty alone. A model can win one and lose the other. The benchmark at the end of this chapter is a table in which exactly that happens.

Calibration and sharpness are different axes

The tempered update, and what it costs

The repair is one scalar, fitted on data the model has never seen.

Three properties make this the first thing to try and the last thing to skip. It has one parameter, so it cannot overfit a validation set of any reasonable size. It is fitted on held-out data, so it corrects the model rather than memorizing the training set. And for a density it does not reorder hypotheses at all — raising a Gaussian likelihood to the power $\kappa$ is identical to scaling its width by $s = 1/\sqrt{\kappa}$, and neither operation moves the mode. It buys calibration without costing accuracy. That is why Guo et al.'s temperature scaling took over the classification literature, and it transposes to a range density without modification.

Learning the model from the model you already have

Thrun et al.'s §9.3.2 recipe is still the right one, and it is worth stating in its modern form because it explains where a learned sensor model's training data comes from when the robot has no labels: sample the forward model, fit the inverse. You have a simulator, a map, and a physical model of the sensor; that is a generator of $(x, z)$ pairs at any volume you like.

The filter is a differentiable program

The reason to care is not that gradient descent is fashionable. It is that hand-tuning $R_t$ and $Q_t$ is the single most common and least principled activity in applied estimation, and there is a loss function that does it properly — one that does not require ground truth.

A worked example you can check by hand

Take the shortest possible trajectory: one step, $a = 1$, $b = 0$, a known initial state ($\mu_0 = 0$, $\sigma_0 = 0$), and unit noise parameters $r = q = 1$, so $\theta = (\log r, \log q) = (0, 0)$. The robot reads $z_1 = 2$.

Run the filter:

The evidence loss is one Gaussian term:

Now the gradients. Since $\sigma_0 = 0$, the prediction $\bar\mu_1$ does not depend on $\theta$ at all, so $\partial \nu_1 = 0$ and only the first bracket of F5 Step 5 survives. With $\partial S_1/\partial \log r = r = 1$ and $\partial S_1/\partial \log q = q = 1$:

Both gradients are negative, so descent grows both variances: the filter was surprised — a $\nu = 2$ innovation against a claimed $S = 2$ is a $1.41\sigma$ event whose square is twice the claimed variance — and the honest response is to admit more uncertainty.

And the two gradients are equal, which is the lesson hiding in the arithmetic. After one step from a known initial state, the loss depends on $\theta$ only through the sum $S_1 = r + q$: no single measurement can tell process noise from measurement noise. What separates them is time. Process noise accumulates across steps and measurement noise does not, so the autocorrelation of the innovation sequence carries the information that identifies the split — which is why the trainer below needs a few hundred steps of trajectory, not a few.

Two things in that widget deserve more than a glance. The trainer starts with a filter that believes the cart is wildly unpredictable and the sensor is nearly perfect — the classic tuning mistake, in both directions at once — and the purple band is far too tight for it. As the loss falls, the band breathes out to a width that contains the truth about 95% of the time, and the mean NEES walks to

1. Nobody told it what 1 was. The evidence loss never saw the gray dashed line.

And the "cheat" toggle is there to make the alternative concrete. With ground truth you can minimize the state NLL directly, which is what a differentiable filter trained in simulation usually does. It converges faster and to a slightly different place, and it is unavailable on any robot that is not in a motion-capture room. When someone reports a learned filter that was trained on state supervision and deployed without it, this toggle is the question to ask them.

Where the gradient dies

Everything above differentiates cleanly because the Kalman filter is a chain of smooth functions. Particle filters are not: they contain a categorical draw, and a categorical draw has no useful derivative.

Learning to act: distributions over actions

The chapter's thesis has a control-side sibling. Chapter 23 already samples action sequences and weights them; a diffusion policy learns the sampler.

The argument for it is the argument for particle filters, transposed. Suppose an obstacle can be passed on the left or on the right, and the demonstrations show both. A policy trained to regress the mean action drives straight into the obstacle — the average of two good plans is not a plan. Learning a sampler preserves the two modes, and the robot commits to one of them.

That is genuinely in the spirit of this book: represent what you cannot summarize. But note the asymmetry, because it is the reason this section is short and has no widget. A diffusion policy is a sampler, not a density: you can draw actions from it, but you cannot cheaply evaluate the probability of an action, so it does not compose with Bayes' rule the way a calibrated sensor model does. Everything else in this chapter can be multiplied into a belief. This cannot. The full derivation of the diffusion ELBO is out of scope here, and the honest reference is Chi et al., whose IJRR version is the standard treatment.

Implementation in Rust

Three modules carry the chapter: the instrument, the learned model, and the trainer.

The instrument

Calibration is model-agnostic on purpose. Grading consumes PIT values and log scores, never a model, so the same code grades a hand-tuned mixture, a network, and a Kalman filter's innovation sequence.

use rand::Rng;

/// Anything that can serve the Bayes filter as p(z | x, m).
///
/// Chapter 10's `BeamModel` and `LikelihoodField` implement it, and so does
/// this chapter's `LearnedBeamModel`. The second and third methods are what
/// make a model *auditable*: a density you cannot integrate is a density you
/// cannot grade, and an atom you do not declare will masquerade as a defect.
pub trait ObservationModel {
    /// log p(z | x, m) for a whole scan — what the filter multiplies in.
    fn log_likelihood(&self, z: &Scan, x: &Se2, m: &OccGrid) -> f64;

    /// Per-beam predictive CDF F(z_k | z*), needed for the PIT.
    fn beam_cdf(&self, z_k: f64, z_star: f64) -> f64;

    /// Mass of the atom at `z_k`, if the model has one. The beam mixture's
    /// max-range spike is an atom; a likelihood field has none.
    fn beam_atom(&self, _z_k: f64) -> f64 {
        0.0
    }
}

pub const N_BINS: usize = 15;

#[derive(Debug, Clone, Copy)]
pub struct ReliabilityBin {
    pub nominal: f64,
    pub empirical: f64,
    pub count: usize,
}

#[derive(Debug, Clone)]
pub struct CalibrationReport {
    pub bins: [ReliabilityBin; N_BINS],
    /// D25.3, in the quantile-binned form: mean |observed − claimed|.
    pub ece: f64,
    /// Largest gap from the diagonal — the Kolmogorov statistic.
    pub max_calibration_error: f64,
    /// Mean negative log score, nats per beam. Accuracy *and* sharpness.
    pub mean_nll: f64,
    /// Pearson χ² of the PIT histogram against uniform, N_BINS − 1 d.o.f.
    pub pit_chi2: f64,
}

/// Randomized PIT: spread each atom uniformly across the jump it makes, so a
/// dropout beam does not pile up at v = 1 and look like miscalibration.
pub fn randomized_pit<M: ObservationModel>(
    model: &M,
    z_k: f64,
    z_star: f64,
    rng: &mut impl Rng,
) -> f64 {
    let upper = model.beam_cdf(z_k, z_star);
    let lower = upper - model.beam_atom(z_k);
    (lower + rng.random::<f64>() * (upper - lower)).clamp(0.0, 1.0)
}

pub fn calibrate<M: ObservationModel>(
    model: &M,
    eval: &[(f64, f64)], // (z*, z) pairs from a held-out set
    log_density: impl Fn(f64, f64) -> f64,
    rng: &mut impl Rng,
) -> CalibrationReport {
    let mut pit: Vec<f64> = eval
        .iter()
        .map(|&(z_star, z)| randomized_pit(model, z, z_star, rng))
        .collect();
    let mean_nll = -eval
        .iter()
        .map(|&(z_star, z)| log_density(z, z_star))
        .sum::<f64>()
        / eval.len() as f64;

    // Sorting turns every empirical-CDF query into an O(1) index lookup.
    pit.sort_by(|a, b| a.partial_cmp(b).unwrap());
    let n = pit.len() as f64;

    let mut bins = [ReliabilityBin { nominal: 0.0, empirical: 0.0, count: 0 }; N_BINS];
    let (mut ece, mut worst) = (0.0, 0.0_f64);
    for b in 0..N_BINS {
        let nominal = (b + 1) as f64 / N_BINS as f64;
        let count = pit.partition_point(|&v| v <= nominal);
        let empirical = count as f64 / n;
        let gap = (empirical - nominal).abs();
        ece += gap / N_BINS as f64;
        worst = worst.max(gap);
        bins[b] = ReliabilityBin { nominal, empirical, count };
    }

    // χ² of the histogram (not the CDF) against a uniform expectation.
    let expected = n / N_BINS as f64;
    let mut counts = [0usize; N_BINS];
    for &v in &pit {
        counts[((v * N_BINS as f64) as usize).min(N_BINS - 1)] += 1;
    }
    let pit_chi2 = counts
        .iter()
        .map(|&c| (c as f64 - expected).powi(2) / expected)
        .sum();

    CalibrationReport { bins, ece, max_calibration_error: worst, mean_nll, pit_chi2 }
}

The learned model

The hand-tuned beam model claims one $\sigma_{hit}$ for every range, which is false on every real LiDAR: a beam returning from 7 m is not the beam returning from 70 cm. The learned model makes the mixture parameters functions of the expected range — a two-layer MLP that emits four mixture logits and a log-width, trained by the NLL that F1 licensed.

The design decision worth recording: training happens natively with candle; the browser never trains a net. Inference is a dependency-free forward pass — two matrix multiplies and a GELU — so the same weights run in the WASM widget with no framework attached.

use candle_core::{DType, Device, Result, Tensor};
use candle_nn::{loss, AdamW, Linear, Module, Optimizer, ParamsAdamW, VarBuilder, VarMap};

/// Learned per-beam range density: MLP(z*) -> (mixture logits, log σ_hit).
///
/// The *shape* of the mixture is kept from Chapter 10 — hit, short, max, rand —
/// because that structure is physics and the network has nothing better to
/// offer. What is learned is how the four weights and the hit width vary with
/// range, which is exactly the part the hand-tuned model got wrong.
pub struct BeamNet {
    l1: Linear,
    l2: Linear,
    max_range: f64,
}

impl BeamNet {
    pub fn new(vb: VarBuilder, hidden: usize) -> Result<Self> {
        Ok(Self {
            l1: candle_nn::linear(2, hidden, vb.pp("l1"))?,
            l2: candle_nn::linear(hidden, 5, vb.pp("l2"))?,
            max_range: 8.0,
        })
    }

    /// Features are deliberately dimensionless: a network trained on metres and
    /// deployed on a sensor with a different range would otherwise extrapolate.
    fn features(&self, z_star: &Tensor) -> Result<Tensor> {
        let u = (z_star / self.max_range)?;
        Tensor::cat(&[&u, &u.sqrt()?], 1)
    }

    /// Negative log-likelihood of a batch under the predicted mixture. This is
    /// the strictly proper log score of D25.1 — nothing else is trained.
    pub fn nll(&self, z_star: &Tensor, z: &Tensor) -> Result<Tensor> {
        let out = self.l2.forward(&self.l1.forward(&self.features(z_star)?)?.gelu()?)?;
        let logits = out.narrow(1, 0, 4)?;
        let log_w = candle_nn::ops::log_softmax(&logits, 1)?; // weights on the simplex
        let sigma = (out.narrow(1, 4, 1)?.exp()? + 1e-3)?; // positivity by construction

        // log of each component density, then log-sum-exp against the weights.
        let comps = Tensor::cat(
            &[
                self.log_hit(z, z_star, &sigma)?,
                self.log_short(z, z_star)?,
                self.log_max(z)?,
                self.log_rand(z)?,
            ],
            1,
        )?;
        let log_mix = (comps + log_w)?.log_sum_exp(1)?;
        log_mix.neg()?.mean_all()
    }
}

/// One epoch of `learn_observation_model`, lines 8–9.
pub fn train(pairs: &[(f64, f64)], epochs: usize, dev: &Device) -> Result<VarMap> {
    let varmap = VarMap::new();
    let vb = VarBuilder::from_varmap(&varmap, DType::F64, dev);
    let net = BeamNet::new(vb, 64)?;
    let mut opt = AdamW::new(varmap.all_vars(), ParamsAdamW { lr: 1e-3, ..Default::default() })?;

    let z_star = Tensor::from_iter(pairs.iter().map(|p| p.0), dev)?.reshape((pairs.len(), 1))?;
    let z = Tensor::from_iter(pairs.iter().map(|p| p.1), dev)?.reshape((pairs.len(), 1))?;

    for _ in 0..epochs {
        let loss = net.nll(&z_star, &z)?;
        opt.backward_step(&loss)?;
    }
    Ok(varmap)
}

The trainer

DiffKf1d is the whole of F5, transcribed. It is about forty lines, it has no dependencies, and it is the engine the trainer widget runs in the browser — the derivation is the source code.

/// A scalar Kalman filter with θ = (log r, log q) exposed and differentiated.
///
/// Logs, not variances, so that descent cannot propose a negative variance:
/// the constraint lives in the parameterization rather than in a clamp.
pub struct DiffKf1d {
    pub log_r: f64,
    pub log_q: f64,
}

pub enum DkfLoss {
    /// F4. Scored on one-step-ahead predictions of measurements the robot took.
    Evidence,
    /// The supervised alternative: available in simulation and motion capture,
    /// and nowhere else.
    StateNll,
}

impl DiffKf1d {
    /// Loss and its exact gradient in a single forward pass (F5).
    pub fn loss_and_grad(&self, log: &TrajLog1d, loss: DkfLoss) -> (f64, f64, f64) {
        let (r, q) = (self.log_r.exp(), self.log_q.exp());
        let (a, b) = (log.a, log.b);

        let (mut mu, mut sigma2) = (log.mu0, log.sigma0 * log.sigma0);
        // Sensitivities ∂(μ, σ²)/∂θ. The initial belief is given, not learned.
        let (mut d_mu, mut d_sig) = ([0.0_f64; 2], [0.0_f64; 2]);
        let (mut total, mut grad) = (0.0, [0.0_f64; 2]);

        for t in 0..log.z.len() {
            let u = log.u.get(t).copied().unwrap_or(0.0);

            let mu_bar = a * mu + b * u;
            let sigma2_bar = a * a * sigma2 + r;
            let d_mu_bar = [a * d_mu[0], a * d_mu[1]];
            // ∂r/∂log r = r, ∂r/∂log q = 0.
            let d_sig_bar = [a * a * d_sig[0] + r, a * a * d_sig[1]];

            let nu = log.z[t] - mu_bar;
            let s = sigma2_bar + q;
            let d_nu = [-d_mu_bar[0], -d_mu_bar[1]];
            let d_s = [d_sig_bar[0], d_sig_bar[1] + q];

            let k = sigma2_bar / s;
            let d_k = [
                (d_sig_bar[0] - k * d_s[0]) / s,
                (d_sig_bar[1] - k * d_s[1]) / s,
            ];

            let mu_new = mu_bar + k * nu;
            let sigma2_new = (1.0 - k) * sigma2_bar;

            match loss {
                DkfLoss::Evidence => {
                    total += 0.5 * ((core::f64::consts::TAU * s).ln() + nu * nu / s);
                    for p in 0..2 {
                        // Negative when ν² > S: descent widens a surprised filter.
                        grad[p] += 0.5 * ((1.0 - nu * nu / s) * (d_s[p] / s) + 2.0 * nu * d_nu[p] / s);
                    }
                }
                DkfLoss::StateNll => {
                    let e = log.x[t] - mu_new;
                    total += 0.5 * ((core::f64::consts::TAU * sigma2_new).ln() + e * e / sigma2_new);
                    for p in 0..2 {
                        let d_mu_new = d_mu_bar[p] + d_k[p] * nu + k * d_nu[p];
                        let d_sig_new = (1.0 - k) * d_sig_bar[p] - d_k[p] * sigma2_bar;
                        grad[p] += 0.5
                            * ((1.0 - e * e / sigma2_new) * (d_sig_new / sigma2_new)
                                - 2.0 * e * d_mu_new / sigma2_new);
                    }
                }
            }

            d_mu = [
                d_mu_bar[0] + d_k[0] * nu + k * d_nu[0],
                d_mu_bar[1] + d_k[1] * nu + k * d_nu[1],
            ];
            d_sig = [
                (1.0 - k) * d_sig_bar[0] - d_k[0] * sigma2_bar,
                (1.0 - k) * d_sig_bar[1] - d_k[1] * sigma2_bar,
            ];
            mu = mu_new;
            sigma2 = sigma2_new;
        }

        (total, grad[0], grad[1])
    }
}

Soft resampling is the same low-variance comb as Thrun's Table 4.4, walked through $q$ instead of $w$, with the importance ratio attached to each survivor:

use rand::Rng;

/// `soft_resample(X_t, λ, rng)` — F6. Additive: `ParticleFilter` opts in by
/// calling this instead of `low_variance_sampler`, and nothing else changes.
pub fn soft_resample<S: Clone>(
    particles: &mut Vec<Particle<S>>,
    lambda: f64,
    rng: &mut impl Rng,
) -> f64 {
    let m = particles.len();
    let total: f64 = particles.iter().map(|p| p.weight).sum();
    let w: Vec<f64> = particles.iter().map(|p| p.weight / total).collect();
    let lambda = lambda.clamp(0.0, 1.0);
    let q: Vec<f64> = w.iter().map(|wi| lambda * wi + (1.0 - lambda) / m as f64).collect();

    let step = 1.0 / m as f64;
    let start = rng.random_range(0.0..step);
    let (mut out, mut acc, mut i, mut transmission) = (Vec::with_capacity(m), q[0], 0usize, 0.0);

    for k in 0..m {
        let u = start + k as f64 * step;
        while u > acc && i < m - 1 {
            i += 1;
            acc += q[i];
        }
        out.push(Particle { state: particles[i].state.clone(), weight: w[i] / q[i] });
        // ∂ log w′ᵢ / ∂ log wᵢ — the fraction of the gradient this survivor passes.
        transmission += ((1.0 - lambda) / m as f64) / q[i];
    }

    *particles = out;
    transmission / m as f64 // ≈ 1 − λ, exactly, in expectation
}

The tests that pin it

#[test]
fn worked_example_ch25_one_step_evidence() {
    // One step, r = q = 1, σ₀ = 0, z₁ = 2 — the example derived in the text.
    let log = TrajLog1d { a: 1.0, b: 0.0, u: vec![0.0], z: vec![2.0], x: vec![0.0],
                          mu0: 0.0, sigma0: 0.0 };
    let kf = DiffKf1d { log_r: 0.0, log_q: 0.0 };

    let (loss, d_log_r, d_log_q) = kf.loss_and_grad(&log, DkfLoss::Evidence);

    // L = ½(ln 4π + 2) = 2.265512…
    assert_relative_eq!(loss, 0.5 * ((4.0 * PI).ln() + 2.0), epsilon = 1e-12);
    // Equal and negative: the filter was surprised, and one measurement cannot
    // separate process noise from measurement noise.
    assert_relative_eq!(d_log_r, -0.25, epsilon = 1e-12);
    assert_relative_eq!(d_log_q, -0.25, epsilon = 1e-12);
}

#[test]
fn analytic_gradient_matches_central_differences() {
    let log = simulate_traj_1d(120, 0.96, 0.1, 0.05, 0.4, &mut SmallRng::seed_from_u64(2025));
    for loss in [DkfLoss::Evidence, DkfLoss::StateNll] {
        for theta in [(-1.5, 0.4), (0.3, -0.8), (-3.0, -3.0)] {
            let kf = DiffKf1d { log_r: theta.0, log_q: theta.1 };
            let (_, g_r, g_q) = kf.loss_and_grad(&log, loss);
            let (fd_r, fd_q) = finite_diff_grad(&kf, &log, loss, 1e-5);
            assert_relative_eq!(g_r, fd_r, max_relative = 1e-5);
            assert_relative_eq!(g_q, fd_q, max_relative = 1e-5);
        }
    }
}

#[test]
fn evidence_training_recovers_the_noises_without_ground_truth() {
    // 1500 steps generated with r = 0.08, q = 0.25; start deliberately wrong
    // in both directions and never look at `log.x`.
    let log = simulate_traj_1d(1500, 1.0, 0.05, 0.08, 0.25, &mut SmallRng::seed_from_u64(1234));
    let mut kf = DiffKf1d { log_r: 0.4_f64.ln(), log_q: 0.02_f64.ln() };
    for _ in 0..12_000 {
        kf.sgd_step(&log, DkfLoss::Evidence, 0.05 / log.z.len() as f64);
    }
    assert!((kf.log_r.exp() / 0.08).ln().abs() < 0.2); // r → 0.085
    assert!((kf.log_q.exp() / 0.25).ln().abs() < 0.1); // q → 0.250
}

Putting it together: the clinic benchmark

The chapter's empirical thesis needs a table, and the table is generated by cargo run --example clinic — the same computation the Calibration Clinic runs live in your browser, so every number here is one you can reproduce by clicking a button in the widget above.

The setup: the true sensor is a beam mixture with $\sigma_{hit} = 0.11$ m and a real unexpected-obstacle component ($z_{short} = 0.08$). The model family available to the engineer is mis-specified — it has no $z_{short}$ at all. Held-out evaluation is 1200 beams from random poses in the Apartment; downstream is one importance-weighting step over twelve 150-particle clouds, each displaced about 42 cm from the truth.

Read the last three rows against each other, because that is the entire chapter in six numbers.

The best likelihood is not the best filter. EM wins the NLL column decisively — it is the only model here that fits the mixture weights as well as the width — and it still loses the posterior-error column to a model with a worse NLL by 0.09 nats. The widened, best-calibrated model gives up score in exchange for honesty, and honesty is what the filter downstream converts into accuracy: a 60% larger effective sample size and 5% less error, from a model that any leaderboard would rank below the one it beats.

Overconfidence is expensive and calibration is cheap. The hand-tuned model is not subtly worse. It is 4.5× the ECE of anything fitted, and it leaves the posterior 0.292 m from the truth when the prior was already only 0.424 m away — it converts a scan into a third of the correction it should have bought. Fixing it required no retraining and no new data: one scalar, fitted on 800 held-out beams, in about fifty function evaluations.

And the honest model still wins on score. The true sensor's NLL of 0.223 is far below everything else, and no amount of temperature scaling closes that gap, because tempering cannot invent the $z_{short}$ component the family is missing. Calibration is a repair, not a substitute for getting the structure right. That is worth saying loudly in a chapter about learning: the largest single improvement available in this table comes from the model family, not from the fitting procedure.

The honesty ledger

The chapter closes where it should — with a line drawn between what is principled and what is alchemy, as of 2026.

Principled. A learned component that emits a calibrated density and is dropped into an unchanged Bayes filter: this is just a better $p(z_t \mid x_t, m)$, and every theorem in Parts II through V still holds. Training noise parameters by maximizing evidence, which is prediction-error identification with a fifty-year pedigree and needs no ground truth. Post-hoc calibration on held-out data, which has one parameter and cannot lie to you in any interesting way. Differentiable filters that keep the filter's structure as an algorithmic prior and learn only its contents.

Still alchemy. End-to-end "filters" whose belief is a vector of activations that nobody has graded — if you cannot draw its reliability diagram, do not call it a posterior. Learned resampling heuristics published without a variance analysis. Sim-to-real transfer of learned sensor models: our own learned model was trained in the simulator that also evaluates it, which is exactly the circularity flagged after learn_observation_model, and no reliability diagram computed inside the simulator says anything about the building. And diffusion policies' uncertainty: they are samplers, not densities, so they do not compose with Bayes' rule, and a policy that produces diverse actions is not thereby a policy that knows what it does not know.

A useful sanity check before adopting any of this: Greenberg, Yannay and Mannor (NeurIPS 2023) showed that a Kalman filter optimized as carefully as the network it is compared against closes much of the reported gap in the differentiable-filtering literature. The technique in this chapter — train the classical filter's parameters properly — is also the baseline that most learned filters should have been compared against.

Chapter 26 puts this to work: its stack exposes a calibrated learned sensor model as a toggle, and its retrospective quotes this table.

Exercises

1. Foundation exerciseFDifficulty 2 of 3Propriety, and a score you can game

   Prove D25.1 for the log score without looking at derivation F1 — expected score equals entropy plus KL, and Gibbs' inequality does the rest. Then construct a scoring rule that a beam model could game: find $S$ and two densities $p \ne q$ with $\E_{z\sim p}[S(q,z)] < \E_{z \sim p}[S(p,z)]$, and say in one sentence what a model trained on it would do to $\sigma_{hit}$. (Hint: reward the density at the observation's mode rather than at the observation.)

2. Foundation exerciseFDifficulty 3 of 3Differentiate the filter yourself

   Derive $\partial \mathcal{L} / \partial \log q$ for the scalar Kalman filter under the evidence loss, carrying the sensitivity recursion of F5 over two steps rather than one. Verify your expression against DiffKf1d::loss_and_grad by central differences at $\theta = (-1.5,\, 0.4)$; the test skeleton is analytic_gradient_matches_central_differences above. Then explain why the one-step worked example has $\partial\mathcal{L}/\partial\log r = \partial\mathcal{L}/\partial\log q$ and the two-step version does not.

3. Foundation exerciseFDifficulty 3 of 3The ESS lemma

   Derive $\mathrm{ESS}/M \approx e^{-\kappa^2 s^2}$ from the Gaussian moment generating function (F3, Steps 3–4), stating precisely where you used the assumption that the per-particle log-likelihoods are Gaussian. Then compute how large a spread $s$ your particle filter can tolerate before $\kappa = 1$ leaves fewer than 5% of $M = 1000$ particles effective, and compare with the spread the Tempering Lab reports for a 16-beam likelihood-field scan.

4. Conceptual exerciseCDifficulty 2 of 3Predict, then temper

   Before touching the Tempering Lab: predict what happens to the kidnap-recovery percentage at $\kappa = 3$, and predict whether $\kappa = 0.5$ helps or hurts tracking accuracy between kidnaps. Run both for at least three kidnap cycles each. Explain both outcomes using F3, and say which of the two effects — precision growth or ESS collapse — dominates the failure you saw.

5. Conceptual exerciseCDifficulty 2 of 3Find the two optima

   In the Calibration Clinic, use the slider to locate the width that minimizes the mean NLL and the width that minimizes the ECE. They are different. Before reading the two right-hand tiles, predict which of those two models gives the lower posterior error — then check. Now explain the gap using the reliability diagram: which end of the PIT histogram is the NLL optimum trading away, and what is it buying with it?

6. Practical exercisePDifficulty 2 of 3Grade the likelihood field

   Implement ObservationModel and calibrate for Chapter 10's LikelihoodField. Its per-beam density is a Gaussian on the distance to the nearest occupied cell plus a uniform floor, so its CDF is available in closed form. Is it calibrated? Produce its PIT histogram over 1200 beams in the Apartment and identify where it bulges — then explain the bulge from the model's structure (hint: what does the likelihood field do with a max-range beam, and with a beam that hits an unmapped obstacle?).

7. Practical exercisePDifficulty 3 of 3Train through the resampler

   Wire soft_resample into Chapter 8's ParticleFilter and train a single parameter — the sensor model's width — end-to-end through three steps of 1-D hallway MCL, using the state NLL of the filter mean as the loss. Sweep $\lambda \in \{0, 0.25, 0.5, 0.75, 1\}$ and plot both the gradient magnitude reaching the parameter and the variance of the resulting estimate. Does the learned width agree with the post-hoc temperature from fit_temperature? If it does not, say which loss each one is optimizing and which of the two you would ship.