docs: re-validate shipped defaults.py hyperparameters with replication (#111)

CHANGELOG.md

@@ -13,6 +13,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ### Changed
 - **Behavior change:** the default (`random_state=None`) now draws fresh entropy on every `fit()` call. Previously, default initialization was silently seeded with a fixed value regardless of configuration, so repeated fits produced identical results without any way to request a different draw. Pass `random_state=<int>` for reproducible runs (#109).
 - `recommend_config()`'s `missingness` parameter now warns (`UserWarning`) when passed anything other than the default `"auto"`, instead of silently ignoring it. Recommendations are still bucketed by `p` only — the Option A trade study's example recommendations are too sparse per (p-bucket, missingness) cell (23 points across 3 p-buckets x 4 missingness categories) to bucket on responsibly without shipping unreplicated values (#110, see also #111).
+- `defaults.py`'s docstring corrected the unsupported "`hp_va` is the dominant lever" claim (the trade study's own marginal sensitivity data doesn't support it) and now documents a real validation: replicated (n_reps=8, seeded) rank_mae for the shipped bucket configs is 28-58% lower than the library default across all three p-buckets, at a 0.4-3.8% cost in holdout RMSE (#111).
 
 ## [0.3.0] - 2026-08-17
 

analysis/trade_study/_world.py

@@ -56,10 +56,26 @@ class VBPCASimulator:
     def __init__(self, *, regime_defaults: dict[str, Any] | None = None) -> None:
         self._regime_defaults: dict[str, Any] = regime_defaults or {}
 
-    def generate(self, config: dict[str, Any]) -> tuple[Any, Any]:
-        """Run a single trial and return (truth, observations)."""
+    def generate(self, config: dict[str, Any], *, rep: int = 0) -> tuple[Any, Any]:
+        """Run a single trial and return (truth, observations).
+
+        ``rep`` opts into trade-study's replicated-trials convention
+        (jcm-sci/trade-study#112): each replicate derives its own data
+        seed and VBPCA ``random_state`` from ``rep``, so ``run_grid(...,
+        n_reps=N)`` produces genuinely independent draws instead of N
+        copies of the same trial. ``rep=0`` (the default) reproduces the
+        previous single-seed behavior for data generation, but note that
+        VBPCA's own default init (``random_state=None`` since #109) draws
+        fresh entropy regardless — this method always passes an explicit
+        ``random_state`` to keep trials reproducible.
+        """
         config = {**FIXED_STRUCTURAL, **self._regime_defaults, **config}
-        rng = np.random.default_rng(config.get("seed", RNG_SEED))
+        base_seed = int(config.get("seed", RNG_SEED))
+        # Distinct large odd multipliers keep the data-generation and
+        # model-init seed sequences independent across replicates.
+        data_seed = base_seed + 1_000_003 * rep
+        init_seed = base_seed + 2_000_003 * rep
+        rng = np.random.default_rng(data_seed)
 
         # ── Data regime ─────────────────────────────────────────
         n: int = config["n"]
@@ -73,7 +89,7 @@ def generate(self, config: dict[str, Any]) -> tuple[Any, Any]:
         train_mask, holdout_mask = holdout_split(obs_mask, HOLDOUT_FRACTION, rng)
 
         # ── Map tunable factors → VBPCA kwargs ─────────────────
-        vbpca_kw: dict[str, Any] = {"verbose": 0}
+        vbpca_kw: dict[str, Any] = {"verbose": 0, "random_state": init_seed}
         for key in (
             "hp_va",
             "hp_vb",

analysis/trade_study/validate_shipped_defaults.py

@@ -0,0 +1,177 @@
+#!/usr/bin/env python
+"""Replicated validation of the exact bucket configs shipped in defaults.py.
+
+Targeted follow-up to the Option A trade study (#111): rather than
+re-running the full Sobol collection, this compares the library's raw
+default configuration against the *exact* configs ``recommend_config()``
+returns today — the values ``v3_compare.py``'s replicated comparison never
+actually validated, since it reads from older per-family JSON artifacts
+instead of the shipped module.
+
+Each regime (training + held-out validation) is evaluated ``--n-reps``
+times per condition, with both the data draw and the VBPCA ``random_state``
+varied per replicate (via ``VBPCASimulator.generate(..., rep=...)``), using
+trade-study's replicated-trials support (jcm-sci/trade-study#112). Both
+conditions share the same base seed per regime, so a given (regime, rep)
+pair draws identical synthetic data and VBPCA init noise in both
+conditions -- a paired (common-random-numbers) design that isolates the
+effect of the config choice itself from other randomness sources.
+
+Usage
+-----
+    python -m analysis.trade_study.validate_shipped_defaults --n-reps 8
+"""
+
+from __future__ import annotations
+
+import argparse
+import json
+import pathlib
+from typing import Any
+
+from trade_study import Direction, Observable, run_grid
+
+from vbpca_py import recommend_config
+
+from ._common import VALIDATION_REGIMES
+from ._world import VBPCAScorer, VBPCASimulator
+from .option_a_pipeline import TRAINING_REGIMES
+
+RESULTS_DIR = pathlib.Path("analysis/results/optionA")
+
+OBSERVABLES: list[Observable] = [
+    Observable("rank_mae", Direction.MINIMIZE),
+    Observable("rank_under", Direction.MINIMIZE),
+    Observable("rank_over", Direction.MINIMIZE),
+    Observable("holdout_rmse", Direction.MINIMIZE),
+    Observable("coverage_95", Direction.MAXIMIZE),
+]
+
+CONDITIONS = ("default", "shipped")
+
+# Display-only bucket label, mirroring recommend_config's own documented
+# p-bucketing (p<=30 smallp, p<=70 trans, else large). Not used for
+# behavior -- the "shipped" condition always calls the public
+# recommend_config(n, p) directly.
+_SMALLP_MAX_P = 30
+_TRANS_MAX_P = 70
+
+
+def _bucket_label(p: int) -> str:
+    """Return the display bucket label for ``p`` features.
+
+    Returns:
+        One of "smallp", "trans", "large".
+    """
+    if p <= _SMALLP_MAX_P:
+        return "smallp"
+    if p <= _TRANS_MAX_P:
+        return "trans"
+    return "large"
+
+
+def _grid_for_condition(
+    regimes: list[dict[str, Any]], condition: str, seed: int
+) -> list[dict[str, Any]]:
+    """Build a run_grid grid: one row per regime, tagged with its bucket.
+
+    The "default" condition supplies only regime keys, so VBPCASimulator
+    falls back to the library's own untuned defaults. The "shipped"
+    condition merges in whatever ``recommend_config(n, p)`` returns today.
+    Each row gets a distinct base ``seed`` (offset by its index) so
+    different regimes don't share a replicate seed sequence.
+
+    Returns:
+        List of config dicts suitable for run_grid.
+    """
+    grid: list[dict[str, Any]] = []
+    for idx, regime in enumerate(regimes):
+        cfg = dict(regime)
+        n, p = int(regime["n"]), int(regime["p"])
+        cfg["_bucket"] = _bucket_label(p)
+        cfg["seed"] = seed + idx
+        if condition == "shipped":
+            cfg.update(recommend_config(n=n, p=p))
+        grid.append(cfg)
+    return grid
+
+
+def run_validation(n_reps: int = 8, seed: int = 42, n_jobs: int = -1) -> dict[str, Any]:
+    """Run the replicated default-vs-shipped comparison.
+
+    Returns:
+        Dict with per-condition, per-bucket aggregated results, saved
+        alongside the return value at
+        ``analysis/results/optionA/shipped_defaults_validation.json``.
+    """
+    regimes = TRAINING_REGIMES + VALIDATION_REGIMES
+    world = VBPCASimulator()
+    scorer = VBPCAScorer()
+
+    per_condition: dict[str, list[dict[str, Any]]] = {}
+    for condition in CONDITIONS:
+        grid = _grid_for_condition(regimes, condition, seed)
+        # run_grid doesn't know about "_bucket"; VBPCASimulator.generate
+        # ignores unrecognised keys, so it's safe to carry through the
+        # grid purely for post-hoc grouping below.
+        table = run_grid(
+            world,
+            scorer,
+            grid,
+            OBSERVABLES,
+            n_jobs=n_jobs,
+            n_reps=n_reps,
+        )
+        agg = table.aggregate_replicates()
+
+        rows: list[dict[str, Any]] = []
+        for i, cfg in enumerate(agg.configs):
+            row = {
+                "condition": condition,
+                "bucket": cfg["_bucket"],
+                "n": cfg["n"],
+                "p": cfg["p"],
+                "true_rank": cfg["true_rank"],
+                "missingness": cfg["missingness"],
+                "n_reps": agg.metadata[i]["n_reps"],
+            }
+            for j, name in enumerate(agg.observable_names):
+                row[name] = float(agg.scores[i, j])
+                row[f"{name}_std"] = float(agg.metadata[i]["score_std"][name])
+            rows.append(row)
+        per_condition[condition] = rows
+
+    RESULTS_DIR.mkdir(parents=True, exist_ok=True)
+    out = RESULTS_DIR / "shipped_defaults_validation.json"
+    out.write_text(json.dumps(per_condition, indent=2))
+
+    _print_summary(per_condition)
+    print(f"\nSaved -> {out}")
+    return per_condition
+
+
+def _print_summary(per_condition: dict[str, list[dict[str, Any]]]) -> None:
+    buckets = sorted({row["bucket"] for row in per_condition["default"]})
+    print(f"\n{'bucket':10s} {'condition':10s} {'rank_mae':>10s} {'rmse':>8s}")
+    print("-" * 42)
+    for bucket in buckets:
+        for condition in CONDITIONS:
+            rows = [r for r in per_condition[condition] if r["bucket"] == bucket]
+            if not rows:
+                continue
+            mae = sum(r["rank_mae"] for r in rows) / len(rows)
+            rmse = sum(r["holdout_rmse"] for r in rows) / len(rows)
+            print(f"{bucket:10s} {condition:10s} {mae:10.3f} {rmse:8.3f}")
+
+
+def main() -> None:
+    parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument("--n-reps", type=int, default=8)
+    parser.add_argument("--seed", type=int, default=42)
+    parser.add_argument("--n-jobs", type=int, default=-1)
+    args = parser.parse_args()
+    run_validation(n_reps=args.n_reps, seed=args.seed, n_jobs=args.n_jobs)
+
+
+if __name__ == "__main__":
+    main()

src/vbpca_py/defaults.py

@@ -5,11 +5,35 @@
 reconstruction error within a small tolerance of the library defaults.
 
 These configurations come from the Option A regime-surrogate trade study
-(``analysis/trade_study``).  The dominant lever is a **strong ARD loadings
-prior** (``hp_va`` ~ 0.7 versus the library default 0.001): it drives correct
-component pruning, whereas the weak default prior under-prunes and recovers the
-true rank only about a third of the time.  Recommendations are bucketed by the
-feature count ``p`` — the study's primary axis of rank-recovery difficulty.
+(``analysis/trade_study``), which fits a random-forest surrogate over the
+full joint factor space per regime and recommends the joint-optimal
+config. The joint-optimal configs use a moderately strong ARD loadings
+prior (``hp_va`` ~ 0.65-0.75 vs. the library default 0.001) alongside a
+small xprobe fraction (~0.01-0.02); the weak default prior under-prunes
+and recovers the true rank only about a third of the time.
+
+**On "dominant lever" claims:** a marginal (univariate Spearman)
+sensitivity analysis of the same trade study's data does *not* single out
+``hp_va`` as dominant — ``xprobe_fraction`` is the strongest, most
+significant, and most consistent per-bucket predictor of ``rank_mae`` in
+that marginal view, and ``hp_va``'s marginal correlation is weak and not
+statistically significant in 2 of the 3 p-buckets (though its tercile
+means do show a real, nonlinear U-shape a monotonic Spearman correlation
+understates). The RF surrogate optimises interactions a marginal view
+can't see, so which single factor (if any) is "dominant" remains
+unreconciled; treat the recommendation as a joint-optimal bundle rather
+than attributing its effect to any one factor.
+
+**Validation status (#111):** the exact configs this module ships *are*
+now validated with real replication and seeded VBPCA initialization
+(``analysis/trade_study/validate_shipped_defaults.py``, n_reps=8 across
+training + held-out regimes) — replicated ``rank_mae`` for the shipped
+config is 28-58% lower than the library default across all three
+p-buckets (smallp 0.34→0.14, trans 1.20→0.74, large 1.44→1.04), at a
+small cost in holdout RMSE (+0.4-3.8%). What remains unvalidated is the
+*attribution* to a specific factor, not whether the shipped bundle helps.
+Recommendations are bucketed by the feature count ``p`` — the study's
+primary axis of rank-recovery difficulty.
 
 The returned dict is intended to be splatted into the estimator, e.g.::