diff --git a/docs/uncertainty/principled-belief.md b/docs/uncertainty/principled-belief.md index 6f065cbb0..61cd3e863 100644 --- a/docs/uncertainty/principled-belief.md +++ b/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. diff --git a/predicators/code_sim_learning/fit_space.py b/predicators/code_sim_learning/fit_space.py index 8dbf77ff3..f1b51d918 100644 --- a/predicators/code_sim_learning/fit_space.py +++ b/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]: diff --git a/predicators/code_sim_learning/orchestrator.py b/predicators/code_sim_learning/orchestrator.py index 2eb667e2a..7f61a860d 100644 --- a/predicators/code_sim_learning/orchestrator.py +++ b/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: diff --git a/predicators/code_sim_learning/parameter_belief.py b/predicators/code_sim_learning/parameter_belief.py index 37afe6d88..b1aabbfd2 100644 --- a/predicators/code_sim_learning/parameter_belief.py +++ b/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,6 +616,16 @@ 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) @@ -606,6 +633,7 @@ def build_parameter_belief( 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 diff --git a/predicators/code_sim_learning/physical_sysid.py b/predicators/code_sim_learning/physical_sysid.py index 81119c0c8..a65c588d8 100644 --- a/predicators/code_sim_learning/physical_sysid.py +++ b/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 diff --git a/tests/code_sim_learning/test_orchestrator.py b/tests/code_sim_learning/test_orchestrator.py index 1582c6d42..748c10585 100644 --- a/tests/code_sim_learning/test_orchestrator.py +++ b/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,7 +258,8 @@ 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 @@ -264,3 +267,6 @@ def residuals(params): 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) diff --git a/tests/code_sim_learning/test_parameter_belief.py b/tests/code_sim_learning/test_parameter_belief.py index e44a5eb0c..793841535 100644 --- a/tests/code_sim_learning/test_parameter_belief.py +++ b/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 = [ diff --git a/tests/code_sim_learning/test_physical_sysid.py b/tests/code_sim_learning/test_physical_sysid.py index 2b6a20a6d..bc8a30cf7 100644 --- a/tests/code_sim_learning/test_physical_sysid.py +++ b/tests/code_sim_learning/test_physical_sysid.py @@ -1537,7 +1537,8 @@ def fake_lm(_env, _trajs, physical_specs, *_args, **_kwargs): def test_global_multistart_leaves_the_guesses_basin(monkeypatch): """From a guess in the worse basin, the multi-start finds the better one and makes it the point estimate; the worse basin joins the belief only when - it lies within twice the model's own misfit.""" + it lies within twice the model's own misfit, which is returned for the + belief's temperature.""" import dataclasses from predicators.code_sim_learning.config import SysIdConfig @@ -1556,16 +1557,24 @@ def run(noise_sse): # The best fit leaves 0.05 above the noise: the 0.2 basin, 0.2 worse, # is distinguishable. - theta, jac, basins = run(10.25) + theta, jac, basins, misfit = run(10.25) assert theta[0] == 0.8 and jac is not None assert [(b[0]["k"], b[1]) for b in basins] == [(0.8, 10.3)] + assert misfit == pytest.approx(0.05) # The best fit leaves 0.3 above the noise: the data cannot tell the # basins apart, so both are kept, the worse one weighted down. - theta, _jac, basins = run(10.0) + theta, _jac, basins, misfit = run(10.0) assert theta[0] == 0.8 assert [b[0]["k"] for b in basins] == [0.8, 0.2] assert basins[0][2] == 1.0 assert basins[1][2] == pytest.approx(np.exp(-0.2 / 0.6)) + assert misfit == pytest.approx(0.3) + # Without a declared noise channel the misfit is unknown, so the + # point estimate is the only basin. + theta, _jac, basins, misfit = run(0.0) + assert theta[0] == 0.8 + assert [(b[0]["k"], b[2]) for b in basins] == [(0.8, 1.0)] + assert misfit is None def test_global_multistart_runs_for_guessed_starting_values_only(monkeypatch): @@ -1587,7 +1596,7 @@ def fake_lm(_env, def fake_multistart(*args, **_kwargs): calls.append(True) - return args[13][0], args[13][1], [({"mu": 0.5}, 1.0, 1.0)] + return args[13][0], args[13][1], [({"mu": 0.5}, 1.0, 1.0)], 0.2 monkeypatch.setattr(physical_sysid, "fit_map_lm_rollout", fake_lm) monkeypatch.setattr(physical_sysid, "_grid_seed_physical_specs", @@ -1604,6 +1613,33 @@ def fake_multistart(*args, **_kwargs): None, [_trajectory([0.0, 0.1, 0.2])], specs, _RESIDUAL_FEATURES) assert bool(calls) is guesses assert (result.basins is not None) is guesses + assert result.misfit == (0.2 if guesses else None) + + +def test_guessed_fit_measures_its_misfit_without_the_multistart(monkeypatch): + """With the whole-box search off, a fit from guessed starting values still + measures its misfit at the local fit for the belief's temperature.""" + from predicators.settings import CFG + + def fake_lm(_env, _trajs, physical_specs, *_args, **_kwargs): + return np.array([s.init_value for s in physical_specs]), None + + monkeypatch.setattr(physical_sysid, "fit_map_lm_rollout", fake_lm) + monkeypatch.setattr(physical_sysid, "_grid_seed_physical_specs", + lambda _e, _t, specs, *_a, **_k: (list(specs), {})) + monkeypatch.setattr(physical_sysid, "compute_rollout_sse", + lambda *_args, **_kwargs: 5.0) + monkeypatch.setattr(physical_sysid, "expected_noise_sse", + lambda *_args, **_kwargs: 4.0) + monkeypatch.setattr(CFG, "code_sim_learning_prior_spans_bounds", True) + monkeypatch.setattr(CFG, "code_sim_learning_interval_belief", True) + monkeypatch.setattr(CFG, "code_sim_learning_rollout_global_design_points", + 0) + result = physical_sysid.fit_params_rollout( + None, [_trajectory([0.0, 0.1, 0.2])], + [ParamSpec("mu", 0.5, lo=0.1, hi=1.0)], _RESIDUAL_FEATURES) + assert result.basins is None + assert result.misfit == pytest.approx(1.0) def test_anchor_ablation_reverts_compensatory_param():