Correct maximin profiled Gaussian score diagnostics

CHANGELOG.md

@@ -2,6 +2,15 @@
 
 All notable changes to Bayesian-ACh will be documented here.
 
+## Unreleased
+
+- Corrected maximin planning diagnostics for the candidate-specific residual
+  variances used by held-out recovery: the population-optimal asymptotic
+  per-trial quantity is the profiled Gaussian log-score gap `0.5 log1p(a^2 R / sigma^2)`, not a
+  fixed-variance expected log Bayes factor.
+- Renamed exported rate/target fields, retained deprecated Python/CLI aliases,
+  and documented/tested affine reparameterization equivalence.
+
 ## 0.7.0 — 2026-08-23
 
 - Added a transparent finite-grid optimizer for prospective discrimination of

docs/optimal_design.md

@@ -61,13 +61,23 @@ y=a x_k+\epsilon,
 \epsilon\sim\mathcal N(0,\sigma^2),
 ```
 
-the expected held-out log-evidence separation from candidate $`l`$, per trial,
-is
+the generating candidate has residual variance $`\sigma^2`$, whereas the
+alternative's population-optimal profiled residual variance is
+$`\sigma^2+a^2R_{k\mid l}(w)`$. Because the recovery code estimates a separate
+training residual variance for every candidate, the corresponding expected
+held-out Gaussian log-score gap per trial is
 
 ```math
-\frac{a^2}{2\sigma^2}R_{k\mid l}(w).
+G_{k\mid l}(w)
+=
+\frac{1}{2}\log\!\left(
+1+\frac{a^2R_{k\mid l}(w)}{\sigma^2}
+\right).
 ```
 
+The earlier linear expression $`a^2R/(2\sigma^2)`$ is only the first-order
+small-residual expansion of this profiled-variance gap; it is not a Bayes
+factor.
 The primary design criterion is therefore
 
 ```math
@@ -76,8 +86,35 @@ The primary design criterion is therefore
 
 This directly optimizes the worst candidate confusion rather than average
 variance or a global determinant that can hide one nearly indistinguishable
-pair.
+pair. For fixed $`a/\sigma`$, $`G`$ is strictly increasing in $`R`$, so the
+maximin allocation and all residual-ratio comparisons are unchanged by the
+profiled-variance correction.
+
+### Affine-equivalence and identification proposition
+
+For $`\mathrm{Var}_w(x_l)>0`$, the ordered residual has the projection
+interpretation
 
+```math
+R_{k\mid l}(w)
+=
+\min_{b,c}\;\mathbb E_w[(x_k-b-cx_l)^2].
+```
+
+Hence $`R_{k\mid l}=0`$ if and only if $`x_k=b+cx_l`$ almost surely on the
+positive-weight design support. If $`x_l`$ is constant, it adds nothing beyond
+the fitted intercept and the implementation sets
+$`R_{k\mid l}=\mathrm{Var}_w(x_k)`$; a constant generator consequently has
+zero residual against every alternative.
+
+Every candidate recovery fit contains an intercept and a free slope. Replacing
+a candidate column by $`b+c x`$ with $`c\ne0`$ therefore leaves its affine
+column space, fitted predictions, candidate-specific residual variance, and
+held-out Gaussian log score unchanged. Independently z-standardizing the
+declared candidate columns also leaves every residual geometry and the optimized
+allocation unchanged under such affine reparameterizations (including sign
+reversal). This invariance does not justify unequal biological amplitudes:
+$`a`$ remains a prespecified effect per standardized candidate unit.
 ## Integer allocation algorithm
 
 `optimize_maximin_design` uses deterministic greedy allocation followed by
@@ -128,25 +165,26 @@ Gaussian response model or any candidate is biologically correct.
 ## Quantitative trial guidance
 
 The geometry also converts a prespecified signal-to-noise ratio into a transparent
-first-order trial target. For desired expected log Bayes factor $`B`$, the
-worst-pair approximation is
+asymptotic trial target. For desired cumulative expected profiled Gaussian
+log-score gap $`B`$, the worst-pair diagnostic is
 
 ```math
 N_{B}
 =
 \left\lceil
-\frac{2\sigma^2 B}
-     {a^2\min_{k\ne l}R_{k\mid l}(w)}
+\frac{B}
+     {\frac12\log\left(
+       1+a^2\min_{k\ne l}R_{k\mid l}(w)/\sigma^2
+     \right)}
 \right\rceil.
 ```
 
 At unit standardized amplitude, unit noise, and $`B=5`$, the default residuals
-correspond to approximately 40 trials for the maximin design, 88 for the seeded
-uniform factorial design, and 1,112 for the coupled-novelty design. These values
-are planning diagnostics, not retrospective power guarantees: serial dependence,
+correspond to 45 trials for the maximin design, 93 for the seeded uniform
+factorial design, and 1,113 for the coupled-novelty design. These values are
+planning diagnostics, not retrospective power guarantees: serial dependence,
 subject variation, sensor convolution, missing trials, and model misspecification
 must be included in a study-specific simulation before animal numbers are fixed.
-
 ## Scaling assumption and sensitivity requirement
 
 Global standardization gives each computational candidate one unit of variation

scripts/check_markdown_math.py

@@ -9,7 +9,6 @@
 
 from pathlib import Path
 
-
 ROOT = Path(__file__).resolve().parents[1]
 SKIP_PARTS = {".git", ".venv", "build", "dist"}
 

src/bayesian_ach/design.py

@@ -4,6 +4,7 @@
     DesignDiagnostics,
     design_diagnostics,
     pairwise_residual_matrix,
+    profiled_gaussian_log_score_gap,
 )
 from bayesian_ach.design_grid import (
     DESIGN_CANDIDATE_NAMES,
@@ -24,5 +25,6 @@
     "generate_transition_design_grid",
     "optimize_maximin_design",
     "pairwise_residual_matrix",
+    "profiled_gaussian_log_score_gap",
     "uniform_factorial_design",
 ]

src/bayesian_ach/design_benchmark.py

@@ -16,6 +16,7 @@
     generate_transition_design_grid,
     optimize_maximin_design,
     pairwise_residual_matrix,
+    profiled_gaussian_log_score_gap,
     uniform_factorial_design,
 )
 from bayesian_ach.design_recovery import DesignRecoveryRow, recover_design
@@ -36,9 +37,31 @@ class DesignBenchmarkConfig:
     test_fraction: float = 0.35
     effect_size: float = 1.0
     noise_std: float = 1.0
-    target_log_bf: float = 5.0
+    target_log_score_gap: float = 5.0
     max_point_fraction: float = 0.15
     seed: int = 7
+    target_log_bf: float | None = None
+
+    def __post_init__(self) -> None:
+        if self.target_log_bf is None:
+            object.__setattr__(
+                self,
+                "target_log_bf",
+                float(self.target_log_score_gap),
+            )
+
+    @property
+    def resolved_target_log_score_gap(self) -> float:
+        if self.target_log_bf is None:
+            return float(self.target_log_score_gap)
+        if self.target_log_score_gap != 5.0 and not np.isclose(
+            self.target_log_score_gap,
+            self.target_log_bf,
+        ):
+            raise ValueError(
+                "target_log_score_gap and deprecated target_log_bf disagree"
+            )
+        return float(self.target_log_bf)
 
     def validate(self) -> None:
         if self.budget < len(DESIGN_CANDIDATE_NAMES) + 1:
@@ -49,8 +72,9 @@ def validate(self) -> None:
             raise ValueError("test_fraction must lie in (0, 1)")
         if self.effect_size <= 0.0 or self.noise_std <= 0.0:
             raise ValueError("effect_size and noise_std must be positive")
-        if self.target_log_bf <= 0.0:
-            raise ValueError("target_log_bf must be positive")
+        target = self.resolved_target_log_score_gap
+        if not np.isfinite(target) or target <= 0.0:
+            raise ValueError("target_log_score_gap must be finite and positive")
         if not 0.0 < self.max_point_fraction <= 1.0:
             raise ValueError("max_point_fraction must lie in (0, 1]")
 
@@ -84,7 +108,7 @@ def run_design_benchmark(
         max_point_fraction=config.max_point_fraction,
         effect_size=config.effect_size,
         noise_std=config.noise_std,
-        target_log_bf=config.target_log_bf,
+        target_log_score_gap=config.resolved_target_log_score_gap,
     )
     allocations = {
         "coupled_novelty": coupled_novelty_design(rows, config.budget),
@@ -138,7 +162,7 @@ def _design_tables(
             counts,
             effect_size=config.effect_size,
             noise_std=config.noise_std,
-            target_log_bf=config.target_log_bf,
+            target_log_score_gap=config.resolved_target_log_score_gap,
         )
         diagnostics_rows.append({"design": name, **diagnostics.as_dict()})
         geometry = pairwise_residual_matrix(standardized, counts)
@@ -153,9 +177,12 @@ def _design_tables(
                         "generator": generator,
                         "alternative": alternative,
                         "residual_variance": residual,
-                        "expected_log_bf_per_trial": (
-                            config.effect_size**2 * residual
-                            / (2.0 * config.noise_std**2)
+                        "expected_profiled_log_score_gap_per_trial": (
+                            profiled_gaussian_log_score_gap(
+                                residual,
+                                effect_size=config.effect_size,
+                                noise_std=config.noise_std,
+                            )
                         ),
                     }
                 )
@@ -197,7 +224,14 @@ def _summary(
     )
     return {
         "experiment": "prospective_maximin_trial_design",
-        "config": asdict(config),
+        "config": {
+            **{
+                key: value
+                for key, value in asdict(config).items()
+                if key not in {"target_log_bf", "target_log_score_gap"}
+            },
+            "target_log_score_gap": config.resolved_target_log_score_gap,
+        },
         "grid_config": asdict(grid_config),
         "grid_point_count": len(rows),
         "candidate_names": list(DESIGN_CANDIDATE_NAMES),

src/bayesian_ach/design_cli.py

@@ -20,7 +20,13 @@ def _parser() -> argparse.ArgumentParser:
     parser.add_argument("--replicates", type=int, default=200)
     parser.add_argument("--effect-size", type=float, default=1.0)
     parser.add_argument("--noise-std", type=float, default=1.0)
-    parser.add_argument("--target-log-bf", type=float, default=5.0)
+    parser.add_argument("--target-log-score-gap", type=float, default=5.0)
+    parser.add_argument(
+        "--target-log-bf",
+        type=float,
+        default=None,
+        help="Deprecated alias for --target-log-score-gap.",
+    )
     parser.add_argument("--max-point-fraction", type=float, default=0.15)
     parser.add_argument("--seed", type=int, default=7)
     return parser
@@ -34,6 +40,7 @@ def main(argv: Sequence[str] | None = None) -> int:
             replicates_per_generator=args.replicates,
             effect_size=args.effect_size,
             noise_std=args.noise_std,
+            target_log_score_gap=args.target_log_score_gap,
             target_log_bf=args.target_log_bf,
             max_point_fraction=args.max_point_fraction,
             seed=args.seed,