diff --git a/CHANGELOG.md b/CHANGELOG.md index 1b00edb..b9ec462 100644 --- a/CHANGELOG.md +++ b/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=` 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 diff --git a/analysis/trade_study/_world.py b/analysis/trade_study/_world.py index 5946416..bcdc6ef 100644 --- a/analysis/trade_study/_world.py +++ b/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", diff --git a/analysis/trade_study/validate_shipped_defaults.py b/analysis/trade_study/validate_shipped_defaults.py new file mode 100644 index 0000000..50efc24 --- /dev/null +++ b/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() diff --git a/src/vbpca_py/defaults.py b/src/vbpca_py/defaults.py index 160fa50..957721a 100644 --- a/src/vbpca_py/defaults.py +++ b/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.::