Refine theoretical justification for alpha estimation.

Document limitations of various $\alpha$-estimation methods, clarifying that tree-based methods require a tree for accurate heterogeneity estimation and that `PartitionDispersion` avoids violating identifiability requirements.

DevDocMD/theory/evolutionary_distances.md

@@ -2467,6 +2467,213 @@ mechanical substitution but weren't independently checked against an
 `ape`-equivalent reference for those three specifically — flagged here, not
 silently assumed correct.
 
+### α-estimation literature survey and design discussion (2026-09-11)
+
+Follow-up discussion after implementing `PartitionDispersion::estimate_alpha`
+(previous section), triggered by the user pointing out that the literature
+conflates two distinct problems: *estimating α from an alignment* versus
+*using a given α in a Jin-Nei-style corrected distance*. Jin & Nei (1990)
+itself is squarely in the second camp — α is a user-supplied parameter
+there, not something their paper estimates (confirmed by reading the paper
+directly, see previous section). This section catalogs the actual
+α-estimation literature and records which approaches were considered and
+why each was or wasn't adopted.
+
+**Bibliography** (methodological papers, not exhaustive users of `+Γ`):
+
+| Method | Reference | Needs a tree? | Estimates "true" α or a task-specific one? | Cost |
+|---|---|---|---|---|
+| ML, continuous Γ (foundational) | **Gu, Fu & Li (1995)**, MBE 12:546–557 — read in full | yes, topology fixed a priori | true α (+ θ jointly, invariant+Γ) | very high, ≤5–6 taxa in practice |
+| ML, discrete Γ | Yang (1993/94), MBE 10:1396; JME 39:306 | yes | true α | high |
+| Parsimony-based ML on inferred changes | Yang & Kumar (1996), MBE 13:650 | yes, equal branch lengths assumed | true α, biased when branches unequal | low |
+| Corrected substitutions/site → NB fit | **Gu & Zhang (1997)**, MBE 14:1106 — read in full | yes, + ancestral reconstruction | true α, close to ML accuracy | low (given a tree) |
+| Topology-optimizing grid search | **Guindon & Gascuel (2002)**, MBE 19:534 — read in full | yes, but built *from* the candidate distances themselves (BIONJ) | **not** true α — deliberately biased *upward* for topological accuracy | low–moderate (grid × tree build) |
+| Free (non-parametric) rate distribution | Susko et al. (2003), Syst. Biol. 52:594 | yes | tests the Γ assumption itself | high |
+| Parsimony counts/site → NB fit | **Wakeley (1993)**, JME 37:613 — read in full | yes, NJ tree + Fitch parsimony ancestral reconstruction | true α, biased *upward* by parsimony (conservative) | low (given a tree) |
+| Poisson-Gamma/negative-binomial, foundational | Uzzell & Corbin (1971), Science 172:1089 — **not read**, only seen via citations | unknown — likely also count-per-site/tree-based, unverified | true α (indirect, via NB overdispersion) | unknown |
+| Bayesian posterior over α | BPP/MrBayes/BEAST2/RevBayes | yes | posterior, not point estimate | very high |
+
+**Correction (2026-09-11, after reading Wakeley in full):** an earlier
+revision of this table filed Wakeley (1993) under "no tree needed" —
+wrong. Wakeley's method: build a tree (neighbor-joining, Saitou & Nei
+1987, on the real mtDNA data; coalescent-simulated for the simulation
+study), reconstruct ancestral states at internal nodes by Fitch (1971)
+parsimony, count the *minimum* number of changes per site this implies,
+then fit either a two-rate Poisson mixture or a gamma-distributed-rates
+model (⟹ negative binomial, parameters by the methods of Cohen 1965 and
+Johnson & Kotz 1969) to the resulting per-site count distribution — the
+same general "tree + ancestral reconstruction + per-site counts + NB fit"
+shape as Gu & Zhang (1997), just using raw parsimony counts directly
+rather than Gu & Zhang's multiple-hit-corrected `k̂`. Confirmed by
+simulation (Wakeley's own Fig. 1) that parsimony underestimates both the
+mean *and* (more severely) the variance of per-site change counts, which
+biases `α` **upward** (toward apparent uniformity) — consistent with Gu &
+Zhang's own citation of Wakeley for exactly this bias, and with the
+general "parsimony overestimates α" pattern noted throughout this
+literature. Real hypervariable-region-1 mtDNA data gave `α̂ ≈ 0.44–0.60`
+per data set, `0.47` combined — despite the conservative bias, still
+solidly in the "strong heterogeneity" range, which is the paper's own
+point: the method is biased but usefully conservative, not useless.
+
+**Where `PartitionDispersion::estimate_alpha` actually sits**: nowhere in
+this table's tree-dependent rows. It shares only the *abstract
+mathematical identity* (Poisson-Gamma mixture ⟹ negative binomial
+over-dispersion) with Wakeley/Uzzell-Corbin/Gu&Zhang, applied via a
+genuinely different, tree-free and ancestral-reconstruction-free route:
+pairwise genome-to-genome substitution counts pooled by *partition*
+instead of per-site counts pooled across an alignment's columns after
+ancestral state reconstruction. As far as this survey has established
+(Uzzell & Corbin not read, so not fully ruled out), **no tree and no
+ancestral reconstruction of any kind** appears to be a genuinely
+distinguishing property of `PartitionDispersion`, not something it
+inherits from prior art — consistent with it being an independent
+derivation rather than a literature implementation (see previous
+section's caveat, now on firmer footing).
+
+**Why a tree-free estimate is even valid — the identifiability question.**
+Raised directly by the user after reading Gu, Fu & Li (1995), which states
+explicitly: *"one cannot estimate the rate heterogeneity when the number
+of sequences is <3."* This is a real non-identifiability result, not a
+practical inconvenience: with only 2 sequences, a site's history is
+summarized by one binary observation (differ / don't differ), collapsed
+across all sites into two aggregate numbers (transition and transversion
+proportions) — divergence time and among-site rate variance are
+confounded in that pair of numbers, with no way to separate them. Every
+method surveyed above (Gu,Fu&Li, Yang, Wakeley, Gu&Zhang) needs ≥3
+sequences *sharing one evolutionary history* (a tree) specifically to get
+multiple independent looks at the *same* site's rate across different
+lineages, which is what identifies its variance.
+
+`PartitionDispersion` doesn't violate this, because it isn't estimating
+the same quantity. The literature above estimates rate variation **across
+sites, within one shared tree**. `PartitionDispersion` estimates rate
+variation **across partitions, pooled over every genome pair in the
+index** — each genome pair stands in for one independent draw of a
+partition's relative rate, the same identifying role multiple lineages
+play in the tree-based methods, just substituting "genome pair" for
+"lineage" and "partition" for "site." With `P` genome pairs (potentially
+in the hundreds here) all contributing to every partition's pooled count,
+the ≥3-samples identifiability requirement is met by the pair count, not
+by tree depth — closer in spirit to Jin & Nei's own original motivating
+scenario (comparing several genomic *regions* across one fixed panel of
+taxa) than to single-tree site-rate estimation.
+
+This substitution carries its own assumption, which must be named rather
+than left implicit: that a partition's *relative* rate (fast/slow) is
+reasonably stable **across genome pairs**, not just across sites within
+one pair — the direct analogue of the "no lineage-specific rate variation,
+only site-specific" assumption every tree-based `+Γ` method already makes,
+just moved from lineage→pair and site→partition. If this breaks — e.g. an
+index mixing very closely related and very divergent genome pairs, where a
+partition's saturation behavior differs qualitatively between the two —
+the estimate could be misled in a way a tree-based method would at least
+have the topology to detect and a tree-free, pool-everything method
+cannot. Not yet tested against real data with strongly heterogeneous
+pairwise divergence; worth keeping in mind as the main open validity
+question for `--gamma-shape auto`, not the citation question (now settled)
+this whole discussion started from.
+
+**Gu, Fu & Li (1995), read in full — the foundational ML method, not a
+candidate to implement.** Models rate variation as invariant+Γ (`θ` =
+proportion of invariant sites, `α` = gamma shape among the variable sites),
+and derives the exact joint likelihood over nucleotide *configurations*
+across all `n` sequences on a **fixed, a-priori-known tree topology** —
+for 3 sequences this is closed-form (Jukes-Cantor), for `n` sequences it's
+a sum over `4^n` configurations weighted by coefficients tied to the
+specific topology's branching structure (their eq. 32–34), with branch
+lengths, `α`, and `θ` jointly optimized by Newton-Raphson (Hessian over all
+free parameters). The site-rate integral has the same closed form
+`E = θ + (1-θ)(1+D/α)^{-α}` that underlies `+Γ` distance corrections
+generally (the same mechanism as `corrected_log`, not a coincidence — Jin &
+Nei's own formula is the `θ=0` special case applied to a pairwise `D`
+rather than a whole-tree one). Explicitly stated by the authors as
+infeasible beyond about 5–6 taxa (their own simulations stop at 5); not a
+candidate for obikmer's typical genome counts. Notable contribution worth
+keeping in mind regardless: `ρ = (1+θα)/(1+α)` is shown to be a more robust
+summary of rate heterogeneity than `α` alone whenever `θ` isn't negligible
+— `α` and `θ` are confounded (nearly-invariant sites can come from either a
+small `α`'s left tail or from a nonzero `θ`), so `α` alone can be very
+unstable while `ρ` stays well-behaved. `PartitionDispersion` has no `θ`
+term at all (no invariant-site component), so this confound doesn't arise
+for it the same way, but it's worth remembering if an invariant+Γ variant
+is ever considered.
+
+**Gu & Zhang (1997), read in full — rejected for now.** Their procedure:
+(1) a tree with least-squares branch lengths must already exist; (2)
+ancestral states are reconstructed at every internal node (they use Zhang &
+Nei 1997's likelihood method, but note plain parsimony would work, just
+less accurately); (3) per site, branches are split into "changed"/"unchanged"
+given the ancestral reconstruction, and the *expected* substitution count
+`k̂` (corrected for multiple hits) is obtained by solving their eq. (7)
+(or eq. 12 for the generalized model) numerically — a root-find per site,
+not a closed form; (4) the site-level `k̂` values (real-valued, not
+integers) are fit to a negative binomial (Uzzell & Corbin 1971) by ML to
+get α — again no closed form, numerical optimization. This is a real
+per-site pipeline requiring a tree *before* distances/α can be computed —
+the reverse of obikmer's current dependency order (`snp_distance` computes
+distances, which `--nj`/`--upgma` only turn into a tree *afterward*). Their
+own answer to this chicken-and-egg problem is an iterative
+distance→tree→α→distance loop, which they themselves flag as expensive.
+Adopting this would mean a new ancestral-reconstruction subsystem (even a
+parsimony/Fitch-only version) plus a per-family nonlinear solve plus an NB
+ML fit — a much bigger addition than `PartitionDispersion`, not something
+to build without a concrete need beyond "closer to Gu&Zhang's own
+simulated accuracy than our moment estimator."
+
+**Guindon & Gascuel (2002), read in full — plausible future addition, not implemented.**
+Their **Q criterion**: for a candidate α, build a tree (they use BIONJ) from
+the α-corrected distance matrix; for every internal branch, group taxa into
+the 4 subtrees it separates (A, B, C, D), compute mean inter-subtree
+distances \(\bar\delta_{AB}\), \(\bar\delta_{CD}\), etc., and let `S ≤ M ≤ L`
+be the three pairings' sums (`{AB,CD}`, `{AC,BD}`, `{AD,BC}`); the branch's
+reliability score is `Q_branch = L − M` (zero when the four points are
+perfectly tree-additive). The whole-tree `Q` is the mean over internal
+branches (negative/zero branches excluded). `α* = argmin_α Q(α)` over a
+grid (their run: ~60–100 points, 0.1 to 5000, finer spacing where
+sensitivity is highest). Complexity: `O(n²l)` once for the base
+frequencies, then `O(n²r)` for the `r` candidate distance matrices, then
+`O(n³r)` worst case for building `r` trees and evaluating `Q` on each
+(same order as the tree-building step itself, so "free" relative to it).
+
+Central finding, **the opposite direction from an initial misreading in
+this discussion**: `α_opt` (Guindon-Gascuel's topology-optimal value) is
+**always ≥ the true α**, not the reverse — underestimating rate
+heterogeneity (i.e., picking a larger α than reality) reduces the variance
+of distance estimates and thereby improves NJ/BIONJ topological accuracy,
+especially when the molecular clock roughly holds (where `α_opt → ∞`, i.e.
+no correction at all is topologically best). This is an *empirical
+regularity* observed across their tested conditions (η ∈ {0.5, 2.0}, true
+α ∈ {0.1, 0.7}, 20-taxon trees) — not a theorem — so treating it as a hard
+guarantee would be another unverified-claim mistake of the kind this whole
+discussion has been correcting.
+
+**Proposed integration (not yet implemented)**: since `α_true ≤ α_opt`
+empirically, `PartitionDispersion::estimate_alpha()`'s output could seed
+the *lower bound* of Guindon-Gascuel's grid (search `[α̂, ∞)` instead of
+`[0.1, 5000]`), narrowing the grid substantially and reusing the existing
+estimator rather than needing an independent "true α" source — this also
+sidesteps Gu & Zhang's tree-first dependency problem entirely, since
+Guindon-Gascuel's own grid already builds a tree per candidate. Still
+requires: (a) a **BIONJ** implementation (see below), (b) the Q-criterion
+subtree-grouping logic, (c) a grid-search driver. Not started.
+
+**BIONJ availability in the Rust ecosystem — checked 2026-09-11, none found.**
+Searched crates.io directly (API query for "bionj": 0 results) and read the
+docs.rs pages of the most plausible candidates: `phylotree` (builds/reads/
+manipulates trees but has no distance-matrix reconstruction — no NJ, no
+BIONJ, no UPGMA despite web-search summaries claiming otherwise, corrected
+after checking the actual docs), `speedytree` (NJ only — Canonical/RapidNJ/
+Hybrid variants, explicitly no BIONJ), `nj` (plain NJ only), `phylo` (no
+distance-based reconstruction at all). `obikphylo::neighbor_joining`
+(`tree.rs`) is plain NJ, not BIONJ — BIONJ (Gascuel 1997) differs only in
+how branch lengths are weighted during agglomeration (accounts for
+distance-estimate variance, not just the sum), a modest delta over an
+existing NJ implementation rather than a rewrite, but would need to be
+written in-house if the Guindon-Gascuel integration above is ever pursued
+(their own paper notes other tree-building methods gave similar results in
+their simulations, so plain NJ may be an acceptable substitute if BIONJ
+proper is judged not worth the effort).
+
 ### Output format: PHYLIP-relaxed by default for the distance matrix
 
 **Implemented.** The primary distance-matrix output