wdm0006/python-skills/threshold-calibration

Calibrating Thresholds and Baselines

A scoring pipeline has three layers: features computed per item, reference constants the features are compared against, and a flag/score layer that turns the comparison into something a user reads. reporting-derived-metrics covers the first layer — what a feature function returns when the sample is too thin to support it. This skill covers the other two: where the constant came from, and whether the flag it drives carries any information at all.

The failure mode is not a crash. It is a pipeline that runs clean, emits confident-looking flags on every input, and nobody notices because no one ever asked how often each flag fires per class.

A committed baseline file is not a measurement

Baselines get checked in as JSON or a module constant and immediately start reading as ground truth, because they have the shape of data:

{
  "name": "general_prose",
  "created": "2023-09-04",
  "sample_size": 750,
  "features": {
    "avg_item_len": {"mean": 12.4, "sd": 5.1},
    "diversity":    {"mean": 0.61, "sd": 0.07},
    "rate_per_word":{"mean": 0.22, "sd": 0.05}
  }
}

Every number there is two significant figures, sample_size is a round 750, and nothing in the repository produces the file. It was typed. Downstream, z = (value - 0.22) / 0.05 is reported to four decimals and quoted in an issue as evidence.

Before trusting a stored baseline, ask four questions:

• Is there a script in the tree that regenerates it? If not, the constants

are an author's estimate. Say so where they are defined, in a comment the next reader cannot miss.

• What precision do they claim? A sd of 0.05 quoted to two significant

figures supports roughly one digit of z-score. Do not print z = -3.5814.

• Is there a second copy? Baselines are routinely duplicated as an in-code

fallback next to the file they mirror. Two copies drift; make one import the other, or add a test that asserts they are equal.

• What population is it supposed to describe, and is it the population you

are scoring? A baseline built from long-form documents applied to 150-word snippets is a different distribution, not a stricter one.

When you do regenerate a baseline from real data, commit the extraction script alongside it and record the corpus, the filter, and the date in the file.

The feature and the constant must be the same quantity

This is the highest-yield bug in the whole pipeline, and it does not look like a bug — the arithmetic is correct end to end, the units just disagree.

# Feature: marks per CHARACTER.
def mark_rate(text: str) -> float:
    return sum(c in MARKS for c in text) / len(text)
# Baseline: {"mean": 0.22, "sd": 0.05} — authored as marks per WORD.
z = (mark_rate(text) - baseline["mean"]) / baseline["sd"]

A per-character rate lands near 0.041 where a per-word rate lands near 0.219. The z-score is then about -3.6 for essentially every input — against -0.02 for the quantity the constant actually describes — and the feature is reported as an extreme outlier on 100% of items. The same defect arrives by a second route — length confounding — where the units match but the normalization does not: a per-item rate that falls as items get longer (distinct values over total values, unique-token rates, any "new things per thing" measure) compared against a constant derived from documents an order of magnitude longer will sit far from the mean for reasons that have nothing to do with the property being detected.

The shared tell is the same in both cases, and it is not a weird number — it is a z-score that is large and nearly constant across every input. A feature that genuinely discriminates has a z-score that moves. Before shipping any new feature-plus-baseline pair, print the distribution of its z-score over a sample of known-normal items and confirm it is centered near zero:

zs = [z_score(feature(item), baseline) for item in known_normal_sample]
print(f"n={len(zs)} median={statistics.median(zs):+.2f} "
      f"min={min(zs):+.2f} max={max(zs):+.2f}")
# median well away from 0 with a narrow spread => units or population mismatch,
# not a detection.

Per-class fire rates are the cheapest calibration probe you have

You do not need a tuned model, a new run, or a metrics framework to find a worthless flag. If the pipeline persists per-item output, the fire rate of every flag per class is a few lines over a file you already have:

from collections import Counter
fires = {"known_normal": Counter(), "known_anomalous": Counter()}
for row in load_scored_items():          # one persisted record per item
    for name in row["warnings"] + row["errors"] + row["indicators"]:
        fires[row["label"]][name] += 1

Read the result as a table of fired / total per class, and apply one rule:

A flag that fires on nearly every item of both classes carries zero information, no matter how severe its name is. rate_per_word at 487/500 normal and 463/500 anomalous is not a strict check; it is a constant, and its presence in every report trains readers to ignore the whole list. The same is true in reverse: a flag that fires 0/1000 has never been evidence for anything, and its threshold has no support.

Two practical notes:

• Read the persisted record's real shape, not the API's. The serialized row

and the public result object often nest differently — flags at the top level of the record but under a flags sub-dict in the response. A lookup against the wrong shape returns nothing and reads exactly like "no flags fired." Print one record's keys before writing the loop.

• Surface these rates in the generated report. If the fire rates only exist

in an ad-hoc script, the next person tuning a threshold will not have them.

Re-render reports without re-scoring the sample

A reporting-only change should consume the persisted scores and their original run metadata. Re-running the full benchmark can change feature values because the scoring code, dependency versions, or model runtime have moved since the recorded run — even when no final classification changes. That mixes a new measurement with a presentation change and makes the report diff hard to assess.

Use the report renderer directly, then verify that its inputs stayed byte-for-byte unchanged. For example, around the project's report-rendering command:

from pathlib import Path
inputs = [Path("results/scores.jsonl"), Path("results/run_metadata.json")]
before = {path: path.read_bytes() for path in inputs}
render_report(scores_path=inputs[0], metadata_path=inputs[1], output="report.md")
assert all(path.read_bytes() == contents for path, contents in before.items())

If a full run already overwrote those files, restore only its generated scores and run metadata to the intended recorded version, then re-render. Keep intentional re-scoring in a separate change with its own runtime and dependency provenance; do not demand that today's full benchmark reproduce old scores unchanged.

Score composition: count each family once

A confidence score assembled by addition quietly weights whichever feature family has the most members:

# Bad — one anomalous category adds 0.1, but the *reason* is appended once,
# so eight anomalous categories silently contribute 0.8 under one label.
for category in categories:
    if is_anomalous(category):
        confidence += 0.1
if any(is_anomalous(c) for c in categories):
    indicators.append("category_anomalies")

The report then shows a single indicator next to a near-maximum score, and the contribution is invisible to anyone reading it. Cap each family's contribution, or emit one reason per increment so the arithmetic in the score matches the narrative beside it. Log the per-indicator contributions when the score is built; a score you cannot decompose cannot be tuned.

Verify saved scores before replaying a flag-layer change

Saved features and z-scores let you compare score-composition changes without re-running expensive feature extraction. But saved flags may come from an older implementation than the base you are comparing against. First replay the base flag function and compare it with the saved output, including reasons and confidence, not just the final boolean:

mismatches = []
for row in scored_rows:
    replayed = base_generate_flags(row["features"], row["z_scores"])
    if replayed != row["flags"]:
        mismatches.append(row["id"])
assert not mismatches, f"saved flags differ from base replay: {mismatches}"

If this check passes, replay the changed function on the same inputs. If it fails, do not compare stored confidence values against newly computed ones and attribute the entire delta to your change. Generate fresh base and changed results with the same corpus, feature extractor, dependencies, and configuration, varying only the flag implementation; keep both outputs in scratch directories. For a score-layer-only change, run only the affected analysis rather than loading unrelated expensive models. Report the comparison's actual provenance, and leave committed benchmark results alone unless refreshing them is part of the change.

Sweep every statistic a threshold reads, in both directions

A threshold sweep that hardcodes its comparison cannot evaluate half your thresholds:

# Bad — only expresses "high is anomalous".
def sweep(scores, labels, thresholds):
    return [(t, confusion(scores, labels, lambda s: s >= t)) for t in thresholds]
# Good — direction is an input, because some statistics are low-is-anomalous.
def sweep(scores, labels, thresholds, direction="high"):
    cmp = (operator.ge if direction == "high" else operator.le)
    return [(t, confusion(scores, labels, lambda s: cmp(s, t))) for t in thresholds]

And sweep every statistic a threshold reads, not only the headline score. A maximum-value cutoff that fires on 0 of 1000 items is usually reported as "the threshold is too loose" — but with no sweep over that statistic there is no evidence about whether any cutoff separates the classes on it. Those are very different findings, and only one of them justifies changing the number.

Check the combined rule, not only its cutoffs

A rule can be inactive even when one of its inputs separates the classes:

flag = score <= score_max and variation < variation_min

If variation_min is below every observed variation, the second predicate is false for every item. No change to score_max can make this rule fire. A zero fire rate does not establish that the score is useless; it may identify a blocking predicate in the conjunction.

On the same labeled rows, report each predicate's fire rate and the combined rule's fire rate, separately for each class. Include the observed range and the configured cutoff for each input. Then sweep each statistic alone in its actual comparison direction before evaluating candidate cutoffs together:

for label in ("normal", "anomalous"):
    rows = [r for r in scored_rows if r["label"] == label]
    assert rows, f"no scored rows for {label}"
    score_hits = [r["score"] <= score_max for r in rows]
    variation_hits = [r["variation"] < variation_min for r in rows]
    joint_hits = [a and b for a, b in zip(score_hits, variation_hits)]
    print(label, len(rows), sum(score_hits), sum(variation_hits), sum(joint_hits))

Use the same rows with both measurements available for this comparison, and report excluded rows rather than treating missing measurements as zero. A useful single-statistic cutoff is evidence to investigate the compound rule, not a reason to silently remove its other condition. Report TPR and FPR for the final rule; label cutoffs selected on the evaluation sample as exploratory rather than claiming validated performance on that same sample.

Know whether the feature carries the decision or only the warning list