Count a guessed fit's misfit once in the belief's temperature

docs/uncertainty/principled-belief.md

@@ -28,6 +28,9 @@ There is no added model-error allowance and no between-segment term; both were m
 
 1. Temperature: evaluate the residual vector at the MAP; `lambda = max(1, SSE_min / (N * sigma_n^2))`, where `N` counts residual terms and `sigma_n` is the fit's Gaussian width on scaled residuals.
    It is the maximum-likelihood variance of the standardized errors if they were independent Gaussian; it does not correct for errors that persist through a segment, and the fit report should show the lag-one autocorrelation of the MAP residuals so that gap can be measured.
+   When the multi-start of step 6 has measured the misfit `m`, the temperature counts it as one error that persists through the recordings: `lambda = max(1, SSE_min / (N * sigma_n^2), m / sigma_n^2)`.
+   A replay that drifts off its recording, or a contact that resolves differently, moves many residuals together, and a parameter change can trade that error for a lower SSE, so SSE differences up to `m` are no evidence.
+   It is the scale that weighs the basins, so the basin factors and their weights share one target.
 2. Lines: for each continuous parameter, evaluate the negative log target at that temperature along its coordinate through the MAP in fit space, expanding outward until it rises by the cutoff or reaches the declared bound, then bisecting the interval holding the most mass.
    For a discrete parameter, the line is its set of values: evaluate the target at each value with the other parameters at the MAP, and normalize.
 3. Line posterior: interpolate the negative log density linearly between probes (a piecewise-exponential density), store it on a finer grid, and normalize; a flat line leaves the prior.

predicators/code_sim_learning/fit_space.py

@@ -123,6 +123,12 @@ class FitResult:
     # mixture weight); the first is the point estimate. The parameter
     # belief mixes them. None = no multi-start ran.
     basins: Optional[List[Tuple[Dict[str, float], float, float]]] = None
+    # The SSE the multi-start's best fit leaves above the declared noise
+    # channel's expected SSE (at least the likelihood floor). It weighs
+    # the basins, and the parameter belief's temperature counts it as one
+    # error that persists through the recordings. None = no multi-start
+    # ran.
+    misfit: Optional[float] = None
 
     @property
     def point_estimate(self) -> Dict[str, float]:

predicators/code_sim_learning/orchestrator.py

@@ -480,7 +480,9 @@ def _build_belief(
     anchor (else its declared init) with the fit's widths, recomputed
     the same way when the result does not carry them. Draws are seeded
     by the global seed and the fit's identity, so a cached fit reuses
-    its draws.
+    its draws. When the fit measured its misfit (the global multi-start
+    did), every basin's factor counts it as one persistent error, the
+    same scale that weighs the basins.
     """
     # pylint: disable-next=import-outside-toplevel
     from predicators.settings import CFG
@@ -508,7 +510,8 @@ def around(point: Dict[str, float], config: BeliefConfig,
             config=config,
             seed=stable_seed(CFG.seed, sorted(point.items()), num_survivors,
                              *seed_parts),
-            batch_residuals=batch_residuals_fn)
+            batch_residuals=batch_residuals_fn,
+            misfit=float(result.misfit or 0.0))
 
     basins = list(result.basins or [])
     if len(basins) <= 1:

predicators/code_sim_learning/parameter_belief.py

@@ -33,7 +33,7 @@
 
 import hashlib
 import logging
-from dataclasses import dataclass, field
+from dataclasses import dataclass, field, replace
 from typing import Any, Callable, Dict, Generator, Iterator, List, Optional, \
     Sequence, Tuple
 
@@ -273,6 +273,10 @@ class ParameterBelief:
     # components' weighted mixtures. Empty for a single basin.
     components: List["ParameterBelief"] = field(default_factory=list)
     weights: List[float] = field(default_factory=list)
+    # The fit's misfit (SSE above the declared noise's expected SSE)
+    # when the temperature counts it as one persistent error
+    # (build_parameter_belief); 0 when it does not.
+    misfit: float = 0.0
 
     @property
     def num_draws(self) -> int:
@@ -379,7 +383,15 @@ def describe(self) -> List[str]:
                 "parameter sets about as well as the model fits them at "
                 "all; the belief mixes them, and each draw comes from one "
                 "of them.")
-        if self.noise_scale > 1.0:
+        if self.noise_scale > 1.0 and self.misfit > 0.0:
+            lines.append(
+                f"The program's best fit leaves an SSE of {self.misfit:.3g} "
+                "above what the declared noise explains, and that error "
+                "persists through each recording, so the posterior counts it "
+                "as one error: parameter values whose SSE is within about "
+                "that much of the best fit's are about as likely "
+                f"(temperature {self.noise_scale:.3g}).")
+        elif self.noise_scale > 1.0:
             lines.append(
                 f"The program misfits the recordings {self.noise_scale:.3g}x "
                 "beyond the declared noise, so the posterior uses that "
@@ -425,6 +437,7 @@ def to_dict(self) -> Dict[str, Any]:
             },
             "components": [c.to_dict() for c in self.components],
             "weights": list(self.weights),
+            "misfit": self.misfit,
         }
 
     @classmethod
@@ -450,7 +463,8 @@ def from_dict(cls, data: Dict[str, Any]) -> ParameterBelief:
                 for n, (v, p) in data.get("discrete", {}).items()
             },
             components=[cls.from_dict(c) for c in data.get("components", [])],
-            weights=[float(w) for w in data.get("weights", [])])
+            weights=[float(w) for w in data.get("weights", [])],
+            misfit=float(data.get("misfit", 0.0)))
 
 
 def _mix_lines(lines: Sequence[LinePosterior],
@@ -531,7 +545,8 @@ def mix_beliefs(components: Sequence[ParameterBelief],
                            evaluations=sum(c.evaluations for c in components),
                            discrete=discrete,
                            components=list(components),
-                           weights=shares)
+                           weights=shares,
+                           misfit=first.misfit)
 
 
 def join_beliefs(first: ParameterBelief,
@@ -570,7 +585,8 @@ def join_beliefs(first: ParameterBelief,
             **second.discrete
         },
         components=[join_beliefs(c, second) for c in first.components],
-        weights=list(first.weights))
+        weights=list(first.weights),
+        misfit=first.misfit)
 
 
 def stable_seed(*parts: Any) -> int:
@@ -589,6 +605,7 @@ def build_parameter_belief(
     config: BeliefConfig,
     seed: int,
     batch_residuals: Optional[BatchResidualsFn] = None,
+    misfit: float = 0.0,
 ) -> ParameterBelief:
     """Build ``q(theta)`` around ``map_params`` and draw from it.
 
@@ -599,13 +616,24 @@ def build_parameter_belief(
     a list of parameter dicts together, each as ``residuals`` would; the
     lines are then traced in lockstep (:func:`_trace_lines`), with the
     same result.
+
+    The temperature is the estimated noise level, at least 1: ``SSE_min
+    / (N noise_sigma^2)`` spreads the misfit over the ``N`` residuals as
+    independent errors. ``misfit`` (SSE units), when given, is the fit's
+    error above the declared noise counted as ONE error that persists
+    through the recordings: a replay that drifts off its recording, or a
+    contact that resolves differently, moves many residuals together,
+    and a parameter change can trade that error for a lower SSE. SSE
+    differences up to the misfit are then no evidence, so the
+    temperature is at least ``misfit / noise_sigma^2``.
     """
     assert noise_sigma > 0.0, noise_sigma
     at_map = np.asarray(residuals(dict(map_params)), dtype=float)
     sse_min = float(np.dot(at_map, at_map))
     noise_scale = 1.0
     if at_map.size:
         noise_scale = max(1.0, sse_min / (at_map.size * noise_sigma**2))
+    noise_scale = max(noise_scale, misfit / noise_sigma**2)
     # Negative log likelihood per unit of summed squared residual, at the
     # estimated noise level.
     k = 1.0 / (2.0 * noise_scale * noise_sigma**2)
@@ -648,12 +676,13 @@ def build_parameter_belief(
         lines[spec.name] = LinePosterior.from_neg_log(grid,
                                                       [line(z) for z in grid])
 
-    belief = _with_draws(specs, map_params, lines, held, noise_scale,
-                         evaluations, config, seed, discrete)
+    belief = replace(_with_draws(specs, map_params, lines, held, noise_scale,
+                                 evaluations, config, seed, discrete),
+                     misfit=max(misfit, 0.0))
     logger.info(
         "Parameter belief: %d draws over %d parameters from %d line "
-        "evaluations (noise level %.3gx the declared one).", belief.num_draws,
-        len(specs), evaluations, noise_scale)
+        "evaluations (temperature %.3g).", belief.num_draws, len(specs),
+        evaluations, noise_scale)
     return belief
 
 

predicators/code_sim_learning/physical_sysid.py

@@ -334,13 +334,21 @@ def fit_params_rollout(
                 "sensitive": False,
             }
     basins: Optional[List[Tuple[Dict[str, float], float, float]]] = None
+    misfit: Optional[float] = None
     if (config.anchors_are_guesses and config.global_design_points > 0
             and config.global_starts > 0 and trajectories):
-        lm_theta, lm_jac, basins = _global_multistart(
+        lm_theta, lm_jac, basins, misfit = _global_multistart(
             base_env, trajectories, all_specs, physical_specs, rule_specs,
             residual_features, rules, latent_init, scaling, center_int,
             prior_sigma, noise_sigma, config, (lm_theta, lm_jac), noise_sse,
             sigma_tol)
+    elif config.anchors_are_guesses and trajectories and noise_sse > 0.0:
+        # The same misfit the multi-start measures, at the local fit.
+        map_sse = compute_rollout_sse(
+            base_env, trajectories,
+            dict(zip(names, (float(v) for v in lm_theta))), residual_features,
+            physical_names, rules, latent_init, scaling)
+        misfit = max(map_sse - noise_sse, sigma_tol)
     result = FitResult(names=names,
                        samples=np.asarray(lm_theta, dtype=float)[None, :],
                        log_probs=np.zeros(1),
@@ -350,7 +358,8 @@ def fit_params_rollout(
                        scales=scales,
                        sensitivity=sensitivity,
                        lm_notes=lm_notes,
-                       basins=basins)
+                       basins=basins,
+                       misfit=misfit)
     n_lm = num_rollouts_run() - n_start - n_grid
     # The ablation pins moved parameters back to their anchors, which
     # only calibrated baselines warrant; a guess has no claim beyond the
@@ -407,8 +416,8 @@ def _global_multistart(
     incumbent: Tuple[np.ndarray, Optional[np.ndarray]],
     noise_sse: float,
     sigma_tol: float,
-) -> Tuple[np.ndarray, Optional[np.ndarray], List[Tuple[Dict[str, float],
-                                                        float, float]]]:
+) -> Tuple[np.ndarray, Optional[np.ndarray], List[Tuple[Dict[
+        str, float], float, float]], Optional[float]]:
     """Search the whole box for the fit when the starting values are guesses.
 
     A scrambled Sobol design of ``config.global_design_points`` over
@@ -425,7 +434,10 @@ def _global_multistart(
     leaves above the declared noise channel's expected SSE, or the
     likelihood floor when that is larger. Recordings the model fits no
     closer than that cannot prefer one of them, so the parameter belief
-    mixes them, each weighted by ``exp(-dSSE / (2 misfit))``.
+    mixes them, each weighted by ``exp(-dSSE / (2 misfit))``. The misfit
+    is returned too: the belief's temperature counts it the same way.
+    Without a declared noise channel (``noise_sse`` 0) the misfit is
+    unknown: the point estimate is the only basin and the misfit None.
     """
     names = [s.name for s in all_specs]
     physical_names = [s.name for s in physical_specs]
@@ -486,13 +498,16 @@ def started_at(specs: Sequence[ParamSpec],
     ]
     ranked = [int(i) for i in np.argsort(np.asarray(costs), kind="stable")]
     best = ranked[0]
-    misfit = max(finals[best] - noise_sse, sigma_tol)
-    band = 2.0 * misfit
+    misfit: Optional[float] = None
+    band = 0.0
+    if noise_sse > 0.0:
+        misfit = max(finals[best] - noise_sse, sigma_tol)
+        band = 2.0 * misfit
     width = hi - lo
     basins: List[Tuple[Dict[str, float], float, float]] = []
     kept: List[np.ndarray] = []
     for i in ranked:
-        if finals[i] > finals[best] + band:
+        if finals[i] > finals[best] + band or (misfit is None and i != best):
             continue
         z = to_fit_space(list(all_specs), thetas[i])
         if any(
@@ -501,7 +516,7 @@ def started_at(specs: Sequence[ParamSpec],
             continue
         kept.append(z)
         gap = finals[i] - finals[best]
-        weight = float(np.exp(-gap / (2.0 * misfit))) if misfit > 0 else 1.0
+        weight = float(np.exp(-gap / (2.0 * misfit))) if misfit else 1.0
         basins.append((points[i], float(finals[i]), weight))
     logger.info(
         "Rollout sysID global multi-start: %d design points (best SSE %.4g), "
@@ -510,7 +525,7 @@ def started_at(specs: Sequence[ParamSpec],
         float(np.min(sses)),
         len(thetas) - 1, finals[0], [round(float(s), 4) for s in finals[1:]],
         len(basins), band)
-    return thetas[best], jacobians[best], basins
+    return thetas[best], jacobians[best], basins, misfit
 
 
 # A MAP within this fit-space distance of its anchor counts as unmoved

tests/code_sim_learning/test_orchestrator.py

@@ -9,6 +9,7 @@
 
 import numpy as np
 import pybullet as p
+import pytest
 
 import predicators.approaches  # noqa: F401  # pylint: disable=unused-import
 from predicators import utils
@@ -238,7 +239,8 @@ def test_draw_shares_follow_the_weights():
 
 def test_belief_mixes_the_fit_basins():
     """A fit with two basins gets a belief over both, its draws split by
-    weight, and the point estimate stays the best basin's."""
+    weight, and the point estimate stays the best basin's; every basin's factor
+    takes the temperature of the misfit that weighed them."""
     spec = ParamSpec("k", 0.5, lo=0.0, hi=1.0)
 
     def residuals(params):
@@ -256,11 +258,15 @@ def residuals(params):
                        }, 10.3, 1.0),
                                ({
                                    "k": 0.2
-                               }, 10.5, float(np.exp(-1.0 / 3.0)))])
+                               }, 10.5, float(np.exp(-1.0 / 3.0)))],
+                       misfit=0.3)
     belief = _build_belief(result, [spec], {"k": 0.5}, residuals,
                            BeliefConfig(num_draws=16, line_max_evals=40), 1)
     assert len(belief.components) == 2
     assert belief.num_draws == 16
     assert belief.map_estimate == {"k": 0.8}
     values = np.sort(belief.draws[:, 0])
     assert values[0] < 0.5 < values[-1]
+    assert [c.noise_scale for c in belief.components] == \
+        pytest.approx([0.3 / 0.05**2] * 2)
+    assert belief.misfit == pytest.approx(0.3)

tests/code_sim_learning/test_parameter_belief.py

@@ -29,7 +29,13 @@ def residuals(params):
     return residuals
 
 
-def _build(specs, map_params, residuals, prior_sigma=1e6, config=None, seed=0):
+def _build(specs,
+           map_params,
+           residuals,
+           prior_sigma=1e6,
+           config=None,
+           seed=0,
+           misfit=0.0):
     return build_parameter_belief(
         specs,
         map_params,
@@ -40,7 +46,8 @@ def _build(specs, map_params, residuals, prior_sigma=1e6, config=None, seed=0):
         prior_sigmas={s.name: prior_sigma
                       for s in specs},
         config=config or BeliefConfig(num_draws=64, line_max_evals=40),
-        seed=seed)
+        seed=seed,
+        misfit=misfit)
 
 
 def test_line_posterior_moments_and_sampling():
@@ -251,6 +258,39 @@ def test_misfit_scales_the_noise_level():
     assert ratio == pytest.approx(np.sqrt(27.0 / 4.0), rel=0.06)
 
 
+def test_a_persistent_misfit_counts_as_one_error():
+    """A misfit the fit measured counts once: the temperature is the misfit
+    over sigma_n^2, where spreading the same 27 sigma_n^2 over four residual
+    terms gives 27/4.
+
+    Mixing, joining and checkpoints keep it, and the report says what it
+    means.
+    """
+    specs = [ParamSpec("a", 0.3, lo=-10, hi=10)]
+    residuals = _gaussian_residuals(["a"], [0.3], [[20.0]],
+                                    offset=[3.0 * _SIGMA_N] * 3)
+    spread = _build(specs, {"a": 0.3}, residuals)
+    once = _build(specs, {"a": 0.3}, residuals, misfit=27.0 * _SIGMA_N**2)
+    assert once.noise_scale == pytest.approx(27.0)
+    assert once.misfit == pytest.approx(27.0 * _SIGMA_N**2)
+    ratio = np.sqrt(once.lines["a"].variance() / spread.lines["a"].variance())
+    assert ratio == pytest.approx(2.0, rel=0.06)
+    assert "persists through each recording" in "\n".join(once.describe())
+    assert "persists" not in "\n".join(spread.describe())
+    mixed = mix_beliefs([once, once], [1.0, 1.0])
+    assert ParameterBelief.from_dict(mixed.to_dict()).misfit == once.misfit
+    extra = prior_belief([ParamSpec("c", 0.5, lo=0.0, hi=1.0)], {"c": 0.0},
+                         {"c": 1.0},
+                         BeliefConfig(num_draws=128),
+                         seed=0)
+    assert join_beliefs(mixed, extra).misfit == once.misfit
+    # A misfit within the noise leaves the ordinary posterior.
+    fits = _gaussian_residuals(["a"], [0.3], [[20.0]])
+    assert _build(specs, {
+        "a": 0.3
+    }, fits, misfit=_SIGMA_N**2).noise_scale == pytest.approx(1.0)
+
+
 def test_draws_are_seeded_bounded_and_sample_discrete_parameters():
     """Draws respect the box, the scale, discreteness and the seed."""
     specs = [