diff --git a/CHANGELOG.md b/CHANGELOG.md index 7de413a..d19be71 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/docs/optimal_design.md b/docs/optimal_design.md index 10b7d0f..b243b5b 100644 --- a/docs/optimal_design.md +++ b/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 diff --git a/scripts/check_markdown_math.py b/scripts/check_markdown_math.py index 24979d1..88c9642 100644 --- a/scripts/check_markdown_math.py +++ b/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"} diff --git a/src/bayesian_ach/design.py b/src/bayesian_ach/design.py index 259cbea..95af0ad 100644 --- a/src/bayesian_ach/design.py +++ b/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", ] diff --git a/src/bayesian_ach/design_benchmark.py b/src/bayesian_ach/design_benchmark.py index 47affab..f372f87 100644 --- a/src/bayesian_ach/design_benchmark.py +++ b/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), diff --git a/src/bayesian_ach/design_cli.py b/src/bayesian_ach/design_cli.py index 67e4b39..249fc59 100644 --- a/src/bayesian_ach/design_cli.py +++ b/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, diff --git a/src/bayesian_ach/design_geometry.py b/src/bayesian_ach/design_geometry.py index a4ddf02..c5353b4 100644 --- a/src/bayesian_ach/design_geometry.py +++ b/src/bayesian_ach/design_geometry.py @@ -1,4 +1,4 @@ -"""Finite-design covariance and expected evidence geometry.""" +"""Finite-design covariance and profiled Gaussian score geometry.""" from __future__ import annotations @@ -20,8 +20,20 @@ class DesignDiagnostics: minimum_pairwise_residual_variance: float covariance_condition_number: float covariance_log_determinant: float - expected_log_bf_per_trial: float - trials_for_expected_log_bf_target: int + expected_profiled_log_score_gap_per_trial: float + trials_for_expected_log_score_gap_target: int + + @property + def expected_log_bf_per_trial(self) -> float: + """Deprecated compatibility alias for the profiled log-score gap.""" + + return self.expected_profiled_log_score_gap_per_trial + + @property + def trials_for_expected_log_bf_target(self) -> int: + """Deprecated compatibility alias for the profiled log-score target.""" + + return self.trials_for_expected_log_score_gap_target def as_dict(self) -> dict[str, float | int]: return asdict(self) @@ -77,9 +89,55 @@ def pairwise_residuals_from_covariance( float(variance[generator]) - float(covariance[generator, alternative]) ** 2 / denominator, ) + else: + # A constant alternative adds nothing beyond the fitted intercept. + result[generator, alternative] = float(variance[generator]) return result +def profiled_gaussian_log_score_gap( + residual_variance: float, + *, + effect_size: float, + noise_std: float, +) -> float: + """Return the population-optimal gap with candidate-specific variance. + + The generating candidate has residual variance sigma squared. An + alternative whose signal leaves projection residual R has profiled + variance sigma squared plus a squared times R. Their expected held-out + Gaussian log-score gap is 0.5 log1p(a squared R / sigma squared). + """ + + residual = float(residual_variance) + amplitude = float(effect_size) + noise = float(noise_std) + if not np.isfinite(residual) or residual < 0.0: + raise ValueError("residual_variance must be finite and nonnegative") + if not np.isfinite(amplitude): + raise ValueError("effect_size must be finite") + if not np.isfinite(noise) or noise <= 0.0: + raise ValueError("noise_std must be finite and positive") + return 0.5 * math.log1p(amplitude**2 * residual / noise**2) + + +def _resolve_log_score_target( + target_log_score_gap: float, + target_log_bf: float | None, +) -> float: + target = float(target_log_score_gap) + if target_log_bf is not None: + alias = float(target_log_bf) + if target != 5.0 and not math.isclose(target, alias): + raise ValueError( + "target_log_score_gap and deprecated target_log_bf disagree" + ) + target = alias + if not np.isfinite(target) or target <= 0.0: + raise ValueError("target_log_score_gap must be finite and positive") + return target + + def diagnostics_from_covariance( covariance: NDArray[np.float64], *, @@ -87,16 +145,20 @@ def diagnostics_from_covariance( support_size: int, effect_size: float, noise_std: float, - target_log_bf: float, + target_log_score_gap: float = 5.0, + target_log_bf: float | None = None, ) -> DesignDiagnostics: - """Summarize identifiability and Gaussian expected-evidence geometry.""" + """Summarize identifiability and population-optimal log-score geometry.""" if trial_count < 1: return DesignDiagnostics( 0, 0, 0.0, 1.0, 0.0, math.inf, -math.inf, 0.0, 2**31 - 1 ) - if noise_std <= 0.0: - raise ValueError("noise_std must be positive") + target = _resolve_log_score_target(target_log_score_gap, target_log_bf) + if not np.isfinite(effect_size): + raise ValueError("effect_size must be finite") + if not np.isfinite(noise_std) or noise_std <= 0.0: + raise ValueError("noise_std must be finite and positive") variance = np.maximum(np.diag(covariance), 0.0) residuals = pairwise_residuals_from_covariance(covariance) off_diagonal = ~np.eye(covariance.shape[0], dtype=bool) @@ -110,11 +172,15 @@ def diagnostics_from_covariance( ) maximum_correlation = float(np.max(np.abs(correlations[off_diagonal]))) eigenvalues = np.maximum(np.linalg.eigvalsh(covariance), 1e-12) - expected_per_trial = effect_size**2 * minimum_residual / (2.0 * noise_std**2) + expected_per_trial = profiled_gaussian_log_score_gap( + minimum_residual, + effect_size=effect_size, + noise_std=noise_std, + ) required = ( 2**31 - 1 if expected_per_trial <= 1e-15 - else int(math.ceil(target_log_bf / expected_per_trial)) + else int(math.ceil(target / expected_per_trial)) ) return DesignDiagnostics( trial_count=trial_count, @@ -124,8 +190,8 @@ def diagnostics_from_covariance( minimum_pairwise_residual_variance=minimum_residual, covariance_condition_number=float(np.max(eigenvalues) / np.min(eigenvalues)), covariance_log_determinant=float(np.sum(np.log(eigenvalues))), - expected_log_bf_per_trial=expected_per_trial, - trials_for_expected_log_bf_target=required, + expected_profiled_log_score_gap_per_trial=expected_per_trial, + trials_for_expected_log_score_gap_target=required, ) @@ -135,7 +201,8 @@ def design_diagnostics( *, effect_size: float = 1.0, noise_std: float = 1.0, - target_log_bf: float = 5.0, + target_log_score_gap: float = 5.0, + target_log_bf: float | None = None, ) -> DesignDiagnostics: """Evaluate one allocation without expanding individual trials.""" @@ -153,6 +220,7 @@ def design_diagnostics( support_size=int(np.sum(allocation > 0)), effect_size=effect_size, noise_std=noise_std, + target_log_score_gap=target_log_score_gap, target_log_bf=target_log_bf, ) diff --git a/src/bayesian_ach/design_optimizer.py b/src/bayesian_ach/design_optimizer.py index 9f5c29e..abac614 100644 --- a/src/bayesian_ach/design_optimizer.py +++ b/src/bayesian_ach/design_optimizer.py @@ -41,7 +41,8 @@ def optimize_maximin_design( exchange_passes: int = 3, effect_size: float = 1.0, noise_std: float = 1.0, - target_log_bf: float = 5.0, + target_log_score_gap: float = 5.0, + target_log_bf: float | None = None, ) -> OptimizedDesign: """Allocate trials greedily, then refine by deterministic one-for-one swaps.""" @@ -71,6 +72,7 @@ def evaluate( support_size=support_size, effect_size=effect_size, noise_std=noise_std, + target_log_score_gap=target_log_score_gap, target_log_bf=target_log_bf, ) diff --git a/tests/test_design.py b/tests/test_design.py index cae5532..38c4454 100644 --- a/tests/test_design.py +++ b/tests/test_design.py @@ -9,9 +9,11 @@ generate_transition_design_grid, optimize_maximin_design, pairwise_residual_matrix, + profiled_gaussian_log_score_gap, uniform_factorial_design, ) from bayesian_ach.design_benchmark import DesignBenchmarkConfig, run_design_benchmark +from bayesian_ach.design_recovery import _fit_and_score def test_design_grid_is_finite_and_dissociates_all_candidates() -> None: @@ -64,6 +66,103 @@ def test_equal_budget_recovery_beats_coupled_novelty() -> None: assert result.summary["optimized_over_novelty_residual_ratio"] > 10.0 +def test_profiled_gaussian_gap_and_default_trial_targets() -> None: + rows, _, standardized = generate_transition_design_grid() + budget = 60 + optimized = optimize_maximin_design(standardized, budget) + diagnostics = { + "maximin_optimized": optimized.diagnostics, + "uniform_factorial": design_diagnostics( + standardized, + uniform_factorial_design(len(rows), budget, seed=7), + ), + "coupled_novelty": design_diagnostics( + standardized, + coupled_novelty_design(rows, budget), + ), + } + + assert profiled_gaussian_log_score_gap( + 0.25, + effect_size=1.0, + noise_std=1.0, + ) == pytest.approx(0.5 * np.log1p(0.25)) + assert { + name: value.trials_for_expected_log_score_gap_target + for name, value in diagnostics.items() + } == { + "maximin_optimized": 45, + "uniform_factorial": 93, + "coupled_novelty": 1113, + } + assert optimized.diagnostics.expected_log_bf_per_trial == ( + optimized.diagnostics.expected_profiled_log_score_gap_per_trial + ) + legacy = design_diagnostics( + standardized, + optimized.counts, + target_log_bf=5.0, + ) + assert legacy.trial_count == optimized.diagnostics.trial_count + assert legacy.trials_for_expected_log_score_gap_target == ( + optimized.diagnostics.trials_for_expected_log_score_gap_target + ) + assert legacy.minimum_pairwise_residual_variance == pytest.approx( + optimized.diagnostics.minimum_pairwise_residual_variance, + abs=1e-12, + ) + assert legacy.expected_profiled_log_score_gap_per_trial == pytest.approx( + optimized.diagnostics.expected_profiled_log_score_gap_per_trial, + abs=1e-12, + ) + assert "expected_log_bf_per_trial" not in legacy.as_dict() + + +def test_pairwise_residual_identifies_affine_equivalence_on_support() -> None: + alternative = np.array([-1.0, 0.0, 1.0, 2.0]) + affine_generator = 2.0 + 3.0 * alternative + nonaffine_generator = np.array([0.0, 1.0, 0.0, 2.0]) + constant_alternative = np.ones(4) + signals = np.column_stack( + (alternative, affine_generator, nonaffine_generator, constant_alternative) + ) + residual = pairwise_residual_matrix(signals, np.ones(4, dtype=np.int64)) + + assert residual[1, 0] == pytest.approx(0.0, abs=1e-12) + assert residual[2, 0] > 0.0 + assert residual[2, 3] == pytest.approx(np.var(nonaffine_generator)) + assert residual[3, 0] == pytest.approx(0.0, abs=1e-12) + +def test_affine_reparameterization_preserves_geometry_and_profiled_scores() -> None: + _, _, standardized = generate_transition_design_grid() + shifts = np.linspace(-4.0, 3.0, standardized.shape[1]) + scales = np.array([0.5, -1.5, 2.0, -0.75, 3.0, -2.5]) + transformed = standardized * scales + shifts + transformed = (transformed - transformed.mean(axis=0)) / transformed.std(axis=0) + baseline = optimize_maximin_design(standardized, 24, exchange_passes=1) + affine = optimize_maximin_design(transformed, 24, exchange_passes=1) + + np.testing.assert_array_equal(affine.counts, baseline.counts) + np.testing.assert_allclose( + pairwise_residual_matrix(transformed, affine.counts), + pairwise_residual_matrix(standardized, baseline.counts), + atol=1e-12, + ) + + trial_signals = standardized[np.repeat(np.arange(len(baseline.counts)), baseline.counts)] + response = 0.7 * trial_signals[:, 2] + np.linspace(-0.2, 0.2, len(trial_signals)) + train = np.arange(0, 16, dtype=np.int64) + test = np.arange(16, len(trial_signals), dtype=np.int64) + winner, margin = _fit_and_score(trial_signals, response, train, test) + affine_winner, affine_margin = _fit_and_score( + trial_signals * scales + shifts, + response, + train, + test, + ) + assert affine_winner == winner + assert affine_margin == pytest.approx(margin, abs=1e-10) + def test_invalid_budget_and_allocation_are_rejected() -> None: _, _, standardized = generate_transition_design_grid( TransitionDesignGridConfig() diff --git a/tests/test_design_cli.py b/tests/test_design_cli.py index 90c821d..8a58a35 100644 --- a/tests/test_design_cli.py +++ b/tests/test_design_cli.py @@ -14,6 +14,8 @@ def test_design_cli_writes_complete_evidence(tmp_path: Path) -> None: "24", "--replicates", "12", + "--target-log-score-gap", + "3.0", "--seed", "5", ] @@ -22,6 +24,29 @@ def test_design_cli_writes_complete_evidence(tmp_path: Path) -> None: summary = json.load(handle) assert summary["experiment"] == "prospective_maximin_trial_design" assert summary["designs"][0]["trial_count"] == 24 + assert summary["config"]["target_log_score_gap"] == 3.0 + assert "target_log_bf" not in summary["config"] + assert "expected_profiled_log_score_gap_per_trial" in summary["designs"][0] assert (output / "design_allocation.csv").is_file() assert (output / "design_pairwise_geometry.csv").is_file() assert (output / "design_optimization_trace.csv").is_file() + + +def test_design_cli_retains_deprecated_target_alias(tmp_path: Path) -> None: + output = tmp_path / "legacy-design" + assert main( + [ + "--output", + str(output), + "--budget", + "12", + "--replicates", + "1", + "--target-log-bf", + "2.5", + ] + ) == 0 + with (output / "summary.json").open(encoding="utf-8") as handle: + summary = json.load(handle) + assert summary["config"]["target_log_score_gap"] == 2.5 + assert summary["designs"][0]["trials_for_expected_log_score_gap_target"] > 0