diff --git a/crates/navigator-analysis/examples/archaic_classify_dump.rs b/crates/navigator-analysis/examples/archaic_classify_dump.rs
new file mode 100644
index 0000000..2150622
--- /dev/null
+++ b/crates/navigator-analysis/examples/archaic_classify_dump.rs
@@ -0,0 +1,38 @@
+//! Dump the archaic **diagnostic** sites (position, derived base, lineage class) as TSV.
+//!
+//! Written to test a different observable for the Tier B HMM. The current model counts *all*
+//! private variants per window, and that signal is weak: measured on a real European, archaic
+//! tracts carry only 2.89x the background density while the background itself varies 5.3x between
+//! its 10th and 90th percentile. Restricting the observable to sites where the derived allele is
+//! actually known to be archaic should be far more specific.
+//!
+//! ```sh
+//! cargo run --release -p navigator-analysis --example archaic_classify_dump -- \
+//!   ~/.decodingus/ancestry/archaic_classify_chm13v2.0.bin chr21 chr22 > classify.tsv
+//! ```
+
+use navigator_analysis::archaic::ArchaicClassify;
+
+fn main() -> Result<(), Box<dyn std::error::Error>> {
+    let mut a = std::env::args().skip(1);
+    let path = a.next().expect("usage: archaic_classify_dump <classify.bin> [contig ...]");
+    let want: Vec<String> = a.collect();
+
+    let cls = ArchaicClassify::from_bytes(&std::fs::read(&path)?).map_err(|e| e.to_string())?;
+    println!("contig\tposition\tderived\tclass");
+    for c in &cls.contigs {
+        let name = &c.positions.contig;
+        if !want.is_empty() && !want.contains(name) {
+            continue;
+        }
+        let mut n = 0usize;
+        for (i, p) in c.positions.iter().enumerate() {
+            let d = c.derived.get(i).copied().unwrap_or(b'N') as char;
+            let k = c.classes.get(i).copied().unwrap_or(2);
+            println!("{name}\t{p}\t{d}\t{k}");
+            n += 1;
+        }
+        eprintln!("{name}: {n} diagnostic sites");
+    }
+    Ok(())
+}
diff --git a/crates/navigator-analysis/examples/archaic_outgroup_density.rs b/crates/navigator-analysis/examples/archaic_outgroup_density.rs
new file mode 100644
index 0000000..78e356e
--- /dev/null
+++ b/crates/navigator-analysis/examples/archaic_outgroup_density.rs
@@ -0,0 +1,45 @@
+//! Per-window counts of African-outgroup segregating sites — a candidate local mutation-rate proxy.
+//!
+//! The Tier B emission model assumes one background rate genome-wide. Measured, the background
+//! private-variant density varies 5.3x between its 10th and 90th percentile and is 14.6x
+//! overdispersed relative to the Poisson it is modelled with, which is larger than the 2.89x
+//! enrichment inside real archaic tracts — so the model calls its own upper tail archaic. hmmix
+//! avoids this with a mutation-rate map; we have no such asset.
+//!
+//! The density of sites segregating in Africans is already in `archaic_outgroup_af_<build>.bin` and
+//! is a direct measure of how variable a region is, for reasons that have nothing to do with
+//! archaic introgression (mutation rate, reference quality, mappability). This dumps it so that
+//! proxy can be tested as a normalizer before an asset is built for the purpose.
+//!
+//! ```sh
+//! cargo run --release -p navigator-analysis --example archaic_outgroup_density -- \
+//!   ~/.decodingus/ancestry/archaic_outgroup_af_chm13v2.0.bin 1000 chr21 chr22 > og_density.tsv
+//! ```
+
+use navigator_analysis::archaic::ArchaicOutgroup;
+
+fn main() -> Result<(), Box<dyn std::error::Error>> {
+    let mut a = std::env::args().skip(1);
+    let path = a.next().expect("usage: archaic_outgroup_density <outgroup.bin> [window_bp] [contig ...]");
+    let window: i64 = a.next().and_then(|s| s.parse().ok()).unwrap_or(1000);
+    let want: Vec<String> = a.collect();
+
+    let og = ArchaicOutgroup::from_bytes(&std::fs::read(&path)?).map_err(|e| e.to_string())?;
+    println!("contig\twindow_start\tn_outgroup_sites");
+    for c in &og.contigs {
+        if !want.is_empty() && !want.contains(&c.contig) {
+            continue;
+        }
+        let mut counts: std::collections::BTreeMap<i64, u32> = Default::default();
+        let mut n = 0u64;
+        for p in c.iter() {
+            *counts.entry(p / window * window).or_insert(0) += 1;
+            n += 1;
+        }
+        eprintln!("{}: {n} outgroup sites in {} non-empty windows", c.contig, counts.len());
+        for (w, k) in counts {
+            println!("{}\t{}\t{}", c.contig, w, k);
+        }
+    }
+    Ok(())
+}
diff --git a/crates/navigator-analysis/examples/archaic_private_dump.rs b/crates/navigator-analysis/examples/archaic_private_dump.rs
new file mode 100644
index 0000000..86c573b
--- /dev/null
+++ b/crates/navigator-analysis/examples/archaic_private_dump.rs
@@ -0,0 +1,50 @@
+//! Dump the **private** variant positions the Tier B HMM actually sees — the subject's derived
+//! variants after the African-outgroup strip — so the input to the model can be checked against an
+//! external truth set independently of the model.
+//!
+//! The segment caller is a density model over exactly these positions. If they are not enriched
+//! inside known archaic tracts, no amount of HMM tuning can help, and the fault is upstream in the
+//! variant calls or the outgroup strip rather than in the model. That question is unanswerable from
+//! the caller's own output, which is why this exists.
+//!
+//! ```sh
+//! cargo run --release -p navigator-analysis --example archaic_private_dump -- \
+//!   calls.json ~/.decodingus/ancestry/archaic_outgroup_af_chm13v2.0.bin > private.tsv
+//! ```
+
+use navigator_analysis::archaic::ArchaicOutgroup;
+use navigator_analysis::caller::SiteGenotype;
+
+fn main() -> Result<(), Box<dyn std::error::Error>> {
+    let mut a = std::env::args().skip(1);
+    let calls_path = a.next().expect("usage: archaic_private_dump <calls.json> <outgroup.bin>");
+    let og_path = a.next().expect("usage: archaic_private_dump <calls.json> <outgroup.bin>");
+
+    let calls: Vec<SiteGenotype> = serde_json::from_str(&std::fs::read_to_string(&calls_path)?)?;
+    let og = ArchaicOutgroup::from_bytes(&std::fs::read(&og_path)?).map_err(|e| e.to_string())?;
+
+    // Group by contig, mirroring what the caller does before it strips.
+    let mut by_contig: std::collections::BTreeMap<String, Vec<&SiteGenotype>> = Default::default();
+    for c in &calls {
+        by_contig.entry(c.contig.clone()).or_default().push(c);
+    }
+
+    // Quality columns come out too: whether the background's excess variance is real biology or
+    // this caller's own error rate varying by region is not answerable without them.
+    println!("contig\tposition\tdosage\tgq\tdepth");
+    for (contig, mut sites) in by_contig {
+        sites.sort_by_key(|s| s.position);
+        let carried: Vec<&SiteGenotype> = sites.iter().copied().filter(|s| s.dosage > 0).collect();
+        let positions: Vec<i64> = carried.iter().map(|s| s.position).collect();
+        let keep: std::collections::HashSet<i64> = og.retain_private(&contig, &positions).into_iter().collect();
+        let mut kept = 0usize;
+        for s in &carried {
+            if keep.contains(&s.position) {
+                println!("{contig}\t{}\t{}\t{}\t{}", s.position, s.dosage, s.gq, s.depth);
+                kept += 1;
+            }
+        }
+        eprintln!("{contig}: {kept} private of {} carried variants", carried.len());
+    }
+    Ok(())
+}
diff --git a/documents/design/ArchaicAncestry_Design.md b/documents/design/ArchaicAncestry_Design.md
index 3174393..851c248 100644
--- a/documents/design/ArchaicAncestry_Design.md
+++ b/documents/design/ArchaicAncestry_Design.md
@@ -1,11 +1,15 @@
 # Archaic Ancestry Report (Neanderthal / Denisovan) — Design
 
 **Status:** **Tier A SHIPPED** (`v0.1.0-alpha.14`). **Tier B GATED OFF** (`ARCHAIC_SEGMENTS_ENABLED
-= false`) after per-individual validation showed the segment caller carries no per-person signal —
-see *Tier B validation* at the end of §10. It shipped enabled in alpha.14 and was withdrawn in the
-next release. Drafted 2026-07-23; plan added 2026-07-26; all three §9 questions resolved.
-**Read *Tier B validation* and *Deviations from the plan* (end of §10) before trusting anything
-below about Tier B, §7's expected percentage, or M3's feature-gate rule.**
+= false`) after per-individual validation showed the segment caller carries no per-person signal, and
+diagnosed as **built on the wrong observable** — §3's choice of a method designed for people who do
+*not* have archaic reference genomes, which we do. It shipped enabled in alpha.14 and was withdrawn
+in the next release. Drafted 2026-07-23; plan added 2026-07-26; all three §9 questions resolved.
+
+> **Read *Tier B validation* and *Why it failed* (both at the end of §10) before anything else in
+> this document about Tier B.** §3's method choice, §5's Tier B pipeline and M3's calibration are all
+> superseded by that diagnosis. *Deviations from the plan* also qualifies §7's expected percentage
+> and M3's feature-gate rule. Tier A (§5 Tier A, M1, M2) stands unaffected.
 **Goal:** Reconstruct a 23andMe-style Neanderthal report — and go beyond it with a Denisovan
 estimate and a true whole-genome introgression map — from public archaic reference genomes and
 recent methods, using the app's existing ancestry/panel/HMM machinery.
@@ -95,6 +99,15 @@ are population-scale; **IBDmix** is reference-based and heavier; **DAIseg** (joi
 Nea/Den, no post-processing) is attractive but an unreviewed 2025 preprint with self-benchmarks — a
 possible future upgrade, not the v1 foundation.
 
+> **THIS CHOICE WAS WRONG, and the one-clause dismissal of IBDmix is the root of it.** hmmix's
+> premise is detecting introgression *without* archaic reference genomes — it infers them from
+> private-mutation density because it assumes you do not have them. We do: all four archaic genomes,
+> and `archaic_classify` with 2,031,406 diagnostic sites derived from them. Measured, the density
+> observable leaves ~1 informative variant per median tract (14.3 % sensitivity at 5 % false
+> positives) where reference-based allele matching leaves ~30 (95.1 %). "Reference-based and
+> heavier" traded the only thing that made the problem tractable for a saving that was never needed.
+> See *Why it failed* at the end of §10.
+
 ---
 
 ## 3a. Marker-list sourcing — RESOLVED: compute our own
@@ -1038,7 +1051,108 @@ intersection with real 23andMe v5 chip content). It remains enabled and is what
 
 ### Re-enabling
 
-Needs a method change, not a threshold sweep — Skov 2020 matches a segment's whole haplotype against
-each archaic genome relative to a background expectation, where ours tests private-variant density
-against pre-classified sites. Re-running the validation harness that produced the numbers above is
-the gate.
+Needs a method change, not a threshold sweep. **The change required is the OBSERVABLE, not the
+model** — see *Why it failed* immediately below. Re-running the validation harness that produced the
+numbers above is the gate.
+
+---
+
+## Why it failed (2026-07-31) — we chose a method for a problem we do not have
+
+The caller was not mis-tuned. It was built on the wrong observable, and no amount of fitting could
+have rescued it. This section records how that was established, because the wrong turn is easy to
+repeat and every intermediate hypothesis here was plausible.
+
+### The chain of diagnosis
+
+Each step was measured on HG00096, chr21+22, against hmmix's own calls for the same individual.
+
+1. **The signal is real and the pipeline upstream is sound.** Private-variant density inside real
+   archaic tracts is **2.89x** background (357.8 vs 123.9/Mb). The outgroup strip and the variant
+   calls both work.
+2. **The hard-coded `archaic_rate_multiple = 6.0` is roughly double the real effect.** At 1 kb
+   windows and 0.124 background variants/window, the Poisson crossover lands at k* ~ 0.35 — *any*
+   window containing a single private variant is called archaic.
+3. **Per-individual Baum-Welch EM (what hmmix does and we do not) makes it worse.** Unconstrained,
+   it learns a 22.4x multiple and a 0.116 switch probability (~9 kb tracts), calling 17.4 Mb in
+   1,452 segments — 7x the truth. It fits the zero-inflation of sparse counts, not spatial structure.
+4. **Oracle parameters do not help either.** Given the *measured* rates and several tract-length
+   priors, precision stays pinned at 4-5% and every configuration sits at or near its
+   random-placement null. **The model class is wrong, not its parameters.**
+5. **The background is not the distribution the emission assumes.** In callable, non-archaic 100 kb
+   bins: p10-p90 spread **5.3x**, **overdispersion 14.6x** against the Poisson's assumed 1.0x, and
+   **5.9%** of pure-background bins already exceed 2.89x the median. The noise is bigger than the
+   signal, so a single-lambda model calls its own upper tail archaic — and background bins outnumber
+   archaic ones ~20:1.
+6. **A mutation-rate map would not have been enough.** African-outgroup site density — the best
+   proxy available without new data — explains only **38%** of the background variance and takes
+   overdispersion 14.6x -> 7.4x, still above the signal. Quality filtering does not help either: it
+   lowers overdispersion only by discarding variants proportionally, and enrichment falls with it.
+7. **It is not our variant calling.** HG00096 is in the 1000G callset, so the same quantity is
+   computable from their calls. Ours are noisier — 2.12x the carried SNVs on chr21 (93,643 vs
+   44,255), 6x the private-variant density, 97% of our private calls unique to us — but the
+   **contrast is identical**: 1.98x for their calls against 2.08x for ours. Noise dilutes counts,
+   not the ratio.
+8. **It is not the truth set or the lift.** Checked in native hg38 with no lifting, using only
+   hmmix's own segments and their own DAV SNP list: their tracts are enriched **1.84x** for their
+   own archaic SNPs against a null of 1.04x (p95 1.39x). The callset is internally consistent and
+   the lift preserved it.
+
+### The actual cause: one bit of evidence per tract
+
+At ~20-120 private variants/Mb, the **median 36 kb tract contains about ONE informative variant**
+(1.3 expected inside a tract against 0.64 in background). That is the entire evidence separating a
+real tract from noise. Detectability, computed directly:
+
+| observable | evidence per 36 kb tract | sensitivity at <=5% FP | at <=1% FP |
+|---|---|---|---|
+| private-variant density (built) | ~1 variant | **14.3 %** | 4.3 % |
+| archaic-allele matching at diagnostic sites | ~30 sites | **95.1 %** | 80.9 % |
+
+To reach 80% sensitivity at 5% false positives, allele matching needs **36 kb** — the real median
+tract. Density does not get there **at 500 kb**, and hmmix's p10 tract is 7 kb. Both observables
+carry the *same* ~3x contrast (2.89x vs 3.04x), so contrast was never the problem: **evidence per
+tract** was.
+
+### The wrong turn, named
+
+§3 chose **Skov 2018 (hmmix)**, whose premise is detecting introgression **without archaic reference
+genomes** — it infers them indirectly from private-mutation density precisely because it assumes you
+do not have them. We adopted it while **holding all four archaic genomes** and already shipping
+`archaic_classify` with 2,031,406 diagnostic sites derived from them. §3 rejected the reference-based
+alternative in one clause — *"IBDmix is reference-based and heavier"* — and that clause is the root
+of every failure above.
+
+An introgressed tract is **a haplotype inherited intact from an archaic ancestor**. The question is
+whether a stretch of the genome *matches an archaic genome*, not whether it is slightly more mutated
+than average. Runs, not counts.
+
+Measured, the reference-based observable gives **39.5%** archaic-allele carrying inside real tracts
+against **13.0%** elsewhere, over ~30 sites per tract instead of ~1 variant.
+
+### What the reframe also fixes
+
+- **The mutation-rate map becomes unnecessary.** Diagnostic sites are the denominator, so their
+  uneven density cancels out — the thing the map existed to correct.
+- **The 5.3x background spread and 14.6x overdispersion stop applying**, because nothing is being
+  modelled as a density any more.
+- **The assets already exist.** `archaic_classify_<build>.bin` carries the sites, derived bases and
+  per-lineage classes; no new offline pipeline is needed to start.
+- **Attribution may come back with it.** The same per-site matching, split by lineage class, *is*
+  Skov's post-hoc annotation — the thing `attribute_lineage` was gated for.
+
+### Open before this becomes a plan
+
+- **The background carrying rate measures 13.0%**, against the 4.3% recorded earlier in this
+  document. Unreconciled; probably our call noise. It sets the contrast, so it matters.
+- **~30 sites per tract assumes adequate coverage at diagnostic sites.** Check against real call
+  rates rather than assuming.
+- **n = 1 for the reframe.** The 3.04x enrichment is HG00096 only. The harness in
+  `scripts/archaic-validation/` runs the cohort cheaply now that the CRAM fix has landed.
+
+### Instruments
+
+The examples these numbers came from, all under `crates/navigator-analysis/examples/`:
+`archaic_private_dump` (the HMM's actual input, with quality columns), `archaic_outgroup_density`
+(the rate-map proxy), `archaic_classify_dump` (diagnostic sites), `archaic_callable_dump` (what the
+caller can see at all), and `cram_query_probe` (the CRAM defect found on the way here).
diff --git a/scripts/archaic-validation/observable/README.md b/scripts/archaic-validation/observable/README.md
new file mode 100644
index 0000000..c1e635d
--- /dev/null
+++ b/scripts/archaic-validation/observable/README.md
@@ -0,0 +1,52 @@
+# Why Tier B failed — the diagnostic scripts
+
+These are the instruments behind
+[`ArchaicAncestry_Design.md`](../../../documents/design/ArchaicAncestry_Design.md) § *Why it failed*.
+They establish that the Tier B segment caller was built on the **wrong observable**, not mis-tuned —
+and, just as importantly, they rule out the explanations that looked right first.
+
+Kept because every one of the rejected hypotheses is plausible enough to be proposed again.
+
+## Run order
+
+All expect a scratch directory containing, for one individual (HG00096 in the recorded run):
+
+| file | produced by |
+|---|---|
+| `HG00096.chr{21,22}.calls.json` | `sqlite3` on `analysis_artifact`, kind `diploid_denovo:chrN` |
+| `private.chr{21,22}.tsv` | `cargo run --example archaic_private_dump` |
+| `classify.tsv` | `cargo run --example archaic_classify_dump` |
+| `callable.bed` | `cargo run --example archaic_callable_dump` (threshold 0.5) |
+| `og_density.tsv` | `cargo run --example archaic_outgroup_density` |
+| `truth_<sample>.chm13.bed` | hmmix segments lifted hg38→CHM13 — see [../README.md](../README.md) |
+
+## What each one answers
+
+| script | question | answer it gave |
+|---|---|---|
+| `background_variation.py` | is the background the flat Poisson the emission assumes? | **No** — 5.3× p10–p90 spread, **14.6× overdispersed**, larger than the 2.89× signal |
+| `mutrate_check.py` | would a mutation-rate map fix that? | **Not enough** — the best available proxy explains 38% of the variance, leaving 7.4× |
+| `quality_effect.py` | is the excess variance our caller's artifacts? | **No** — filtering lowers overdispersion only by discarding variants; enrichment falls with it |
+| `callset_compare.py` | is our variant calling diluting the signal? | **No** — 1000G's own calls for the same person give the *same* contrast (1.98× vs 2.08×), though ours are 6× noisier |
+| `truth_selfcheck.py` | is the truth set (or my lift) wrong? | **No** — hmmix's tracts are enriched 1.84× for their own archaic SNPs in native hg38, null 1.04× |
+| `observable_compare.py` | is a more specific observable available? | private ∩ diagnostic is unusable (7 sites); the African strip removes almost all diagnostic sites |
+| `haplotype_match.py` | does archaic-allele matching separate tracts? | **Yes** — 39.5% carrying inside tracts vs 13.0% elsewhere, over ~30 sites per tract |
+| `discriminability.py` | how detectable is ONE tract under each observable? | density **14.3%** vs matching **95.1%** sensitivity at 5% false positives |
+
+`discriminability.py` is the one to read first. It needs no inputs and is the whole argument:
+both observables carry ~3× contrast, so contrast was never the problem — **evidence per tract** was.
+
+## The trap these scripts exist to prevent
+
+Two harness bugs in here produced confident wrong answers before being caught, both worth knowing:
+
+- **`haplotype_match.py` originally conditioned on the subject having a variant call** at a
+  diagnostic site, which samples only sites where he already has a variant and reported an
+  impossible ~80% carrying rate against a known 4.3% background. Every diagnostic site in callable
+  territory belongs in the denominator; no call means hom-reference.
+- **Binning at 100 kb to measure in-tract contrast** dilutes a 36 kb tract with 64 kb of background
+  and reported 1.14× where the correct measurement is 1.98×. Measure contrast at tract boundaries,
+  not in fixed bins.
+
+More generally: sensitivity alone is gameable by calling more sequence, because the
+random-placement null rises with it. Always report the null at the extent actually called.
diff --git a/scripts/archaic-validation/observable/background_variation.py b/scripts/archaic-validation/observable/background_variation.py
new file mode 100644
index 0000000..214d26c
--- /dev/null
+++ b/scripts/archaic-validation/observable/background_variation.py
@@ -0,0 +1,106 @@
+"""Is the private-variant background uniform, as the caller's emission model assumes?
+
+The model says: background windows emit Poisson(lambda) with ONE genome-wide lambda, so a window
+with several private variants is evidence of an archaic tract. That inference only holds if the
+non-archaic background is actually flat. hmmix does not assume this — it requires a mutation-rate
+map and scales the expected density per window by it. We have no such asset.
+
+This measures the spread of background density directly, in regions that are callable and NOT in
+hmmix's archaic tracts, and asks how a Poisson model would fare against it.
+"""
+
+import bisect
+import collections
+import math
+
+BIN = 100_000
+CONTIGS = ('chr21', 'chr22')
+
+
+def merge(iv, tol=0):
+    iv = sorted(iv)
+    out, cs, ce = [], *iv[0]
+    for s, e in iv[1:]:
+        if s > ce + tol:
+            out.append((cs, ce))
+            cs, ce = s, e
+        else:
+            ce = max(ce, e)
+    out.append((cs, ce))
+    return out
+
+
+def load_bed(path, tol=0):
+    d = collections.defaultdict(list)
+    for line in open(path):
+        c, s, e = line.split()[:3]
+        d[c].append((int(s), int(e)))
+    return {c: merge(v, tol) for c, v in d.items()}
+
+
+def main():
+    priv = collections.defaultdict(list)
+    for c in CONTIGS:
+        with open(f'private.{c}.tsv') as f:
+            next(f)
+            for line in f:
+                p = line.split('\t')
+                priv[p[0]].append(int(p[1]))
+    for c in priv:
+        priv[c].sort()
+
+    callable_r = load_bed('callable.bed')
+    truth = load_bed('truth_HG00096.chm13.bed', 1000)
+
+    def covered(regions, s, e):
+        tot = 0
+        for s2, e2 in regions:
+            lo, hi = max(s, s2), min(e, e2)
+            if hi > lo:
+                tot += hi - lo
+        return tot
+
+    dens = []
+    for c in CONTIGS:
+        a = priv[c]
+        lo = min(s for s, _ in callable_r[c])
+        hi = max(e for _, e in callable_r[c])
+        for b in range(lo, hi, BIN):
+            e = b + BIN
+            cal = covered(callable_r[c], b, e)
+            if cal < BIN * 0.5:
+                continue                        # mostly uncallable: not a background sample
+            if covered(truth.get(c, []), b, e) > 0:
+                continue                        # overlaps a real archaic tract: not background
+            n = bisect.bisect_left(a, e) - bisect.bisect_left(a, b)
+            dens.append(n / (cal / 1e6))
+    dens.sort()
+    n = len(dens)
+    mean = sum(dens) / n
+    var = sum((d - mean) ** 2 for d in dens) / (n - 1)
+    q = lambda f: dens[int(f * (n - 1))]
+
+    print(f'BACKGROUND private-variant density, {BIN // 1000} kb bins, callable and NOT archaic')
+    print(f'  n bins {n}')
+    print(f'  mean {mean:.0f}/Mb   median {q(.5):.0f}   p10 {q(.1):.0f}   p90 {q(.9):.0f}'
+          f'   max {dens[-1]:.0f}')
+    print(f'  p90/p10 spread = {q(.9) / max(q(.1), 1e-9):.1f}x')
+    print()
+    print(f'  enrichment inside real archaic tracts, for comparison: 2.89x')
+    print(f'  fraction of BACKGROUND bins already above 2.89x the median: '
+          f'{sum(1 for d in dens if d > 2.89 * q(.5)) / n * 100:.1f}%')
+    print()
+    # Poisson would predict variance == mean for the per-bin COUNT; overdispersion is the degree to
+    # which a single-lambda model is simply the wrong distribution.
+    counts = [d * BIN / 1e6 for d in dens]
+    cm = sum(counts) / len(counts)
+    cv = sum((x - cm) ** 2 for x in counts) / (len(counts) - 1)
+    print(f'  per-bin counts: mean {cm:.1f}, variance {cv:.1f} -> overdispersion {cv / cm:.1f}x')
+    print('  (Poisson assumes variance = mean, i.e. 1.0x. Anything well above that means a')
+    print('   single-lambda background will call its own upper tail archaic.)')
+    print()
+    print(f'  sd/mean of background density = {math.sqrt(var) / mean:.2f}')
+
+
+if __name__ == '__main__':
+    main()
diff --git a/scripts/archaic-validation/observable/callset_compare.py b/scripts/archaic-validation/observable/callset_compare.py
new file mode 100644
index 0000000..5c05da7
--- /dev/null
+++ b/scripts/archaic-validation/observable/callset_compare.py
@@ -0,0 +1,100 @@
+"""Same individual, same outgroup, same build: 1000G's call set versus ours.
+
+The Tier B observable is "variants this person carries that no African carries". Its contrast
+inside real archaic tracts is 2.89x for our calls -- too weak to separate from a background that
+varies 5.3x. hmmix's documented example emission rates imply a background near 40/Mb and a ~10x
+contrast, and ours is 124/Mb at 2.89x.
+
+Our outgroup track was just verified complete against its source, so if the observable is diluted
+the dilution must come from the variant calls. HG00096 is in the 1000G callset, so the same
+quantity can be computed from their calls and ours and compared directly.
+"""
+
+import bisect
+import collections
+
+
+def merge(iv, tol=0):
+    iv = sorted(iv)
+    out, cs, ce = [], *iv[0]
+    for s, e in iv[1:]:
+        if s > ce + tol:
+            out.append((cs, ce))
+            cs, ce = s, e
+        else:
+            ce = max(ce, e)
+    out.append((cs, ce))
+    return out
+
+
+def load_bed(path, tol=0):
+    d = collections.defaultdict(list)
+    for line in open(path):
+        c, s, e = line.split()[:3]
+        d[c].append((int(s), int(e)))
+    return {c: merge(v, tol) for c, v in d.items()}
+
+
+def covered(regions, s, e):
+    tot = 0
+    for s2, e2 in regions:
+        lo, hi = max(s, s2), min(e, e2)
+        if hi > lo:
+            tot += hi - lo
+    return tot
+
+
+def stats(positions, callable_r, truth, label):
+    a = sorted(positions)
+    cnt = lambda s, e: bisect.bisect_left(a, e) - bisect.bisect_left(a, s)
+    in_bp = sum(covered(callable_r, s, e) for s, e in truth)
+    in_n = sum(cnt(s, e) for s, e in truth)
+    cb = sum(e - s for s, e in callable_r)
+    cn = sum(cnt(s, e) for s, e in callable_r)
+    out_n, out_bp = cn - in_n, cb - in_bp
+    enrich = (in_n / in_bp) / (out_n / out_bp)
+    print(f'{label:26s} {len(a):7d} private   in-tract {in_n / (in_bp / 1e6):6.1f}/Mb   '
+          f'background {out_n / (out_bp / 1e6):6.1f}/Mb   CONTRAST {enrich:5.2f}x')
+    return enrich
+
+
+def main():
+    afr = set(int(x) for x in open('afr_sites.chr21.txt'))
+    kgp_carried = [int(x) for x in open('kgp_carried.chr21.txt')]
+    ours_carried, ours_private = [], []
+    with open('private.chr21.tsv') as f:
+        next(f)
+        for line in f:
+            p = line.split('\t')
+            ours_private.append(int(p[1]))
+
+    callable_r = load_bed('callable.bed')['chr21']
+    truth = load_bed('truth_HG00096.chm13.bed', 1000)['chr21']
+
+    kgp_private = [p for p in kgp_carried if p not in afr]
+
+    print(f'HG00096 chr21, carried SNVs:  1000G {len(kgp_carried):7d}   ours {93643:7d}'
+          f'   ratio {93643 / len(kgp_carried):.2f}x')
+    print(f'African outgroup: {len(afr)} segregating sites (652 unrelated individuals)\n')
+
+    print('observable                    count      density in tract / background      contrast')
+    e_kgp = stats(kgp_private, callable_r, truth, "1000G's calls")
+    e_ours = stats(ours_private, callable_r, truth, 'our calls')
+    print()
+    print(f'  contrast ratio 1000G/ours = {e_kgp / e_ours:.2f}x')
+    print()
+    print('  If 1000G\'s calls give a much higher contrast on the SAME person, the observable is')
+    print('  diluted by our variant calling, not by the outgroup or the model -- and the fix is')
+    print('  upstream of the HMM entirely.')
+
+    # Where do the extra calls sit? If ours are mostly at positions 1000G did not call at all,
+    # they are ours to explain.
+    kgp_set = set(kgp_carried)
+    extra = [p for p in ours_private if p not in kgp_set]
+    shared = len(ours_private) - len(extra)
+    print(f'\n  of our {len(ours_private)} private calls: {shared} also called by 1000G, '
+          f'{len(extra)} ours alone ({len(extra) / max(len(ours_private), 1) * 100:.0f}%)')
+
+
+if __name__ == '__main__':
+    main()
diff --git a/scripts/archaic-validation/observable/discriminability.py b/scripts/archaic-validation/observable/discriminability.py
new file mode 100644
index 0000000..b548439
--- /dev/null
+++ b/scripts/archaic-validation/observable/discriminability.py
@@ -0,0 +1,73 @@
+"""How detectable is a single tract under each observable?
+
+Contrast alone does not decide this. What decides it is how much evidence one tract carries, and the
+two observables differ 30-fold there:
+
+  private-variant density : ~1 informative variant per 36 kb tract (Poisson 0.6 vs 1.3)
+  archaic-allele matching : ~30 diagnostic sites  (Binomial 30 x 13% vs 30 x 39.5%)
+
+Same ~3x contrast, wildly different decidability. This computes, for each, the sensitivity
+achievable at a fixed false-positive rate -- i.e. what fraction of real tracts a perfect classifier
+could call while keeping background calls rare.
+"""
+
+import math
+
+
+def binom_pmf(k, n, p):
+    return math.comb(n, k) * p ** k * (1 - p) ** (n - k)
+
+
+def pois_pmf(k, lam):
+    return math.exp(-lam) * lam ** k / math.factorial(k)
+
+
+def roc(bg, sig, kmax):
+    """Sensitivity at the threshold where background false-positive rate first drops below 5%/1%."""
+    out = {}
+    for target in (0.05, 0.01):
+        for t in range(kmax + 1):
+            fp = sum(bg(k) for k in range(t, kmax + 1))
+            if fp <= target:
+                out[target] = (t, sum(sig(k) for k in range(t, kmax + 1)), fp)
+                break
+        else:
+            out[target] = (None, 0.0, 0.0)
+    return out
+
+
+print('DETECTABILITY OF ONE TRACT (36 kb, the measured median)\n')
+
+# --- private-variant density ---------------------------------------------------------------
+lam_bg, lam_sig = 0.64, 1.30
+r = roc(lambda k: pois_pmf(k, lam_bg), lambda k: pois_pmf(k, lam_sig), 30)
+print(f'private-variant density   background Poisson({lam_bg})  tract Poisson({lam_sig})')
+for tgt, (t, sens, fp) in r.items():
+    print(f'   at <= {tgt:.0%} false positives: threshold k >= {t}, sensitivity {sens:6.1%}')
+
+# --- archaic-allele matching ---------------------------------------------------------------
+n, p_bg, p_sig = 30, 0.130, 0.395
+r = roc(lambda k: binom_pmf(k, n, p_bg), lambda k: binom_pmf(k, n, p_sig), n)
+print(f'\narchaic-allele matching   background Binom({n}, {p_bg})  tract Binom({n}, {p_sig})')
+for tgt, (t, sens, fp) in r.items():
+    print(f'   at <= {tgt:.0%} false positives: threshold k >= {t}, sensitivity {sens:6.1%}')
+
+# What tract size does each need to become usable?
+print('\nTRACT SIZE NEEDED FOR 80% SENSITIVITY AT 5% FALSE POSITIVES')
+for label, kind in (('private density', 'pois'), ('allele matching', 'binom')):
+    for kb in (10, 20, 36, 50, 100, 200, 500):
+        scale = kb / 36
+        if kind == 'pois':
+            bg, sig = lam_bg * scale, lam_sig * scale
+            rr = roc(lambda k: pois_pmf(k, bg), lambda k: pois_pmf(k, sig), 60)
+        else:
+            nn = max(int(round(n * scale)), 1)
+            rr = roc(lambda k: binom_pmf(k, nn, p_bg), lambda k: binom_pmf(k, nn, p_sig), nn)
+        if rr[0.05][1] >= 0.80:
+            print(f'  {label:18s} {kb:4d} kb  -> sensitivity {rr[0.05][1]:.0%}')
+            break
+    else:
+        print(f'  {label:18s} not reached by 500 kb')
+
+print('\n  hmmix p10 tract is 7 kb and the median 31-36 kb, so a method needing hundreds of kb')
+print('  cannot report tracts at the resolution the feature claims.')
diff --git a/scripts/archaic-validation/observable/haplotype_match.py b/scripts/archaic-validation/observable/haplotype_match.py
new file mode 100644
index 0000000..e6033aa
--- /dev/null
+++ b/scripts/archaic-validation/observable/haplotype_match.py
@@ -0,0 +1,121 @@
+"""Should the observable be archaic-ALLELE MATCHING rather than private-mutation density?
+
+The current caller counts private variants per window: measured contrast inside real archaic tracts
+is ~2x, and at ~20-120 private variants/Mb a median 31 kb tract carries about ONE informative
+variant. That is why no parameter setting works.
+
+But an introgressed tract is a haplotype inherited intact from an archaic ancestor, and we HOLD the
+archaic genomes -- 2,031,406 diagnostic sites where archaics carry a derived allele. A real tract
+should carry the archaic allele at a large fraction of the diagnostic sites it spans, while a
+non-introgressed region carries it only at the background rate (~4%, per the design's own
+measurement).
+
+Diagnostic sites are the denominator here, not megabases, which controls for their uneven density
+for free -- the thing a mutation-rate map was going to have to correct.
+"""
+
+import bisect
+import collections
+import json
+import random
+
+CONTIGS = ('chr21', 'chr22')
+
+
+def merge(iv, tol=0):
+    iv = sorted(iv)
+    out, cs, ce = [], *iv[0]
+    for s, e in iv[1:]:
+        if s > ce + tol:
+            out.append((cs, ce))
+            cs, ce = s, e
+        else:
+            ce = max(ce, e)
+    out.append((cs, ce))
+    return out
+
+
+def load_bed(path, tol=0):
+    d = collections.defaultdict(list)
+    for line in open(path):
+        c, s, e = line.split()[:3]
+        d[c].append((int(s), int(e)))
+    return {c: merge(v, tol) for c, v in d.items()}
+
+
+def main():
+    # Diagnostic sites: position -> (derived base, class)
+    diag = collections.defaultdict(dict)
+    with open('classify.tsv') as f:
+        next(f)
+        for line in f:
+            c, p, d, k = line.rstrip('\n').split('\t')
+            diag[c][int(p)] = (d, int(k))
+
+    # The subject's calls, so we can ask whether he carries the archaic allele.
+    carries = collections.defaultdict(dict)
+    for c in CONTIGS:
+        for rec in json.load(open(f'HG00096.{c}.calls.json')):
+            ref = rec['reference_allele']
+            alt = rec['alternate_allele']
+            dos = rec['dosage']
+            if dos is None or dos < 0:
+                continue
+            alleles = set()
+            if dos < 2:
+                alleles.add(ref)
+            if dos > 0:
+                alleles.add(alt)
+            carries[rec['contig']][rec['position']] = alleles
+
+    callable_r = load_bed('callable.bed')
+    truth = load_bed('truth_HG00096.chm13.bed', 1000)
+
+    def in_regions(regions, p):
+        i = bisect.bisect_right([s for s, _ in regions], p) - 1
+        return i >= 0 and regions[i][1] > p
+
+    stats = {'in': [0, 0], 'out': [0, 0]}
+    for c in CONTIGS:
+        cal = callable_r[c]
+        tr = truth.get(c, [])
+        for p, (derived, _k) in diag[c].items():
+            if not in_regions(cal, p):
+                continue
+            # EVERY diagnostic site in callable territory is in the denominator. The caller emits
+            # only variant records, so no record means hom-reference — i.e. NOT carrying the
+            # archaic allele, given the panel orients derived as ALT. Conditioning on "has a call"
+            # instead samples only sites where he already has a variant, which is why that version
+            # reported an impossible ~80% carrying rate against a known 4.3% background.
+            called = carries[c].get(p)
+            bucket = 'in' if in_regions(tr, p) else 'out'
+            stats[bucket][1] += 1
+            if called is not None and derived in called:
+                stats[bucket][0] += 1
+
+    print('ARCHAIC-ALLELE CARRYING RATE at diagnostic sites the subject has a call for')
+    for b, label in (('in', "inside hmmix's tracts"), ('out', 'elsewhere')):
+        hit, tot = stats[b]
+        rate = hit / tot * 100 if tot else 0
+        print(f'  {label:26s} {hit:6d} / {tot:6d} = {rate:5.1f}%')
+    ri = stats['in'][0] / max(stats['in'][1], 1)
+    ro = stats['out'][0] / max(stats['out'][1], 1)
+    print(f'  ENRICHMENT {ri / ro if ro else 0:.2f}x')
+    print()
+    print('  Compare: private-variant DENSITY gives 2.89x with ~1 informative variant per tract.')
+    print('  Sites here are the denominator, so uneven diagnostic-site density cancels out.')
+
+    # How much evidence does a typical tract actually carry under this observable?
+    spans = [e - s for c in CONTIGS for s, e in truth.get(c, [])]
+    spans.sort()
+    med = spans[len(spans) // 2]
+    dens = sum(len(diag[c]) for c in CONTIGS) / sum(
+        e - s for c in CONTIGS for s, e in merge([(min(x for x, _ in callable_r[c]),
+                                                   max(y for _, y in callable_r[c]))]))
+    print(f'\n  median tract {med / 1000:.0f} kb; diagnostic sites ~{dens * 1e6:.0f}/Mb'
+          f' -> ~{med * dens:.0f} informative sites per tract')
+    print('  (the density model gets ~1 private variant for the same tract)')
+
+
+if __name__ == '__main__':
+    main()
diff --git a/scripts/archaic-validation/observable/mutrate_check.py b/scripts/archaic-validation/observable/mutrate_check.py
new file mode 100644
index 0000000..05e6105
--- /dev/null
+++ b/scripts/archaic-validation/observable/mutrate_check.py
@@ -0,0 +1,132 @@
+"""Does African-outgroup site density explain the background variation the emission model ignores?
+
+If it does, it is a usable mutation-rate normalizer and the fix is an asset build. If it does not,
+a different proxy is needed and building this one would waste the effort.
+
+Two questions, in order:
+  1. Correlation between outgroup density and private-variant density in BACKGROUND regions
+     (callable, non-archaic). This is the normalizer's whole job.
+  2. Whether normalizing by it actually flattens the background — the overdispersion should fall
+     from the measured 14.6x toward 1.0x. Correlation alone does not guarantee that.
+"""
+
+import bisect
+import collections
+import math
+
+BIN = 100_000
+CONTIGS = ('chr21', 'chr22')
+
+
+def merge(iv, tol=0):
+    iv = sorted(iv)
+    out, cs, ce = [], *iv[0]
+    for s, e in iv[1:]:
+        if s > ce + tol:
+            out.append((cs, ce))
+            cs, ce = s, e
+        else:
+            ce = max(ce, e)
+    out.append((cs, ce))
+    return out
+
+
+def load_bed(path, tol=0):
+    d = collections.defaultdict(list)
+    for line in open(path):
+        c, s, e = line.split()[:3]
+        d[c].append((int(s), int(e)))
+    return {c: merge(v, tol) for c, v in d.items()}
+
+
+def covered(regions, s, e):
+    tot = 0
+    for s2, e2 in regions:
+        lo, hi = max(s, s2), min(e, e2)
+        if hi > lo:
+            tot += hi - lo
+    return tot
+
+
+def pearson(x, y):
+    n = len(x)
+    mx, my = sum(x) / n, sum(y) / n
+    num = sum((a - mx) * (b - my) for a, b in zip(x, y))
+    dx = math.sqrt(sum((a - mx) ** 2 for a in x))
+    dy = math.sqrt(sum((b - my) ** 2 for b in y))
+    return num / (dx * dy) if dx and dy else 0.0
+
+
+def main():
+    priv = collections.defaultdict(list)
+    for c in CONTIGS:
+        with open(f'private.{c}.tsv') as f:
+            next(f)
+            for line in f:
+                p = line.split('\t')
+                priv[p[0]].append(int(p[1]))
+    for c in priv:
+        priv[c].sort()
+
+    og = collections.defaultdict(collections.Counter)
+    with open('og_density.tsv') as f:
+        next(f)
+        for line in f:
+            c, w, k = line.split('\t')
+            og[c][int(w)] = int(k)
+
+    callable_r = load_bed('callable.bed')
+    truth = load_bed('truth_HG00096.chm13.bed', 1000)
+
+    xs, ys, weights = [], [], []
+    for c in CONTIGS:
+        a = priv[c]
+        lo = min(s for s, _ in callable_r[c])
+        hi = max(e for _, e in callable_r[c])
+        for b in range(lo, hi, BIN):
+            e = b + BIN
+            cal = covered(callable_r[c], b, e)
+            if cal < BIN * 0.5 or covered(truth.get(c, []), b, e) > 0:
+                continue
+            n = bisect.bisect_left(a, e) - bisect.bisect_left(a, b)
+            ogn = sum(og[c].get(w, 0) for w in range(b, e, 1000))
+            if ogn == 0:
+                continue
+            xs.append(ogn / (cal / 1e6))     # outgroup sites per callable Mb
+            ys.append(n / (cal / 1e6))       # our private variants per callable Mb
+            weights.append(cal)
+
+    r = pearson(xs, ys)
+    print(f'BACKGROUND bins ({BIN // 1000} kb, callable, non-archaic): n = {len(xs)}')
+    print(f'  outgroup-site density vs private-variant density:  Pearson r = {r:+.3f}'
+          f'   (r^2 = {r * r:.2f})')
+    print(f'  -> outgroup density explains {r * r * 100:.0f}% of the background variance')
+    print()
+
+    # Does dividing by the proxy actually flatten it?
+    mean_x = sum(xs) / len(xs)
+    raw = ys
+    norm = [y / (x / mean_x) for x, y in zip(xs, ys)]
+
+    def stats(v, label):
+        v = sorted(v)
+        n = len(v)
+        m = sum(v) / n
+        q = lambda f: v[int(f * (n - 1))]
+        # overdispersion of the implied per-bin count
+        cnt = [d * BIN / 1e6 for d in v]
+        cm = sum(cnt) / len(cnt)
+        cv = sum((z - cm) ** 2 for z in cnt) / (len(cnt) - 1)
+        print(f'  {label:22s} p10 {q(.1):6.0f}  median {q(.5):6.0f}  p90 {q(.9):6.0f}'
+              f'   p90/p10 {q(.9) / max(q(.1), 1e-9):5.1f}x   overdispersion {cv / cm:5.1f}x')
+
+    print('BACKGROUND FLATNESS (the emission model assumes this is flat):')
+    stats(raw, 'raw density')
+    stats(norm, 'normalized by proxy')
+    print()
+    print('  The archaic signal is 2.89x. For the model to separate it, the background spread')
+    print('  must be well below that. Overdispersion near 1.0x is what a Poisson emission needs.')
+
+
+if __name__ == '__main__':
+    main()
diff --git a/scripts/archaic-validation/observable/observable_compare.py b/scripts/archaic-validation/observable/observable_compare.py
new file mode 100644
index 0000000..e5f7799
--- /dev/null
+++ b/scripts/archaic-validation/observable/observable_compare.py
@@ -0,0 +1,123 @@
+"""Compare candidate observables for the Tier B HMM on signal-to-noise.
+
+The current observable is "any private variant". Measured, that gives 2.89x enrichment inside real
+archaic tracts against a background that varies 5.3x (p10-p90) and is 14.6x overdispersed — the
+noise is bigger than the signal, which is why no parameter setting works.
+
+A more specific observable should trade count for contrast: fewer observations, but a much higher
+ratio inside tracts. What matters is whether the enrichment clears the background spread, because
+that ratio is what any two-state model has to separate.
+"""
+
+import bisect
+import collections
+
+BIN = 100_000
+CONTIGS = ('chr21', 'chr22')
+
+
+def merge(iv, tol=0):
+    iv = sorted(iv)
+    out, cs, ce = [], *iv[0]
+    for s, e in iv[1:]:
+        if s > ce + tol:
+            out.append((cs, ce))
+            cs, ce = s, e
+        else:
+            ce = max(ce, e)
+    out.append((cs, ce))
+    return out
+
+
+def load_bed(path, tol=0):
+    d = collections.defaultdict(list)
+    for line in open(path):
+        c, s, e = line.split()[:3]
+        d[c].append((int(s), int(e)))
+    return {c: merge(v, tol) for c, v in d.items()}
+
+
+def covered(regions, s, e):
+    tot = 0
+    for s2, e2 in regions:
+        lo, hi = max(s, s2), min(e, e2)
+        if hi > lo:
+            tot += hi - lo
+    return tot
+
+
+def assess(pos_by_contig, callable_r, truth, label):
+    in_n = in_bp = out_n = out_bp = 0
+    dens = []
+    for c in CONTIGS:
+        a = sorted(pos_by_contig.get(c, []))
+        cnt = lambda s, e: bisect.bisect_left(a, e) - bisect.bisect_left(a, s)
+        in_bp += sum(covered(callable_r[c], s, e) for s, e in truth[c])
+        in_n += sum(cnt(s, e) for s, e in truth[c])
+        cb = sum(e - s for s, e in callable_r[c])
+        out_bp += cb - sum(covered(callable_r[c], s, e) for s, e in truth[c])
+        out_n += sum(cnt(s, e) for s, e in callable_r[c]) - sum(cnt(s, e) for s, e in truth[c])
+        lo = min(s for s, _ in callable_r[c])
+        hi = max(e for _, e in callable_r[c])
+        for b in range(lo, hi, BIN):
+            e = b + BIN
+            cal = covered(callable_r[c], b, e)
+            if cal < BIN * 0.5 or covered(truth.get(c, []), b, e) > 0:
+                continue
+            dens.append(cnt(b, e) / (cal / 1e6))
+    if not in_bp or not out_n:
+        print(f'{label:38s} (insufficient data)')
+        return
+    enrich = (in_n / in_bp) / (out_n / out_bp)
+    dens.sort()
+    n = len(dens)
+    q = lambda f: dens[int(f * (n - 1))]
+    spread = q(.9) / q(.1) if q(.1) > 0 else float('inf')
+    print(f'{label:38s} n={in_n + out_n:6d}  in-tract {in_n / (in_bp / 1e6):6.1f}/Mb  '
+          f'bg {out_n / (out_bp / 1e6):6.1f}/Mb  ENRICH {enrich:5.2f}x  bg p90/p10 {spread:6.1f}x')
+
+
+def main():
+    priv = collections.defaultdict(list)
+    for c in CONTIGS:
+        with open(f'private.{c}.tsv') as f:
+            next(f)
+            for line in f:
+                p = line.rstrip('\n').split('\t')
+                priv[p[0]].append((int(p[1]), int(p[2]), int(p[3]), int(p[4])))
+
+    diag = collections.defaultdict(dict)
+    with open('classify.tsv') as f:
+        next(f)
+        for line in f:
+            c, p, d, k = line.rstrip('\n').split('\t')
+            diag[c][int(p)] = (d, int(k))
+
+    callable_r = load_bed('callable.bed')
+    truth = load_bed('truth_HG00096.chm13.bed', 1000)
+
+    print('observable                             count     density in/out         contrast')
+    assess({c: [r[0] for r in v] for c, v in priv.items()}, callable_r, truth,
+           'all private variants (current)')
+
+    # Private variants that also sit at a known archaic-diagnostic site.
+    d_all = {c: [r[0] for r in v if r[0] in diag.get(c, {})] for c, v in priv.items()}
+    assess(d_all, callable_r, truth, 'private AND archaic-diagnostic')
+
+    for k, name in ((0, 'Neanderthal-diagnostic'), (2, 'shared-archaic')):
+        f = {c: [r[0] for r in v if diag.get(c, {}).get(r[0], (None, -1))[1] == k]
+             for c, v in priv.items()}
+        assess(f, callable_r, truth, f'  ...of which {name}')
+
+    # Control: diagnostic sites the subject does NOT carry should show no enrichment.
+    carried = {c: {r[0] for r in v} for c, v in priv.items()}
+    notc = {c: [p for p in diag.get(c, {}) if p not in carried.get(c, set())] for c in CONTIGS}
+    assess(notc, callable_r, truth, 'diagnostic sites NOT carried (control)')
+
+    print()
+    print('  Contrast is what matters: the enrichment has to clear the background spread for a')
+    print('  two-state model to separate the states. "all private" fails that (2.89x vs 5.3x).')
+
+
+if __name__ == '__main__':
+    main()
diff --git a/scripts/archaic-validation/observable/quality_effect.py b/scripts/archaic-validation/observable/quality_effect.py
new file mode 100644
index 0000000..23f008c
--- /dev/null
+++ b/scripts/archaic-validation/observable/quality_effect.py
@@ -0,0 +1,123 @@
+"""Is the background's excess variance real, or is it this caller's error rate varying by region?
+
+The emission model needs a flat background. Measured raw, it is 14.6x overdispersed; normalizing by
+African-outgroup density gets it to 7.4x, still above the 2.89x archaic signal. If the remainder is
+low-confidence calls -- artifacts clustering in hard regions -- then filtering should flatten the
+background AND raise the archaic enrichment, because artifacts dilute real tracts too.
+
+Reports both numbers for each filter, since a filter that flattens the background by discarding the
+signal has bought nothing.
+"""
+
+import bisect
+import collections
+
+BIN = 100_000
+CONTIGS = ('chr21', 'chr22')
+
+
+def merge(iv, tol=0):
+    iv = sorted(iv)
+    out, cs, ce = [], *iv[0]
+    for s, e in iv[1:]:
+        if s > ce + tol:
+            out.append((cs, ce))
+            cs, ce = s, e
+        else:
+            ce = max(ce, e)
+    out.append((cs, ce))
+    return out
+
+
+def load_bed(path, tol=0):
+    d = collections.defaultdict(list)
+    for line in open(path):
+        c, s, e = line.split()[:3]
+        d[c].append((int(s), int(e)))
+    return {c: merge(v, tol) for c, v in d.items()}
+
+
+def covered(regions, s, e):
+    tot = 0
+    for s2, e2 in regions:
+        lo, hi = max(s, s2), min(e, e2)
+        if hi > lo:
+            tot += hi - lo
+    return tot
+
+
+def load_private():
+    rows = collections.defaultdict(list)
+    for c in CONTIGS:
+        with open(f'private.{c}.tsv') as f:
+            next(f)
+            for line in f:
+                p = line.rstrip('\n').split('\t')
+                rows[p[0]].append((int(p[1]), int(p[2]), int(p[3]), int(p[4])))
+    return rows
+
+
+def assess(rows, callable_r, truth, label):
+    pos = {c: sorted(r[0] for r in v) for c, v in rows.items()}
+    in_n = in_bp = out_n = out_bp = 0
+    dens = []
+    for c in CONTIGS:
+        a = pos.get(c, [])
+        cnt = lambda s, e: bisect.bisect_left(a, e) - bisect.bisect_left(a, s)
+        tb = sum(covered(callable_r[c], s, e) for s, e in truth[c])
+        tn = sum(cnt(s, e) for s, e in truth[c])
+        cb = sum(e - s for s, e in callable_r[c])
+        cn = sum(cnt(s, e) for s, e in callable_r[c])
+        in_n += tn
+        in_bp += tb
+        out_n += cn - tn
+        out_bp += cb - tb
+        lo = min(s for s, _ in callable_r[c])
+        hi = max(e for _, e in callable_r[c])
+        for b in range(lo, hi, BIN):
+            e = b + BIN
+            cal = covered(callable_r[c], b, e)
+            if cal < BIN * 0.5 or covered(truth.get(c, []), b, e) > 0:
+                continue
+            dens.append(cnt(b, e) / (cal / 1e6))
+    enrich = (in_n / in_bp) / (out_n / out_bp) if out_n and in_bp else 0
+    dens.sort()
+    n = len(dens)
+    q = lambda f: dens[int(f * (n - 1))]
+    counts = [d * BIN / 1e6 for d in dens]
+    cm = sum(counts) / len(counts)
+    cv = sum((z - cm) ** 2 for z in counts) / (len(counts) - 1)
+    kept = sum(len(v) for v in rows.values())
+    print(f'{label:26s} kept {kept:6d}  enrich {enrich:4.2f}x  '
+          f'bg p90/p10 {q(.9) / max(q(.1), 1e-9):5.1f}x  overdisp {cv / max(cm, 1e-9):5.1f}x')
+
+
+def main():
+    raw = load_private()
+    callable_r = load_bed('callable.bed')
+    truth = load_bed('truth_HG00096.chm13.bed', 1000)
+
+    print('filter                      variants   archaic signal   background flatness')
+    assess(raw, callable_r, truth, 'none (current)')
+    for gq in (20, 30, 50):
+        f = {c: [r for r in v if r[2] >= gq] for c, v in raw.items()}
+        assess(f, callable_r, truth, f'GQ >= {gq}')
+    for dp in (10, 15, 20):
+        f = {c: [r for r in v if r[3] >= dp] for c, v in raw.items()}
+        assess(f, callable_r, truth, f'depth >= {dp}')
+    # Homozygous-derived only: an introgressed tract is usually heterozygous, so this should HURT
+    # if the signal is real -- a useful control that the enrichment is not an artifact of genotype.
+    f = {c: [r for r in v if r[1] == 1] for c, v in raw.items()}
+    assess(f, callable_r, truth, 'het only (control)')
+    f = {c: [r for r in v if r[1] == 2] for c, v in raw.items()}
+    assess(f, callable_r, truth, 'hom-alt only (control)')
+    for gq, dp in ((30, 15), (50, 20)):
+        f = {c: [r for r in v if r[2] >= gq and r[3] >= dp] for c, v in raw.items()}
+        assess(f, callable_r, truth, f'GQ>={gq} & depth>={dp}')
+    print()
+    print('  Wanted: enrichment UP and overdispersion DOWN together. Overdispersion must fall well')
+    print('  below the enrichment for a Poisson emission to separate the two states.')
+
+
+if __name__ == '__main__':
+    main()
diff --git a/scripts/archaic-validation/observable/truth_selfcheck.py b/scripts/archaic-validation/observable/truth_selfcheck.py
new file mode 100644
index 0000000..97be171
--- /dev/null
+++ b/scripts/archaic-validation/observable/truth_selfcheck.py
@@ -0,0 +1,120 @@
+"""Is my lifted truth actually where hmmix says it is?
+
+Everything downstream rests on it: I concluded the segment caller's locations are below chance, and
+a shipped feature was gated off on that basis. But measuring private-variant density inside those
+tracts gives only ~2x enrichment even in 1000G's OWN calls -- and hmmix found these tracts BY that
+density, so it should be far higher. Either their method is weaker than advertised, or my lift put
+the tracts in the wrong place.
+
+This checks the truth in native hg38, with no lifting anywhere, using only hmmix's own two files:
+their segment calls and their DAV (derived-in-archaic) SNP list. If their segments are enriched for
+their own archaic SNPs in hg38, the callset is self-consistent and any failure is in my lift.
+"""
+
+import bisect
+import collections
+import os
+import random
+
+SEG = os.path.expanduser('~/.decodingus/ancestry-build/tmp/hmmix_segments_chr21_22.tsv')
+SNPS = os.path.expanduser('~/.decodingus/ancestry-build/raw/hmmix/hg38_1000g_SNPS.txt')
+SAMPLE = 'HG00096'
+CONTIGS = ('chr21', 'chr22')
+
+
+def merge(iv, tol=0):
+    iv = sorted(iv)
+    out, cs, ce = [], *iv[0]
+    for s, e in iv[1:]:
+        if s > ce + tol:
+            out.append((cs, ce))
+            cs, ce = s, e
+        else:
+            ce = max(ce, e)
+    out.append((cs, ce))
+    return out
+
+
+def main():
+    # hmmix's own segments for this sample, in hg38, unioned across haplotypes.
+    seg = collections.defaultdict(list)
+    seg_all = collections.defaultdict(list)
+    with open(SEG) as f:
+        next(f)
+        for line in f:
+            p = line.rstrip('\n').split('\t')
+            if p[4] not in CONTIGS:
+                continue
+            seg_all[p[4]].append((int(p[5]), int(p[6])))
+            if p[0] == SAMPLE:
+                seg[p[4]].append((int(p[5]), int(p[6])))
+    S = {c: merge(v) for c, v in seg.items()}
+
+    # hmmix's own archaic SNP list, in hg38.
+    dav = collections.defaultdict(list)
+    with open(SNPS) as f:
+        next(f)
+        for line in f:
+            p = line.split('\t')
+            if p[0] in CONTIGS:
+                dav[p[0]].append(int(p[1]))
+    for c in dav:
+        dav[c].sort()
+
+    print(f'hmmix segments for {SAMPLE} (hg38, unioned): '
+          f'{sum(len(v) for v in S.values())} tracts, '
+          f'{sum(e - s for v in S.values() for s, e in v) / 1e6:.3f} Mb')
+    print(f'hmmix DAV archaic SNPs on chr21+22: {sum(len(v) for v in dav.values())}')
+    print()
+
+    # Span over which everyone's tracts fall -- the region the sampling actually covers.
+    bounds = {c: (min(s for s, _ in v), max(e for _, e in v)) for c, v in seg_all.items()}
+
+    tot_in = tot_bp = 0
+    tot_out = tot_out_bp = 0
+    for c in CONTIGS:
+        a = dav[c]
+        cnt = lambda s, e: bisect.bisect_left(a, e) - bisect.bisect_left(a, s)
+        lo, hi = bounds[c]
+        in_n = sum(cnt(s, e) for s, e in S[c])
+        in_bp = sum(e - s for s, e in S[c])
+        all_n = cnt(lo, hi)
+        all_bp = hi - lo
+        tot_in += in_n
+        tot_bp += in_bp
+        tot_out += all_n - in_n
+        tot_out_bp += all_bp - in_bp
+
+    d_in = tot_in / (tot_bp / 1e6)
+    d_out = tot_out / (tot_out_bp / 1e6)
+    print('IN NATIVE hg38, NO LIFTING:')
+    print(f'  archaic SNPs inside {SAMPLE}\'s own hmmix tracts : {d_in:7.1f} /Mb  ({tot_in} in '
+          f'{tot_bp / 1e6:.2f} Mb)')
+    print(f'  archaic SNPs elsewhere                        : {d_out:7.1f} /Mb')
+    print(f'  ENRICHMENT                                    : {d_in / d_out:6.2f}x')
+    print()
+    print('  hmmix called these tracts from archaic-variant density, so a strong enrichment here')
+    print('  means their callset is internally consistent and my LIFT is what to distrust.')
+    print('  A weak enrichment here means the truth was never as sharp as assumed.')
+
+    # Null: same tract lengths placed at random in the same span.
+    rnd = random.Random(0)
+    vals = []
+    for _ in range(200):
+        n = bp = 0
+        for c in CONTIGS:
+            a = dav[c]
+            lo, hi = bounds[c]
+            for s, e in S[c]:
+                L = e - s
+                p = rnd.randint(lo, max(lo, hi - L))
+                n += bisect.bisect_left(a, p + L) - bisect.bisect_left(a, p)
+                bp += L
+        vals.append(n / (bp / 1e6) / d_out)
+    vals.sort()
+    print(f'  null (same tracts placed at random): mean {sum(vals) / len(vals):.2f}x  '
+          f'p95 {vals[189]:.2f}x')
+
+
+if __name__ == '__main__':
+    main()