`sue-probit-msa` — Probit SUE by MSA, certified by pinned Monte Carlo¶

What. Probit SUE draws perceived route costs from correlated normals instead of the logit’s independent Gumbels, so — unlike logit — it has no closed-form loading. sue-probit-msa reaches the probit fixed point by MSA over a Monte-Carlo STOCH loading, and the harness certifies it with a pinned Monte-Carlo residual (fixed sample count r_cert), reported with a standard error and a noise floor ([daganzo1977stochastic], ADR-003, docs/REFERENCES.md).

Why it is in the benchmark. It is the benchmark’s stochastic-certificate case: the scored number is itself an estimate, so it carries uncertainty. See the model compendium and docs/ARCHITECTURE.md (P1).

Scope. Runs on the built-in two-route probit anchor (β = 0.1). The budget is capped at 50 iterations for tutorial runtime — so the run is deliberately short of the fixed point, and the certificate says so honestly (residual above the MC floor).

How this notebook is graded¶

A notebook never claims a number it does not compute in that cell. Every scored quantity below is recomputed live by the P1 Evaluator from the flows the model emitted, in the cell where it is claimed. Model self-reports are shown only as provenance and diffed against the certificate as an honesty check, exactly as the harness treats them (README, Certified, not self-reported).

# Setup. `sue-probit-msa` is a core model: a plain `pip install -e .` suffices — no
# optional extra, so no guard cell. The inline backend is Agg-based (headless CI
# renders into the notebook); NEVER matplotlib.use("Agg") in-kernel — it silently
# suppresses inline figure capture.
%matplotlib inline
import numpy as np

from tabench import (
    Budget,
    Evaluator,
    RngBundle,
    SueProbitMsaModel,
    Trace,
    two_route_scenario,
    viz,
)

The scenario¶

The same two disjoint routes as the logit anchor, but the probit task: perceived route costs are correlated normals with perception variance β = 0.1 per unit free-flow time. A different equilibrium definition → a different content hash (P2).

scenario = two_route_scenario(sue_theta=0.1, sue_family="probit")
net = scenario.network

print(f"scenario      : {scenario.name}")
print(f"content hash  : {scenario.content_hash()[:16]}…")
print(f"SUE family    : {scenario.sue_family} (beta = {scenario.sue_theta})")
print(f"total demand  : {scenario.demand.total}")

scenario      : tworoute
content hash  : ab27373c328c9ecc…
SUE family    : probit (beta = 0.1)
total demand  : 4.0

Solve¶

The model contract (CONTRIBUTING.md): a model receives (scenario, budget, rng, trace), records checkpoints, and respects the budget. Budgets are hardware-free (iterations / shortest-path calls; wall-clock is recorded but never the ranking axis, P7). Whatever the model writes into self_report is provenance, not a score.

# Budget capped at 50 iterations: every certified checkpoint costs a full Monte-Carlo
# loading, so this keeps the tutorial fast (the model is then short of the fixed point).
model = SueProbitMsaModel()
bundle = model.solve(scenario, Budget(iterations=50), RngBundle(0), Trace())

final = bundle.final
print(f"model              : {model.name}")
print(f"budget spent       : {final.coords.iterations} iterations")
print(f"emitted flows      : {np.round(final.link_flows, 4)}")
print(f"self-reported resid: {final.self_report['sue_fixed_point_residual']:.4f}  (provenance only)")

model              : sue-probit-msa
budget spent       : 50 iterations
emitted flows      : [2.6122 2.6122 1.3878 1.3878]
self-reported resid: 2.6122  (provenance only)

Certify (P1)¶

Probit SUE has no closed form, so the harness certifies with a pinned Monte-Carlo loading (r_cert = 2000 samples): the residual is an estimate with a standard error se and a noise floor below which it is indistinguishable from zero (ADR-003). The model’s own residual estimate uses different draws, so it legitimately DIFFERS from the certificate — the sharpest form of never trust the self-report. We recompute the analytic probit fixed point in-cell with brentq.

from scipy.optimize import brentq
from scipy.stats import norm

# Pinned Monte-Carlo certificate: r_cert samples, root_seed fixes the draws (ADR-003).
evaluator = Evaluator(scenario, root_seed=0, r_cert=2000)
metrics = evaluator.evaluate(final.link_flows)
residual = metrics["sue_fixed_point_residual"]
se = metrics["sue_residual_se"]
floor = metrics["sue_residual_floor"]
print(f"certified probit residual : {residual:.4f}  ± {se:.4f} (MC se)")
print(f"MC noise floor            : {floor:.4f}")
print(f"feasible                  : {metrics['feasible']:.0f}")

assert metrics["feasible"] == 1.0
assert np.isfinite(residual) and residual >= 0.0

# NEVER trust the self-report (P1): the model's own residual estimate differs from the
# pinned-MC certificate — different draws, different counts. The certificate is scored.
self_resid = final.self_report["sue_fixed_point_residual"]
print(f"model self-report resid   : {self_resid:.4f}  (provenance only — NOT the score)")

# Analytic anchor RECOMPUTED: disjoint routes -> independent normal route costs ->
# P(A) = Phi((c_B - c_A) / sqrt(3.5 beta)); the fixed point via brentq.
demand = scenario.demand.total
beta = scenario.sue_theta

def _probit_residual(f_a):
    c_a = 2.0 + f_a
    c_b = 1.5 + 2.0 * (demand - f_a)
    return f_a - demand * norm.cdf((c_b - c_a) / np.sqrt(3.5 * beta))

f_a = brentq(_probit_residual, 1e-9, demand - 1e-9)
ref_flows = np.array([f_a, f_a, demand - f_a, demand - f_a])
print(f"probit fixed point f_A    : {f_a:.4f}  (recomputed);  emitted f_A: {final.link_flows[0]:.4f}")

# At this capped budget the run is short of the fixed point: the certified residual
# sits ABOVE the MC floor, and the certificate reports that honestly.
assert residual > floor
assert 0.0 < f_a < demand

certified probit residual : 1.4862  ± 0.0402 (MC se)
MC noise floor            : 0.0321
feasible                  : 1
model self-report resid   : 2.6122  (provenance only — NOT the score)
probit fixed point f_A    : 2.4444  (recomputed);  emitted f_A: 2.6122

Visualize¶

Both figures come from tabench.viz. Left/top: the emitted probit-SUE flows on the two-route network. Right/bottom: the emitted flows against the brentq-recomputed probit fixed point — points sit OFF the y = x guide because the 50-iteration budget stopped short of the fixed point, exactly what the certified residual reports.

display(viz.plot_network_flows(net, final.link_flows))
display(viz.plot_flow_scatter(("probit SUE (brentq)", ref_flows), {"sue-probit-msa": final.link_flows}))

../../_images/6e4283c84bb44d33e38c9b176c5813659da904d0d57a594e3cec38574efbc331.png

../../_images/c1f51dc6da9b8b2806b0b8e783be48ed16b98bee19762e4fba2dc9a55d8cd555.png

Takeaways & pointers¶

The score is an estimate. The probit residual is Monte-Carlo — the harness reports it with a standard error and a floor; it never claims more precision than the sampling affords (ADR-003).
Self-report ≠ certificate. The model’s own residual differs from the pinned-MC certificate; only the certificate is scored.
Honest about budget. Capped at 50 iterations, the run is short of the fixed point and the residual sits above the floor — reported, not hidden.
Where next. The closed-form logit variant: sue-msa; ADR-003 in the model compendium.

sue-probit-msa — Probit SUE by MSA, certified by pinned Monte Carlo¶