From c8c6c50e81021cd626a2459f1a3e5716b7b5e054 Mon Sep 17 00:00:00 2001 From: ***** <721466+soodoku@users.noreply.github.com> Date: Sat, 8 Aug 2026 23:37:39 -0700 Subject: [PATCH 1/4] Keep the interval endpoints, and gate width and power Coverage cannot see a vacuous interval. `assert_coverage` is satisfied by one so wide it always covers, because a rate of 1.0 sits inside the binomial band at any study smaller than about sixty replicates, and `coverage > 0.9` written by hand is satisfied by it always. Three consuming repositories hit this independently and each left a comment where a test should have been; one of them after an inflation heuristic drove a reported standard error to 3e7 times the estimation error while coverage stayed high. All three are reproduced as tests here. The runner computed the endpoints, used them once and dropped them, so the evidence was being thrown away. It now keeps them, and `MonteCarloResult` gained widths. `assert_intervals_informative` fails only when the width exceeds what the study can still see fail *and* the study never saw it fail. The conjunction is the point: a t interval at n=5 is 1.33 times the oracle width and correct, so width alone cannot separate conservatism from vacuity. The reference width is measured by the study rather than chosen; what remains a convention is written down in `vacuous_width_ratio`'s docstring rather than presented as a derivation. `assert_power` and `assert_more_powerful` fill the hole under the package's own claim to answer whether a test has power. Power is a floor, so the first is one-sided; the second replaces `a.rejection_rate > b.rejection_rate`, which is satisfied by a gap of one replicate in four hundred. `assert_se_calibrated`'s tolerance was 0.15, the one number here chosen by hand. The sampling distribution of the ratio gives it: three sigma is 0.21 at 100 replicates and 0.05 at 2000, so the fixed value was tight enough to fail correct estimators in a fast tier and loose enough to certify a 12% error in a deep one. It is now derived by default. geoinference is the only consumer relying on the default and its suite still passes: 43 passed, 4 subtests passed. Co-Authored-By: Claude Opus 5 (1M context) --- CHANGELOG.md | 68 +++- README.md | 94 ++++- src/simcheck/__init__.py | 20 + src/simcheck/gates.py | 536 ++++++++++++++++++++++++- src/simcheck/results.py | 105 ++++- src/simcheck/runner.py | 11 +- tests/test_against_known_statistics.py | 163 ++++++++ tests/test_negative.py | 478 +++++++++++++++++++++- tests/test_runner.py | 33 ++ 9 files changed, 1486 insertions(+), 22 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 1b759e4..5617d46 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] +## [0.1.0] - 2026-08-08 + ### Added - Initial extraction from `incline/tests/_statistics.py`, generalised so it @@ -30,16 +32,78 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Tier helpers `reps_for`, `deep_tier`, `FAST_REPS` and `DEEP_REPS`, driven by `SIMCHECK_DEEP` and `SIMCHECK_REPS`. +- **Interval endpoints are kept, not just whether each interval covered.** + `MonteCarloResult` gained `lowers`, `uppers`, `widths`, `mean_width` and + `median_width`; `monte_carlo` passes the endpoints through instead of computing + them, using them once and dropping them. Supplying one endpoint without the + other, an upper endpoint below its lower one, or a `covered` array that + contradicts the endpoints are all rejected; supplying endpoints without + `covered` fills `covered` in, since the endpoints and the truth determine it. + +- `assert_intervals_informative`, which fails an interval so wide that the study + never saw it miss. `assert_coverage` cannot: a coverage rate of 1.0 sits inside + the binomial band at any study smaller than about 60 replicates, and a + hand-written `coverage > 0.9` is satisfied by an interval that always covers. + Three consuming repositories hit this independently and left a comment where a + test should have been — one of them after an inflation heuristic drove a + reported standard error to 3e7 times the estimation error while coverage stayed + high. The gate fails only when the width exceeds `vacuous_width_ratio(nominal, + reps)` *and* the study observed fewer than three misses, so a correct but + conservative procedure — a t interval at n=5 is 1.33 times the oracle width — + is not flagged. + +- `vacuous_width_ratio` and `width_ratio`, the threshold and the measured + quantity behind that gate. The reference width, `2 * z * sampling_sd`, is + measured by the study, so no absolute width appears anywhere. How many expected + misses per study counts as "could not have failed" is a convention rather than + a derivation, and `vacuous_width_ratio`'s docstring says so and says why: correct + procedures occupy the whole range of widths above the oracle, so no sampling + distribution separates conservatism from vacuity on width alone. + +- `assert_narrower`, for the efficiency half of an interval comparison, banded by + the Monte Carlo standard error of the difference in mean width. + +- `assert_power` and `assert_more_powerful`. The package documented power as one + of the four questions it answers and had no gate for it; consumers were + reaching for `assert_proportion`, which is two-sided and needs a nominal you + already know analytically, or hand-rolling a two-sample standard error. + `assert_power` is one-sided, because power is a floor and a two-sided band + fails a test for being better than claimed. + - **Negative tests for every gate.** Each is exercised on input that satisfies its property, where it must stay silent, and on input that violates it, where it must raise. Plus gates run against estimators whose behaviour is known analytically: the sample mean trips nothing, the `ddof=0` variance is caught with its textbook `-sigma^2/n` bias, a 1.96 interval at n=5 is caught - under-covering at about 0.875, and a false-positive check confirms correct - estimators are essentially never flagged. + under-covering at about 0.875, a two-sided z test shows the power its formula + gives (0.323 at n=25 and 0.851 at n=100 for delta=0.3), an interval built at a + known scale comes out at a width ratio of exactly one, and a false-positive + check confirms correct estimators are essentially never flagged. + +### Changed + +- **`assert_se_calibrated`'s tolerance now comes from the replicate count.** It + was `0.15`, the one number in the package chosen by hand rather than derived, + and it was wrong in both directions: `se_ratio` divides a mean of `reps` + reported standard errors by a sample standard deviation of `reps` estimates, so + its Monte Carlo spread is `sqrt(cv^2/reps + 1/(2(reps-1)))`, and three of those + is 0.21 at 100 replicates and 0.05 at 2000. The fixed value was therefore tight + enough to fail correct estimators in a fast tier and loose enough to certify a + 12% error in a deep one. Passing `tolerance=` explicitly still overrides it, + and `se_ratio_tolerance(result)` returns the derived band. + + This changes behaviour for callers that relied on the default. Of the six + consuming repositories only `geoinference` does, and its suite was re-run + against this branch: 43 passed, 4 subtests passed. Every other consumer passes + `tolerance=` explicitly or does not call the gate. ### Fixed +- **`assert_se_calibrated` diagnosed a missing standard error as a constant + estimator.** `Estimate` documents that leaving `standard_error` as NaN means + the gate "will have nothing to check and will say so"; it said the estimator + did not vary across replicates, which is a different and false diagnosis. + - **The extracted `assert_rate` could report the worst possible result as the best possible one.** It took either a count or a rate and guessed which: `observed = successes / reps if successes > 1 else float(successes)`. For a diff --git a/README.md b/README.md index e9aceee..a6d30e8 100644 --- a/README.md +++ b/README.md @@ -9,6 +9,8 @@ under the null, and does it have power under an alternative. ```bash pip install simcheck +# until the first PyPI release lands: +pip install "simcheck @ git+https://github.com/finite-sample/simcheck@v0.1.0" ``` Python 3.11+. Depends on numpy, and nothing else. @@ -86,12 +88,74 @@ rate. |---|---| | `assert_unbiased` | the mean estimate is more than 3 Monte Carlo standard errors from the truth | | `assert_coverage` | interval coverage falls outside the binomial band around the nominal level | +| `assert_intervals_informative` | the intervals are so wide the study never saw one miss, so their coverage measures the width | +| `assert_narrower` | one method's intervals are not measurably narrower than another's | | `assert_se_calibrated` | the reported standard error misstates the spread actually observed | -| `assert_proportion` | an observed **rate** — size, power — is inconsistent with the claimed one | +| `assert_power` | a test rejects less often than claimed under an alternative (one-sided) | +| `assert_more_powerful` | one test does not reject measurably more often than another at the same alternative | +| `assert_proportion` | an observed **rate** — size, coverage — is inconsistent with the claimed one | | `assert_count_rate` | the same, given a **count** of successes | -`binomial_band(nominal, reps)` gives the band directly if you want to assert -something else against it. +`binomial_band(nominal, reps)`, `vacuous_width_ratio(nominal, reps)` and +`se_ratio_tolerance(result)` give the three thresholds directly, if you want to +assert something else against them. `width_ratio(result, nominal)` is the +measured quantity behind the first of those. + +## Coverage cannot see a vacuous interval + +`assert_coverage` is satisfied by an interval so wide it always covers, whenever +the study is small enough that a rate of 1.0 still sits inside the binomial band +— and `assert coverage > 0.9`, written by hand, is satisfied by it always. +Three repositories hit this independently and each left a comment where a test +should have been; one had shipped an inflation heuristic that drove the reported +standard error to 3×10⁷ times the estimation error while coverage stayed high, +because a vacuous interval covers everything. + +The fix is to keep the endpoints, not just the hit: + +```python +result = monte_carlo(replicate, truth=2.0, reps=2000, seed=11) +assert_coverage(result, 0.95, "t interval") +assert_intervals_informative(result, 0.95, "t interval") +``` + +Two things must both be true before that fails, and the conjunction is the +point. The interval must be far wider than the width its own level requires +against the spread the estimator actually has — `2 * z * sampling_sd`, measured +by the study, so no absolute width is written down anywhere. *And* the study must +never have seen it miss. A Student t interval at n=5 is 1.33 times the normal +oracle width and an anytime-valid interval is wider still; both are correct, both +miss at their nominal rate, and the study watches them do it. Width alone cannot +tell conservatism from vacuity. Width plus a study that never saw a failure can. + +`vacuous_width_ratio(0.95, reps)` is the width multiple at which a study of that +size stops being able to observe a miss: 1.78 at 100 replicates, 1.96 at 400, +2.15 at 2000. It rises with the replicate count, which inverts the usual +direction and is meant to — more replicates resolve rarer failures. + +The one thing here that is a judgement rather than a derivation is how many +expected misses per study counts as "could not have failed". Its docstring says +so, and says why no sampling distribution fixes it: correct procedures occupy the +whole range of widths above the oracle. + +## Power + +The package claims to answer whether a test has power under an alternative, so +there is a gate for it: + +```python +assert_power(result, 0.80, "score test at delta=0.3") +assert_more_powerful(robust, naive, "robust against naive at the same alternative") +``` + +`assert_power` is one-sided, unlike `assert_proportion`: power is a floor, and a +two-sided band would fail a test for being *better* than promised. Size, which is +a target rather than a floor, still belongs in `assert_proportion`. + +`assert_more_powerful` bands the gap between two rejection rates by the standard +error of the difference. The assertion it replaces — `a.rejection_rate > +b.rejection_rate` — passes on a gap of one replicate in four hundred and reports +whichever method the seed favoured as the winner. ## Two failure modes this is built to prevent @@ -118,6 +182,22 @@ reported as the best possible one. That is why counts and rates are separate functions here, and why `assert_proportion` raises rather than guesses when it is handed something outside `[0, 1]`. +## The rule applies to simcheck too + +`assert_se_calibrated` used to take `tolerance=0.15`, which was the one number in +the package chosen by hand rather than derived — and it was wrong in both +directions at once. `se_ratio` is `mean(reported se) / sd(estimates)`, and both +halves are estimated from the same replicates, so it is noisy even when the +estimator is perfect: the numerator's relative standard error is `cv/sqrt(reps)` +and the denominator's is `1/sqrt(2(reps-1))`. Added in quadrature and taken at +three sigma, that is 0.21 at 100 replicates and 0.05 at 2000. A fixed 0.15 was +therefore tight enough to fail correct estimators in a fast tier and loose enough +to certify a 12% error in a deep one. + +The tolerance is now derived from `reps` by default; passing a number still +overrides it, which is worth doing when the claim really is about a fixed +accuracy at a fixed sample size. `se_ratio_tolerance(result)` returns the band. + ## Negative tests Every gate has one: an input that violates the property, and a check that the @@ -129,9 +209,11 @@ itself as tested. The gates are also run against estimators whose behaviour is known analytically (`tests/test_against_known_statistics.py`): the sample mean must trip nothing, the uncorrected `ddof=0` variance must be caught with its textbook bias of -`-σ²/n`, and a 1.96 interval at n=5 must be caught under-covering at ≈0.875. -There is also a false-positive check, because a gate that fires on 5% of correct -code gets disabled within a week. +`-σ²/n`, a 1.96 interval at n=5 must be caught under-covering at ≈0.875, a +two-sided z test must show the power its formula gives (0.323 at n=25, 0.851 at +n=100 for δ=0.3), and an interval built at a *known* scale must come out at a +width ratio of exactly one. There is also a false-positive check, because a gate +that fires on 5% of correct code gets disabled within a week. ## Tiers diff --git a/src/simcheck/__init__.py b/src/simcheck/__init__.py index 0085d1c..7e16583 100644 --- a/src/simcheck/__init__.py +++ b/src/simcheck/__init__.py @@ -26,6 +26,12 @@ helper that silently passes everything is worse than no helper, because it converts an untested codebase into one that reports itself as tested. +**A gate that a vacuous answer satisfies.** Coverage is satisfied by an interval +so wide it always covers, which is why the endpoints are kept and +:func:`~simcheck.assert_intervals_informative` exists; and one-sided power was +being asserted by hand as ``a > b``, which is satisfied by a gap of one +replicate, which is why :func:`~simcheck.assert_more_powerful` exists. + Examples: >>> import numpy as np >>> from simcheck import MonteCarloResult, assert_coverage, assert_unbiased @@ -51,10 +57,17 @@ GATE_SIGMAS, assert_count_rate, assert_coverage, + assert_intervals_informative, + assert_more_powerful, + assert_narrower, + assert_power, assert_proportion, assert_se_calibrated, assert_unbiased, binomial_band, + se_ratio_tolerance, + vacuous_width_ratio, + width_ratio, ) from .results import MonteCarloResult from .runner import Estimate, monte_carlo @@ -74,6 +87,10 @@ "__version__", "assert_count_rate", "assert_coverage", + "assert_intervals_informative", + "assert_more_powerful", + "assert_narrower", + "assert_power", "assert_proportion", "assert_se_calibrated", "assert_unbiased", @@ -81,4 +98,7 @@ "deep_tier", "monte_carlo", "reps_for", + "se_ratio_tolerance", + "vacuous_width_ratio", + "width_ratio", ] diff --git a/src/simcheck/gates.py b/src/simcheck/gates.py index 79c7ba7..007a29b 100644 --- a/src/simcheck/gates.py +++ b/src/simcheck/gates.py @@ -23,10 +23,24 @@ whole product is assertions would silently pass everything under optimisation -- the exact failure mode it exists to prevent, in itself. ``test_negative.py`` runs the gates in an ``-O`` subprocess to keep it that way. + +**Where each threshold comes from is written down.** Three shapes appear here: + +* A rate against a claimed rate -- coverage, size, power. The band is binomial in + ``reps`` and there is nothing to choose: :func:`binomial_band`. +* A quantity against its own Monte Carlo standard error -- bias, the + claimed-to-actual standard error ratio, the gap between two methods. The + tolerance is again ``reps``, through the sampling distribution of the quantity: + :func:`se_ratio_tolerance`. +* An interval's *width*. Here the reference width is derived from the study + (:func:`vacuous_width_ratio`), but the multiple of it at which conservatism + becomes vacuity is a judgement, and that function's docstring says so plainly + rather than presenting a chosen number as a derived one. """ from __future__ import annotations +import math from typing import TYPE_CHECKING import numpy as np @@ -38,10 +52,17 @@ "GATE_SIGMAS", "assert_count_rate", "assert_coverage", + "assert_intervals_informative", + "assert_more_powerful", + "assert_narrower", + "assert_power", "assert_proportion", "assert_se_calibrated", "assert_unbiased", "binomial_band", + "se_ratio_tolerance", + "vacuous_width_ratio", + "width_ratio", ] # How many standard errors a result may sit from nominal before it counts as @@ -114,6 +135,175 @@ def _band_failure( ) +def _normal_tail_quantile(tail: float) -> float: + """The ``z`` with ``P(Z > z) == tail`` for a standard normal ``Z``. + + numpy has no normal quantile function and simcheck depends on numpy alone, so + this inverts ``math.erfc`` by bisection rather than acquiring scipy for one + call. It works in the upper tail rather than in the cumulative probability + because the tail is where it is used: ``1 - Phi(z)`` has lost every + significant digit by ``z = 9``, while ``erfc`` has not. + + Args: + tail: Upper-tail probability, in ``(0, 0.5]``. + + Returns: + float: The quantile. + + Raises: + ValueError: If ``tail`` is outside ``(0, 0.5]``. + """ + if not 0.0 < tail <= 0.5: + raise ValueError(f"tail probability must be in (0, 0.5], got {tail}") + low, high = 0.0, 40.0 + for _ in range(120): + mid = 0.5 * (low + high) + if 0.5 * math.erfc(mid / math.sqrt(2.0)) > tail: + low = mid + else: + high = mid + return 0.5 * (low + high) + + +def vacuous_width_ratio(nominal: float, reps: int) -> float: + """How many times the calibrated width an interval may reach before it is vacuous. + + **The reference width is derived, and contains no chosen number.** For an + estimator whose sampling distribution is approximately normal with spread + ``sigma``, the shortest interval that contains the truth at rate ``1 - alpha`` + has width ``2 * z_{1 - alpha/2} * sigma``. The study measures ``sigma`` + itself, as ``sampling_sd``, so the width an interval *should* have is a + measurement rather than a threshold. :func:`width_ratio` reports the observed + mean width in units of it. + + **The multiple returned here is derived from ``reps``, with one convention in + it that is stated rather than hidden.** Widen a calibrated interval by a + factor ``r`` and its miss rate falls to ``q(r) = 2 * (1 - Phi(r * z))``. This + function returns the ``r`` at which the *whole study* expects fewer than + ``alpha`` misses -- that is, at which observing a single failure would take + ``1 / alpha`` studies of this size. Past that point the coverage a study + reports is a property of the width rather than of the estimator: the interval + could not have failed, so its covering says nothing. Solving + ``reps * q(r) = alpha`` gives ``r = z_{1 - alpha/(2*reps)} / z_{1 - alpha/2}``. + + The threshold therefore *loosens* as the study grows -- 1.78 at 100 + replicates, 1.96 at 400, 2.15 at 2000 for a nominal 0.95 -- which inverts the + usual direction and is meant to: more replicates resolve rarer failures, so + an interval must be wider before a study of that size can no longer see it + fail. + + **What is not derived.** How many expected misses per study counts as "could + not have failed" -- ``alpha`` of one, here, taken from the interval's own + claim rather than invented -- is a convention. No sampling distribution fixes + it, because correct procedures occupy the whole range above one: a Student t + interval at ``n = 5`` is 1.33 times the normal oracle width and an + anytime-valid interval is wider still, both of them right. That is why + :func:`assert_intervals_informative` does not fail on width alone, but only + on width *together with* a study that never once saw the interval miss. A + procedure that fails at its nominal rate is not vacuous however wide it is, + and no threshold on width can be asked to know that. + + Args: + nominal: The level the intervals claim, for instance 0.95. + reps: Number of replicates in the study. + + Returns: + float: The width multiple at which the study loses the ability to + observe the interval failing. + + Raises: + ValueError: If ``nominal`` is not strictly inside ``(0, 1)`` or ``reps`` + is not positive. + + Examples: + >>> round(vacuous_width_ratio(0.95, 400), 3) + 1.957 + >>> round(vacuous_width_ratio(0.95, 100), 3) + 1.776 + """ + if not 0.0 < nominal < 1.0: + raise ValueError(f"nominal must be strictly inside (0, 1), got {nominal}") + if reps <= 0: + raise ValueError(f"reps must be positive, got {reps}") + alpha = 1.0 - nominal + return _normal_tail_quantile(alpha / (2.0 * reps)) / _normal_tail_quantile( + alpha / 2.0 + ) + + +def width_ratio(result: MonteCarloResult, nominal: float = 0.95) -> float: + """Mean interval width over the width a calibrated interval would have. + + One means the interval is as narrow as its level allows against the spread + this estimator actually has; two means it is twice as wide as it needs to be. + The denominator, ``2 * z_{1 - alpha/2} * sampling_sd``, comes from the study, + so no absolute width is written down anywhere. + + Args: + result: A completed Monte Carlo study that recorded interval endpoints. + nominal: The level the intervals claim. + + Returns: + float: The ratio. + + Raises: + ValueError: If the study recorded no interval endpoints, if ``nominal`` + is not strictly inside ``(0, 1)``, or if the estimator did not vary + across replicates, which leaves no spread to measure the width + against. + """ + if not 0.0 < nominal < 1.0: + raise ValueError(f"nominal must be strictly inside (0, 1), got {nominal}") + spread = result.sampling_sd + if not spread: + raise ValueError( + f"the estimator did not vary across {result.reps} replicates, so " + "there is no sampling spread to compare its interval width against" + ) + calibrated = 2.0 * _normal_tail_quantile((1.0 - nominal) / 2.0) * spread + return result.mean_width / calibrated + + +def _mean_gap_se(first: np.ndarray, second: np.ndarray) -> float: + """Monte Carlo standard error of the difference between two means. + + The two studies are treated as independent, which is what they are when they + were run under different seeds. Running them under the *same* seed pairs the + replicates and makes the true standard error of the difference smaller, so + this is the conservative choice: a paired comparison passing this gate would + also pass a paired one. + + Args: + first: Per-replicate values from one study. + second: Per-replicate values from the other. + + Returns: + float: The standard error of ``mean(first) - mean(second)``. + """ + if len(first) < 2 or len(second) < 2: + return 0.0 + return math.sqrt( + float(np.var(first, ddof=1)) / len(first) + + float(np.var(second, ddof=1)) / len(second) + ) + + +def _gap_is_unresolved(gap: float, standard_error: float, sigmas: float) -> bool: + """Whether a difference is too small for the studies to resolve. + + Args: + gap: The observed difference, signed so that positive is the claim. + standard_error: Monte Carlo standard error of the difference. + sigmas: How many of them the gap must exceed. + + Returns: + bool: True when the claim is not established. + """ + if standard_error == 0.0: + return gap <= 0.0 + return gap <= sigmas * standard_error + + def assert_proportion( observed: float, reps: int, @@ -228,10 +418,64 @@ def assert_coverage( raise AssertionError(f"{f'{label} coverage'.strip()}: {problem}") +def se_ratio_tolerance(result: MonteCarloResult, sigmas: float = GATE_SIGMAS) -> float: + """How far ``se_ratio`` can sit from one on Monte Carlo noise alone. + + ``se_ratio`` is ``mean(reported standard errors) / sd(estimates)``, and both + halves are estimated from the same ``reps`` replicates, so it is noisy even + when the estimator is perfect. Its sampling distribution is available: + + * The numerator is a mean of ``reps`` reported standard errors, so its + relative standard error is ``cv / sqrt(reps)``, where ``cv`` is their + coefficient of variation across replicates. An estimator that reports the + same standard error every time contributes nothing here. + * The denominator is a sample standard deviation of ``reps`` draws, whose + relative standard error is ``1 / sqrt(2 * (reps - 1))`` -- exactly for + normal estimates and closely for anything with a finite fourth moment. + + Adding them in quadrature and multiplying by ``sigmas`` gives the band. The + two are in fact positively correlated for most estimators -- a replicate that + produces a large estimate often reports a large standard error too -- and + ignoring that overstates the variance, which makes this the lenient choice. + + At 100 replicates the band is about 0.21 and at 2000 about 0.05, so the same + call is a sanity check in a fast tier and a real test in a deep one. That is + the whole point: 0.15 was the one number in this package chosen by hand + rather than derived, and it was simultaneously too loose to catch a 12% + error in a 2000-replicate study and tight enough to fail a correct estimator + roughly one time in ten at 50. + + Args: + result: A completed Monte Carlo study. + sigmas: How many Monte Carlo standard errors of slack to allow. + + Returns: + float: The largest deviation of ``se_ratio`` from one that this study + cannot distinguish from noise. + + Raises: + ValueError: If the study has fewer than two replicates, which leaves the + spread -- and so the ratio -- undefined. + """ + if result.reps < 2: + raise ValueError( + "a single replicate has no spread, so its reported standard error " + "cannot be checked against anything" + ) + reported = np.asarray(result.standard_errors, dtype=float) + mean_reported = float(np.mean(reported)) + variation = ( + float(np.std(reported, ddof=1)) / mean_reported if mean_reported else 0.0 + ) + relative = math.sqrt(variation**2 / result.reps + 1.0 / (2.0 * (result.reps - 1))) + return sigmas * relative + + def assert_se_calibrated( result: MonteCarloResult, label: str = "", - tolerance: float = 0.15, + tolerance: float | None = None, + sigmas: float = GATE_SIGMAS, ) -> None: """Fail if the reported standard error misstates the estimator's spread. @@ -239,29 +483,305 @@ def assert_se_calibrated( errors cancel -- an inflated standard error paired with a bias, say. This checks the standard error directly against the spread actually observed. + The tolerance comes from the replicate count by default, through the sampling + distribution of the ratio: see :func:`se_ratio_tolerance`. Passing a number + overrides it, which is worth doing only when the claim being tested is about + a fixed accuracy -- "this sandwich estimator is within 5% at this sample + size" -- rather than about the standard error being right. + Args: result: A completed Monte Carlo study. label: Included in the failure message. tolerance: Largest permitted relative deviation of ``se_ratio`` from one. + Derived from ``reps`` when omitted. + sigmas: How many Monte Carlo standard errors of slack the derived + tolerance allows. Ignored when ``tolerance`` is given. Raises: ValueError: If ``tolerance`` is not positive. - AssertionError: If the ratio falls outside ``1 +- tolerance``, or the - estimator did not vary at all across replicates. + AssertionError: If the ratio falls outside ``1 +- tolerance``, if the + estimator did not vary at all across replicates, or if it reported no + standard error to check. """ - if tolerance <= 0: + if tolerance is not None and tolerance <= 0: raise ValueError(f"tolerance must be positive, got {tolerance}") - ratio = result.se_ratio - if not np.isfinite(ratio): + if not result.sampling_sd: raise AssertionError( f"{label}: the estimator did not vary across {result.reps} " "replicates, so its reported standard error cannot be checked " "against anything" ) - if abs(ratio - 1.0) > tolerance: + ratio = result.se_ratio + if not np.isfinite(ratio): + raise AssertionError( + f"{label}: the estimator reported no usable standard error over " + f"{result.reps} replicates (mean of the reported values is " + f"{result.reported_se}), so there is nothing to check against its " + f"observed spread of {result.sampling_sd:.6f}" + ) + derived = tolerance is None + band = se_ratio_tolerance(result, sigmas) if derived else float(tolerance) + if abs(ratio - 1.0) > band: + source = ( + f"{sigmas:g} Monte Carlo standard errors of the ratio at " + f"{result.reps} replicates" + if derived + else "supplied by the caller" + ) raise AssertionError( f"{label}: reported standard error is {ratio:.3f} times the " f"observed spread ({result.reported_se:.6f} against " f"{result.sampling_sd:.6f}) over {result.reps} replicates; " - f"tolerance is {tolerance:.2f}" + f"tolerance is {band:.3f}, {source}" + ) + + +def assert_intervals_informative( + result: MonteCarloResult, + nominal: float = 0.95, + label: str = "", + max_ratio: float | None = None, + sigmas: float = GATE_SIGMAS, +) -> None: + """Fail if the intervals are so wide that their coverage means nothing. + + :func:`assert_coverage` is satisfied by an interval that always covers, + whenever the study is small enough that a rate of 1.0 still sits inside the + binomial band -- and it is *always* satisfied by ``coverage > 0.9`` written + by hand. Three separate repositories worked around this with a comment + saying so; one of them had shipped an inflation heuristic that drove the + reported standard error to 3e7 times the estimation error while coverage + stayed high, because a vacuous interval covers everything. + + Two things must both be true before this fails, and the conjunction is the + point: + + 1. **The interval is far wider than it needs to be.** ``width_ratio``, the + mean width over the width a calibrated interval would have against this + estimator's own spread, exceeds ``max_ratio`` by more than Monte Carlo + noise. The default ``max_ratio`` is :func:`vacuous_width_ratio`, derived + from ``nominal`` and ``reps``. + 2. **The study never once saw the interval fail.** Fewer than ``sigmas`` + misses in ``reps`` replicates: by the rule of three, a study observing no + failures bounds the miss rate only at ``sigmas / reps``, so its coverage + number is censored rather than measured. + + Requiring both is what keeps the gate off correct code. A Student t interval + at ``n = 5`` is 1.33 times the normal oracle width, and an anytime-valid + interval more, but both miss at their nominal rate, which the study sees, so + neither is vacuous. Width alone cannot tell conservatism from vacuity; + width plus a study that never saw a failure can. + + Args: + result: A completed Monte Carlo study that recorded interval endpoints. + nominal: The level the intervals claim. + label: Included in the failure message. + max_ratio: Override for the derived width multiple. + sigmas: Monte Carlo slack on the width ratio, and the miss count below + which the study is treated as never having seen a failure. + + Raises: + ValueError: If the study recorded no interval endpoints, so there is no + width to check, or if ``max_ratio`` is not positive. + AssertionError: If the intervals are vacuous, or if the estimator did not + vary at all across replicates, which leaves nothing to compare their + width against. + """ + lowers, uppers = result.lowers, result.uppers + if lowers is None or uppers is None: + raise ValueError( + f"{label or 'this study'} recorded no interval endpoints, so the " + "width of its intervals cannot be checked. Have the estimator " + "report `lower` and `upper` on every replicate." ) + if max_ratio is not None and max_ratio <= 0: + raise ValueError(f"max_ratio must be positive, got {max_ratio}") + if not result.sampling_sd: + raise AssertionError( + f"{label}: the estimator did not vary across {result.reps} " + "replicates, so there is no spread to compare its interval width " + "against" + ) + + ratio = width_ratio(result, nominal) + threshold = ( + vacuous_width_ratio(nominal, result.reps) if max_ratio is None else max_ratio + ) + widths = result.widths + mean_width = result.mean_width + relative = math.sqrt( + (float(np.var(widths, ddof=1)) / result.reps) / mean_width**2 + + 1.0 / (2.0 * (result.reps - 1)) + if mean_width + else 0.0 + ) + if ratio - sigmas * ratio * relative <= threshold: + return + + misses = int(np.count_nonzero((lowers > result.truth) | (result.truth > uppers))) + if misses >= sigmas: + return + calibrated = 2.0 * _normal_tail_quantile((1.0 - nominal) / 2.0) * result.sampling_sd + raise AssertionError( + f"{label}: intervals are vacuous. Mean width {mean_width:.6g} is " + f"{ratio:.3g} times the {calibrated:.6g} " + f"a calibrated {nominal:.2f} interval needs against this estimator's " + f"spread ({result.sampling_sd:.6g}), and the study saw {misses} misses " + f"in {result.reps} replicates, so its coverage of " + f"{1 - misses / result.reps:.3f} " + f"measures the width rather than the estimator. The width at which a " + f"study of this size stops being able to see a miss is " + f"{threshold:.3g} times calibrated." + ) + + +def assert_narrower( + narrow: MonteCarloResult, + wide: MonteCarloResult, + label: str = "", + sigmas: float = GATE_SIGMAS, +) -> None: + """Fail unless one method's intervals are measurably narrower than another's. + + The efficiency half of an interval comparison. Width without coverage is not + a virtue -- the narrowest interval of all is the empty one -- so this is + meant to be run *after* :func:`assert_coverage` on both studies, and it says + nothing about either one's calibration. + + The tolerance is the Monte Carlo standard error of the difference in mean + width over the two studies, so a gap this study cannot resolve does not pass, + and the same call becomes stricter as the studies grow. + + Args: + narrow: The study claimed to produce the narrower intervals. + wide: The study it is claimed to beat. + label: Included in the failure message. + sigmas: How many Monte Carlo standard errors the gap must exceed. + + Raises: + ValueError: If either study recorded no interval endpoints. + AssertionError: If the narrower study's intervals are not measurably + narrower. + """ + for name, study in (("narrow", narrow), ("wide", wide)): + if study.lowers is None or study.uppers is None: + raise ValueError( + f"{label or 'this comparison'}: the `{name}` study recorded no " + "interval endpoints, so the two widths cannot be compared" + ) + gap = wide.mean_width - narrow.mean_width + standard_error = _mean_gap_se(wide.widths, narrow.widths) + if not _gap_is_unresolved(gap, standard_error, sigmas): + return + resolved = ( + f"{gap / standard_error:+.2f} Monte Carlo standard errors" + if standard_error + else "an exactly zero Monte Carlo standard error" + ) + raise AssertionError( + f"{label}: mean width {narrow.mean_width:.6g} over {narrow.reps} " + f"replicates is not measurably below {wide.mean_width:.6g} over " + f"{wide.reps}: the gap of {gap:+.6g} is {resolved}, against a gate of " + f"{sigmas:g}" + ) + + +def assert_power( + result: MonteCarloResult, + minimum: float, + label: str = "", + sigmas: float = GATE_SIGMAS, +) -> None: + """Fail if a test rejects less often than claimed under an alternative. + + One-sided, unlike :func:`assert_proportion`: power is a floor, not a target. + Rejecting *more* often than the claim is not a defect of the test, and a + two-sided band would fail an estimator for being better than promised. Size, + which is a target, still belongs in :func:`assert_proportion`. + + The floor is the lower end of the binomial band around ``minimum``, so the + claim is "power is at least ``minimum``, and this study can say so" rather + than "the observed rate happened to clear ``minimum``". Under a claim that is + exactly true the gate fires about once in 740 studies. + + Args: + result: A completed Monte Carlo study run under the alternative. + minimum: The power being claimed, in ``[0, 1]``. + label: Included in the failure message. + sigmas: Slack, in binomial standard errors. + + Raises: + ValueError: If ``minimum`` is not in ``[0, 1]``, or if the study recorded + no reject/accept decisions and so has no power to check. + AssertionError: If the rejection rate falls below the floor. + """ + if not 0.0 <= minimum <= 1.0: + raise ValueError(f"minimum must be a probability, got {minimum}") + if result.rejected is None: + raise ValueError( + f"{label or 'this study'} recorded no reject/accept decisions, so " + "its power cannot be checked. Have the estimator report `rejected` " + "on every replicate." + ) + floor, _ = binomial_band(minimum, result.reps, sigmas) + observed = result.rejection_rate + if observed >= floor: + return + raise AssertionError( + f"{label}: power {observed:.4f} is below the one-sided {sigmas:g}-sigma " + f"floor {floor:.4f} for a claimed minimum of {minimum:.4f} over " + f"{result.reps} replicates" + ) + + +def assert_more_powerful( + more: MonteCarloResult, + less: MonteCarloResult, + label: str = "", + sigmas: float = GATE_SIGMAS, +) -> None: + """Fail unless one test rejects measurably more often than another. + + The comparison a method paper actually makes: at the same alternative, does A + detect it more often than B. Asserting ``a.rejection_rate > b.rejection_rate`` + instead -- which is what this replaces -- passes on a gap of one replicate in + four hundred, which is noise, and so certifies whichever method the seed + happened to favour. + + Both studies must be run at the same alternative, which this cannot check. + Comparing rejection rates under *different* alternatives compares the + alternatives, not the tests. + + Args: + more: The study claimed to be more powerful. + less: The study it is claimed to beat. + label: Included in the failure message. + sigmas: How many standard errors of the difference the gap must exceed. + + Raises: + ValueError: If either study recorded no reject/accept decisions. + AssertionError: If the gap is not measurably positive. + """ + for name, study in (("more", more), ("less", less)): + if study.rejected is None: + raise ValueError( + f"{label or 'this comparison'}: the `{name}` study recorded no " + "reject/accept decisions, so the two cannot be compared" + ) + strong, weak = more.rejection_rate, less.rejection_rate + gap = strong - weak + standard_error = math.sqrt( + strong * (1.0 - strong) / more.reps + weak * (1.0 - weak) / less.reps + ) + if not _gap_is_unresolved(gap, standard_error, sigmas): + return + resolved = ( + f"{gap / standard_error:+.2f} standard errors of the difference" + if standard_error + else "an exactly zero standard error" + ) + raise AssertionError( + f"{label}: rejection rate {strong:.4f} over {more.reps} replicates is " + f"not measurably above {weak:.4f} over {less.reps}: the gap of " + f"{gap:+.4f} is {resolved}, against a gate of {sigmas:g}" + ) diff --git a/src/simcheck/results.py b/src/simcheck/results.py index 03e1216..fb7a369 100644 --- a/src/simcheck/results.py +++ b/src/simcheck/results.py @@ -12,6 +12,12 @@ size while adding almost no information, and makes a badly calibrated estimator look precisely measured. Keep one binomial per point and the arithmetic stays honest. + +**The endpoints are kept, not just the hit.** Coverage on its own cannot tell a +calibrated interval from one so wide it could not have failed, and reducing each +replicate to a boolean throws away the only evidence that would. Keeping +``lowers`` and ``uppers`` is what lets :func:`~simcheck.assert_intervals_informative` +measure the width against the spread the estimator actually has. """ from __future__ import annotations @@ -38,9 +44,17 @@ class MonteCarloResult: rejected: Whether each replicate rejected the null, or None when the estimator performs no test. truth: The true value of the quantity being estimated. + lowers: Lower endpoint of each replicate's interval, or None when the + estimator reported no intervals. Supplying the endpoints is what + makes the interval's *width* checkable, and coverage alone cannot + distinguish a calibrated interval from a vacuous one. + uppers: Upper endpoint of each replicate's interval. Both endpoints or + neither; one without the other is not an interval. Raises: - ValueError: If the arrays disagree in length or are empty. + ValueError: If the arrays disagree in length or are empty, if one + endpoint array is given without the other, if any interval runs + backwards, or if ``covered`` contradicts the endpoints. """ estimates: npt.NDArray[np.float64] @@ -48,18 +62,35 @@ class MonteCarloResult: covered: npt.NDArray[np.bool_] | None rejected: npt.NDArray[np.bool_] | None truth: float + lowers: npt.NDArray[np.float64] | None = None + uppers: npt.NDArray[np.float64] | None = None def __post_init__(self) -> None: - """Reject ragged or empty input. + """Reject ragged, empty or self-contradictory input. + + When the endpoints are given and ``covered`` is not, ``covered`` is + filled in from them. That is not a guess: the endpoints and the truth + determine coverage exactly, and leaving the field None would make a + study that plainly measured coverage report that it had not. Raises: - ValueError: If the arrays disagree in length or are empty. + ValueError: If the arrays disagree in length or are empty, if one + endpoint array is given without the other, if any interval runs + backwards, or if ``covered`` contradicts the endpoints. """ + if (self.lowers is None) != (self.uppers is None): + raise ValueError( + "an interval needs both endpoints: got " + f"lowers={'None' if self.lowers is None else 'an array'} and " + f"uppers={'None' if self.uppers is None else 'an array'}" + ) lengths = { len(self.estimates), len(self.standard_errors), *([] if self.covered is None else [len(self.covered)]), *([] if self.rejected is None else [len(self.rejected)]), + *([] if self.lowers is None else [len(self.lowers)]), + *([] if self.uppers is None else [len(self.uppers)]), } if len(lengths) > 1: raise ValueError( @@ -67,6 +98,28 @@ def __post_init__(self) -> None: ) if lengths == {0}: raise ValueError("a Monte Carlo study needs at least one replicate") + if self.lowers is None or self.uppers is None: + return + + backwards = int(np.count_nonzero(self.uppers < self.lowers)) + if backwards: + raise ValueError( + f"{backwards} of {len(self.lowers)} intervals have an upper " + "endpoint below the lower one, so their width is negative" + ) + hit = (self.lowers <= self.truth) & (self.truth <= self.uppers) + if self.covered is None: + object.__setattr__(self, "covered", hit) + return + disagreements = int(np.count_nonzero(np.asarray(self.covered) != hit)) + if disagreements: + raise ValueError( + f"`covered` disagrees with the endpoints on {disagreements} of " + f"{len(hit)} replicates. One of the two is measuring something " + "the other is not -- an interval on a transformed scale, say -- " + "and pooling them would report a coverage rate for one and a " + "width for the other." + ) @property def reps(self) -> int: @@ -141,6 +194,52 @@ def coverage(self) -> float: ) return float(np.mean(self.covered)) + @property + def widths(self) -> npt.NDArray[np.float64]: + """Width of each replicate's interval. + + Returns: + numpy.ndarray: One width per replicate. + + Raises: + ValueError: If the study did not record interval endpoints. Reporting + zeros here would read as an infinitely precise estimator, which + is the opposite of an unmeasured one. + """ + if self.lowers is None or self.uppers is None: + raise ValueError( + "this study recorded no interval endpoints, so it has no " + "widths. Have the estimator report `lower` and `upper` on every " + "replicate." + ) + return np.asarray(self.uppers - self.lowers, dtype=float) + + @property + def mean_width(self) -> float: + """Mean interval width, in the units of the estimand. + + Reading this on a study that recorded no endpoints raises ``ValueError`` + through :attr:`widths`, rather than reporting a width of zero. + + Returns: + float: The mean width. + """ + return float(np.mean(self.widths)) + + @property + def median_width(self) -> float: + """Median interval width. + + Worth having alongside the mean: a procedure that returns an enormous + interval on a few replicates has a mean width dominated by those, and the + median says what the typical replicate produced. Reading it on a study + that recorded no endpoints raises ``ValueError`` through :attr:`widths`. + + Returns: + float: The median width. + """ + return float(np.median(self.widths)) + @property def rejection_rate(self) -> float: """Fraction of replicates that rejected the null. diff --git a/src/simcheck/runner.py b/src/simcheck/runner.py index 3e238d0..35b45a9 100644 --- a/src/simcheck/runner.py +++ b/src/simcheck/runner.py @@ -81,7 +81,10 @@ def monte_carlo( ``reps`` adds replicates rather than changing the existing ones. Returns: - MonteCarloResult: The recorded sampling behaviour. + MonteCarloResult: The recorded sampling behaviour. The interval + endpoints are kept, not only whether each interval covered, because + coverage on its own cannot tell a calibrated interval from one so wide + it could not have failed. Raises: ValueError: If ``reps`` is not positive, or if the estimator reported an @@ -127,8 +130,10 @@ def monte_carlo( flags[i] = bool(outcome.rejected) covered = None + endpoints: tuple[np.ndarray, np.ndarray] | tuple[None, None] = (None, None) if with_interval.all(): covered = (lowers <= truth) & (truth <= uppers) + endpoints = (lowers, uppers) elif with_interval.any(): raise ValueError( f"{int(with_interval.sum())} of {reps} replicates reported an " @@ -147,4 +152,6 @@ def monte_carlo( "a mixture is not a size or a power." ) - return MonteCarloResult(values, errors, covered, rejected, truth) + return MonteCarloResult( + values, errors, covered, rejected, truth, endpoints[0], endpoints[1] + ) diff --git a/tests/test_against_known_statistics.py b/tests/test_against_known_statistics.py index 601a7fa..dd306f5 100644 --- a/tests/test_against_known_statistics.py +++ b/tests/test_against_known_statistics.py @@ -14,21 +14,90 @@ from __future__ import annotations +import math + import numpy as np import pytest from simcheck import ( + Estimate, MonteCarloResult, assert_coverage, + assert_intervals_informative, + assert_more_powerful, + assert_narrower, + assert_power, + assert_proportion, assert_se_calibrated, assert_unbiased, binomial_band, + monte_carlo, + width_ratio, ) REPS = 2000 N = 25 TRUTH = 3.0 SIGMA = 2.0 +# The two-sided 95% normal quantile, and the t quantile on 4 degrees of freedom, +# from tables. simcheck has no scipy dependency and should not gain one so a test +# can call `norm.ppf`. +Z_95 = 1.959964 +T_95_DF4 = 2.776445 + + +def _normal_cdf(x: float) -> float: + """Standard normal cumulative distribution function. + + Args: + x: Where to evaluate it. + + Returns: + float: ``P(Z <= x)``. + """ + return 0.5 * math.erfc(-x / math.sqrt(2.0)) + + +def _z_test_power(delta: float, n: int, sigma: float = 1.0) -> float: + """Textbook power of a two-sided z test of a zero mean. + + Args: + delta: True mean under the alternative. + n: Observations per replicate. + sigma: Known standard deviation. + + Returns: + float: The probability of rejecting. + """ + shift = delta * math.sqrt(n) / sigma + return _normal_cdf(-Z_95 + shift) + _normal_cdf(-Z_95 - shift) + + +def _z_test_study(delta: float, n: int, reps: int = REPS, seed: int = 0): + """Monte Carlo study of a two-sided z test with a known scale. + + Args: + delta: True mean under the alternative. + n: Observations per replicate. + reps: Replicates. + seed: Seed for the replicate stream. + + Returns: + MonteCarloResult: The completed study. + """ + error = 1.0 / math.sqrt(n) + + def replicate(rng: np.random.Generator) -> Estimate: + mean = float(rng.normal(delta, 1.0, n).mean()) + return Estimate( + mean, + error, + mean - Z_95 * error, + mean + Z_95 * error, + rejected=abs(mean) > Z_95 * error, + ) + + return monte_carlo(replicate, delta, reps, seed=seed) def _study_sample_mean(reps: int = REPS, n: int = N, seed: int = 0): @@ -159,6 +228,100 @@ def test_the_gate_false_positive_rate_is_about_what_three_sigma_implies(): assert tripped <= 1, f"{tripped} of 60 correct studies were flagged as biased" +def test_the_power_of_a_z_test_is_the_one_the_formula_gives(): + """Power has a closed form, so the gate can be checked against it. + + For a two-sided z test at the 5% level, ``Phi(-1.96 + delta*sqrt(n)/sigma) + + Phi(-1.96 - delta*sqrt(n)/sigma)``. At delta=0.3 and n=25 that is 0.323. A + gate for power that could not recover a number this well known would be + measuring something else. + """ + expected = _z_test_power(0.3, 25) + assert expected == pytest.approx(0.3230, abs=1e-4) + + study = _z_test_study(0.3, 25, seed=31) + assert_proportion(study.rejection_rate, study.reps, expected, "z test at n=25") + assert_power(study, expected, "z test at n=25") + + # And it is genuinely underpowered against a claim of one half. + with pytest.raises(AssertionError, match="below the one-sided"): + assert_power(study, 0.50, "z test at n=25") + + +def test_a_sixteen_fold_sample_is_caught_as_more_powerful(): + """Power rises with n by a known amount, and the paired gate must see it. + + From 0.323 at n=25 to 0.851 at n=100, both from the formula. Asserting + ``strong > weak`` would also pass here -- and would equally have passed on a + gap of one replicate, which is the failure this replaces. + """ + assert _z_test_power(0.3, 100) == pytest.approx(0.8508, abs=1e-4) + + strong = _z_test_study(0.3, 100, seed=32) + weak = _z_test_study(0.3, 25, seed=33) + assert_more_powerful(strong, weak, "n=100 against n=25") + + # Two studies of the same test differ only by noise, and must not be ranked. + with pytest.raises(AssertionError, match="not measurably above"): + assert_more_powerful( + _z_test_study(0.3, 25, seed=34), _z_test_study(0.3, 25, seed=35), "same" + ) + + +def test_an_interval_at_the_known_scale_is_exactly_as_wide_as_it_must_be(): + """With sigma known the z interval is the oracle, so its width ratio is one. + + This pins the denominator of :func:`width_ratio`: if the reference width were + off by a factor, this number would not come out at one. + """ + study = _z_test_study(0.0, 25, seed=36) + assert width_ratio(study, 0.95) == pytest.approx(1.0, abs=0.05) + assert_coverage(study, 0.95, "z interval at the known scale") + assert_intervals_informative(study, 0.95, "z interval at the known scale") + + +def test_the_normal_interval_is_narrower_than_the_t_interval_and_that_is_the_defect(): + """Narrower is not better, and the textbook case says why. + + At n=5 the 1.96 interval is ``1.96 / 2.776 = 0.706`` times the width of the t + interval, which :func:`assert_narrower` confirms -- and it covers at 0.875 + rather than 0.95, which :func:`assert_coverage` catches. An efficiency + comparison run without a coverage gate on both sides would report the broken + interval as the better one. + """ + rng = np.random.default_rng(37) + n = 5 + draws = rng.normal(TRUTH, SIGMA, size=(REPS, n)) + means = draws.mean(axis=1) + errors = draws.std(axis=1, ddof=1) / np.sqrt(n) + + def study(critical: float) -> MonteCarloResult: + return MonteCarloResult( + estimates=means, + standard_errors=errors, + covered=None, + rejected=None, + truth=TRUTH, + lowers=means - critical * errors, + uppers=means + critical * errors, + ) + + normal, student = study(Z_95), study(T_95_DF4) + assert normal.mean_width / student.mean_width == pytest.approx( + Z_95 / T_95_DF4, rel=1e-12 + ) + assert_narrower(normal, student, "1.96 against t at n=5") + assert_coverage(student, 0.95, "t interval at n=5") + with pytest.raises(AssertionError, match="outside the 3-sigma band"): + assert_coverage(normal, 0.95, "1.96 interval at n=5") + + # The t interval is 1.33 times the oracle width and still not vacuous: the + # study watches it miss at its nominal rate, which is what tells the two + # apart. + assert width_ratio(student, 0.95) == pytest.approx(1.33, rel=0.1) + assert_intervals_informative(student, 0.95, "t interval at n=5") + + def test_the_band_matches_the_textbook_binomial_interval(): """The band is nominal +- sigmas * sqrt(p(1-p)/n), not something invented.""" nominal, reps, sigmas = 0.95, 400, 3.0 diff --git a/tests/test_negative.py b/tests/test_negative.py index d91da66..0cedb4e 100644 --- a/tests/test_negative.py +++ b/tests/test_negative.py @@ -21,6 +21,7 @@ from __future__ import annotations +import math import subprocess import sys import textwrap @@ -33,12 +34,23 @@ MonteCarloResult, assert_count_rate, assert_coverage, + assert_intervals_informative, + assert_more_powerful, + assert_narrower, + assert_power, assert_proportion, assert_se_calibrated, assert_unbiased, binomial_band, + se_ratio_tolerance, + vacuous_width_ratio, + width_ratio, ) +# Two-sided normal quantile for 95%, from tables. simcheck has no scipy +# dependency and should not gain one so a test can call `norm.ppf`. +Z_95 = 1.959964 + def _study( reps: int = 400, @@ -87,6 +99,85 @@ def _study( return MonteCarloResult(estimates, errors, covered, rejected, truth) +def _interval_study( + reps: int = 400, + sd: float = 0.1, + truth: float = 1.0, + bias: float = 0.0, + half_width: float | None = None, + width_cv: float = 0.0, + seed: int = 0, +) -> MonteCarloResult: + """Build a study whose intervals have a known width and a known coverage. + + Parameters + ---------- + reps + Replicates. + sd + Spread of the estimates. + truth + True value. + bias + Constant added to every estimate. + half_width + Half-width of every interval. Defaults to the calibrated 95% one, + ``1.96 * sd``, which makes ``width_ratio`` exactly one. + width_cv + Relative spread of the width across replicates. Zero gives every + replicate the same width, as a procedure with a known scale would; a + real procedure estimates the scale and so has a width that varies. + seed + Generator seed. + + Returns + ------- + MonteCarloResult + The constructed study. ``covered`` is left for the result object to + derive from the endpoints, which is also what exercises that path. + """ + rng = np.random.default_rng(seed) + estimates = rng.normal(truth + bias, sd, reps) + half = Z_95 * sd if half_width is None else half_width + if width_cv: + half = half * np.abs(1.0 + width_cv * rng.standard_normal(reps)) + return MonteCarloResult( + estimates=estimates, + standard_errors=np.full(reps, sd), + covered=None, + rejected=None, + truth=truth, + lowers=estimates - half, + uppers=estimates + half, + ) + + +def _decision_study(rejection: float, reps: int = 400) -> MonteCarloResult: + """A study with an exact rejection rate and nothing else of interest. + + Parameters + ---------- + rejection + Exact fraction of replicates to mark rejected. + reps + Replicates. + + Returns + ------- + MonteCarloResult + The constructed study. + """ + rejected = np.zeros(reps, dtype=bool) + rejected[: round(rejection * reps)] = True + return MonteCarloResult( + estimates=np.linspace(0.9, 1.1, reps), + standard_errors=np.full(reps, 0.05), + covered=None, + rejected=rejected, + truth=1.0, + ) + + # -------------------------------------------------------------------------- # The bug this package exists to not repeat. # -------------------------------------------------------------------------- @@ -195,6 +286,363 @@ def test_assert_se_calibrated_refuses_a_degenerate_study(): assert_se_calibrated(constant, "constant estimator") +def test_the_se_tolerance_narrows_as_the_study_grows(): + """The derived tolerance must be a function of the replicate count. + + A constant 0.15 is simultaneously too loose to catch a 12% error in a deep + study and tight enough to fail a correct estimator in a shallow one. The + derived band is 0.21 at 100 replicates and 0.05 at 2000. + """ + tolerances = [ + se_ratio_tolerance(_study(reps=reps)) for reps in (100, 400, 2000, 20000) + ] + assert all(a > b for a, b in pairwise(tolerances)), tolerances + # Four times the replicates halves the band, since the spread of a sample + # standard deviation goes as 1/sqrt(2(reps-1)). + assert tolerances[0] / tolerances[2] == pytest.approx(4.47, rel=0.02) + + +def test_the_same_se_error_passes_a_shallow_study_and_fails_a_deep_one(): + """The point of deriving the tolerance, in one test. + + A reported standard error 12% too large sits inside the hand-picked 0.15 at + every replicate count -- a 20000-replicate study could resolve it to better + than 1% and was told not to look. Derived, the same call passes at 100 + replicates, where the study genuinely cannot tell, and fails at 2000. + """ + assert_se_calibrated(_study(reps=100, se_scale=1.12), "too small to resolve 12%") + with pytest.raises(AssertionError, match="times the observed spread"): + assert_se_calibrated(_study(reps=2000, se_scale=1.12), "deep enough") + # And the hand-picked constant would have passed both. + assert_se_calibrated(_study(reps=2000, se_scale=1.12), "0.15", tolerance=0.15) + + +def test_assert_se_calibrated_says_so_when_no_standard_error_was_reported(): + """A NaN standard error is an absent measurement, not a zero spread. + + ``Estimate`` documents that leaving ``standard_error`` as NaN means the gate + "will have nothing to check and will say so". It used to say the estimator + did not vary, which is a different -- and false -- diagnosis. + """ + study = MonteCarloResult( + estimates=np.linspace(0.9, 1.1, 100), + standard_errors=np.full(100, np.nan), + covered=None, + rejected=None, + truth=1.0, + ) + with pytest.raises(AssertionError, match="reported no usable standard error"): + assert_se_calibrated(study, "no standard error") + + +# -------------------------------------------------------------------------- +# Vacuity: an interval so wide it cannot fail. +# -------------------------------------------------------------------------- + + +def test_assert_intervals_informative_passes_a_calibrated_interval(): + """A 1.96-sigma interval is exactly as wide as its level requires.""" + study = _interval_study() + assert width_ratio(study) == pytest.approx(1.0, rel=0.05) + assert_intervals_informative(study, 0.95, "calibrated") + + +def test_assert_intervals_informative_fails_an_interval_that_cannot_miss(): + """The defect `assert_coverage` cannot see: coverage of 1.0 by construction.""" + study = _interval_study(half_width=10 * 0.1) + with pytest.raises(AssertionError, match="intervals are vacuous"): + assert_intervals_informative(study, 0.95, "ten-sigma interval") + + +def test_a_conservative_but_honest_interval_is_not_called_vacuous(): + """The false-positive guard, and the reason the gate has two conditions. + + A Student t interval at n=5 is 1.33 times the width the normal oracle needs, + which is above the naive width threshold at any replicate count -- and it is + correct. It misses at its nominal rate, the study sees those misses, so it is + conservative rather than vacuous. A gate that fired here would be switched + off within a week. + """ + reps, n = 2000, 5 + rng = np.random.default_rng(4) + draws = rng.normal(1.0, 1.0, size=(reps, n)) + means = draws.mean(axis=1) + errors = draws.std(axis=1, ddof=1) / np.sqrt(n) + half = 2.776445 * errors # t_{0.975, 4}, from tables. + study = MonteCarloResult( + estimates=means, + standard_errors=errors, + covered=None, + rejected=None, + truth=1.0, + lowers=means - half, + uppers=means + half, + ) + assert width_ratio(study) == pytest.approx(1.33, rel=0.1) + assert width_ratio(study) > vacuous_width_ratio(0.95, reps) / 2 + assert_coverage(study, 0.95, "t interval at n=5") + assert_intervals_informative(study, 0.95, "t interval at n=5") + + +def test_width_alone_does_not_make_an_interval_vacuous(): + """A wide interval the study watched fail is a different defect. + + Here the estimator is biased by five standard deviations and its intervals + are three times as wide as they need to be, so it misses about a fifth of the + time. That is a bias, which `assert_unbiased` reports; calling it vacuity as + well would send the reader after the wrong thing. + """ + study = _interval_study(bias=0.5, half_width=3 * Z_95 * 0.1) + assert width_ratio(study) == pytest.approx(3.0, rel=0.05) + assert_intervals_informative(study, 0.95, "wide, biased, and caught") + with pytest.raises(AssertionError, match="standard errors from zero"): + assert_unbiased(study, "wide, biased, and caught") + + +@pytest.mark.parametrize( + ("case", "half_width", "reps"), + [ + # geoinference: "an interval so wide it covers every single time, which + # is a defect the assertion cannot see". + ("geoinference: coverage 1.0 at 300 replicates", 8 * 0.1, 300), + # alsgls: "an interval covering 100% of the time, which means it is far + # too wide, passed it". + ("alsgls: pooled points, coverage 1.0", 6 * 0.1, 300), + # incline: "an inflation heuristic drove this ratio to 3e7 while + # coverage stayed high, because a vacuous interval covers everything". + ("incline: reported se 3e7 times the error", 3e7 * Z_95 * 0.1, 400), + ], +) +def test_the_three_workarounds_in_the_wild_are_caught(case, half_width, reps): + """Each sibling repo hand-rolled a comment where this gate should have been. + + Three repositories independently hit the same hole in `assert_coverage` and + wrote a comment about it instead of a test. If the gate cannot reproduce all + three failures it is the wrong gate, so all three run here. + """ + study = _interval_study(reps=reps, half_width=half_width) + assert study.coverage == 1.0, case + with pytest.raises(AssertionError, match="intervals are vacuous"): + assert_intervals_informative(study, 0.95, case) + + +def test_the_vacuity_gate_needs_endpoints(): + """Coverage alone cannot answer the question, and saying so beats guessing.""" + with pytest.raises(ValueError, match="recorded no interval endpoints"): + assert_intervals_informative(_study(), 0.95, "no endpoints") + + +def test_the_vacuity_threshold_comes_from_the_replicate_count(): + """More replicates resolve rarer misses, so the width limit rises with reps.""" + ratios = [vacuous_width_ratio(0.95, reps) for reps in (100, 400, 2000, 10000)] + assert all(a < b for a, b in pairwise(ratios)), ratios + # And it is the width whose implied miss rate leaves `alpha` misses expected + # in the whole study: reps * 2 * (1 - Phi(r * z)) == alpha. + reps = 400 + implied = math.erfc(ratios[1] * Z_95 / math.sqrt(2.0)) + assert reps * implied == pytest.approx(0.05, rel=1e-6) + + +# -------------------------------------------------------------------------- +# Widths across two methods. +# -------------------------------------------------------------------------- + + +def test_assert_narrower_passes_when_one_method_really_is_narrower(): + """The efficiency comparison, when the gap is real.""" + tight = _interval_study(half_width=Z_95 * 0.1, width_cv=0.2, seed=1) + loose = _interval_study(half_width=1.5 * Z_95 * 0.1, width_cv=0.2, seed=2) + assert_narrower(tight, loose, "tight against loose") + + +def test_assert_narrower_fails_when_the_gap_is_noise(): + """Two methods of the same width must not be ranked by the seed. + + ``a.mean_width < b.mean_width`` is true of one of the two whatever they are, + which is exactly the assertion this replaces. + """ + first = _interval_study(half_width=Z_95 * 0.1, width_cv=0.2, seed=1) + second = _interval_study(half_width=Z_95 * 0.1001, width_cv=0.2, seed=1) + assert first.mean_width < second.mean_width + with pytest.raises(AssertionError, match="not measurably below"): + assert_narrower(first, second, "a hair narrower") + + +def test_a_width_that_never_varies_needs_only_to_be_smaller(): + """With a known scale the width is not random, so any gap is a real gap. + + Both studies here produce exactly one width, so the Monte Carlo standard + error of the difference is exactly zero and the gate reduces to a strict + inequality. Demanding a multiple of a zero standard error would make the gate + unsatisfiable for the one case where the answer is certain. + """ + + def fixed(width: float) -> MonteCarloResult: + reps = 200 + return MonteCarloResult( + estimates=np.linspace(0.4, 0.6, reps), + standard_errors=np.full(reps, 0.05), + covered=None, + rejected=None, + truth=0.5, + lowers=np.zeros(reps), + uppers=np.full(reps, width), + ) + + assert_narrower(fixed(0.19), fixed(0.20), "known scale") + with pytest.raises(AssertionError, match="exactly zero Monte Carlo"): + assert_narrower(fixed(0.20), fixed(0.19), "known scale, backwards") + + +def test_assert_narrower_fails_when_the_arguments_are_the_wrong_way_round(): + """A wired-up-backwards comparison must not pass.""" + tight = _interval_study(half_width=Z_95 * 0.1, seed=1) + loose = _interval_study(half_width=2 * Z_95 * 0.1, seed=2) + with pytest.raises(AssertionError, match="not measurably below"): + assert_narrower(loose, tight, "backwards") + + +def test_assert_narrower_needs_endpoints_on_both_studies(): + """Comparing a width against an unmeasured one is not a comparison.""" + with pytest.raises(ValueError, match="`wide` study recorded no interval"): + assert_narrower(_interval_study(), _study(), "one side missing") + + +# -------------------------------------------------------------------------- +# Power. +# -------------------------------------------------------------------------- + + +def test_assert_power_passes_a_test_that_meets_its_claim(): + """A test rejecting 80% of the time meets a claim of 0.75.""" + assert_power(_decision_study(0.80), 0.75, "adequately powered") + + +def test_assert_power_fails_an_underpowered_test(): + """Half the claimed power over 400 replicates cannot be missed.""" + with pytest.raises(AssertionError, match="below the one-sided"): + assert_power(_decision_study(0.40), 0.80, "underpowered") + + +def test_assert_power_is_one_sided(): + """Rejecting more often than promised is not a defect of the test. + + ``assert_proportion`` bands both sides, which is right for size and wrong for + power: it would fail a test for being better than claimed. + """ + strong = _decision_study(0.95) + assert_power(strong, 0.50, "better than claimed") + with pytest.raises(AssertionError, match="outside the 3-sigma band"): + assert_proportion(strong.rejection_rate, strong.reps, 0.50, "two-sided") + + +def test_assert_power_tightens_with_replicates(): + """The floor is a binomial band, so it closes on the claim as reps grow.""" + assert_power(_decision_study(0.74, reps=100), 0.80, "cannot resolve 6 points") + with pytest.raises(AssertionError): + assert_power(_decision_study(0.74, reps=4000), 0.80, "can resolve 6 points") + + +def test_assert_power_needs_decisions(): + """A study that recorded no rejections has no power, rather than zero power.""" + with pytest.raises(ValueError, match="recorded no reject/accept decisions"): + assert_power(_interval_study(), 0.8, "no decisions") + + +def test_assert_power_rejects_an_impossible_claim(): + """A power outside [0, 1] is a caller error.""" + with pytest.raises(ValueError, match="must be a probability"): + assert_power(_decision_study(0.5), 1.4) + + +def test_assert_more_powerful_passes_on_a_real_gap(): + """Half again the rejection rate over 400 replicates is not noise.""" + assert_more_powerful(_decision_study(0.80), _decision_study(0.50), "real gap") + + +def test_assert_more_powerful_fails_on_a_gap_that_is_noise(): + """Two points of rejection rate at 400 replicates is a coin flip. + + This is the assertion the gate replaces: ``strong > weak`` is satisfied by + whichever method the seed favoured, and reports it as a finding. + """ + strong, weak = _decision_study(0.52), _decision_study(0.50) + assert strong.rejection_rate > weak.rejection_rate + with pytest.raises(AssertionError, match="not measurably above"): + assert_more_powerful(strong, weak, "two points of rate") + + +def test_assert_more_powerful_fails_when_the_arguments_are_the_wrong_way_round(): + """A comparison wired up backwards must not pass.""" + with pytest.raises(AssertionError, match="not measurably above"): + assert_more_powerful(_decision_study(0.30), _decision_study(0.80), "backwards") + + +def test_assert_more_powerful_needs_decisions_on_both_studies(): + """One side without decisions is not a comparison.""" + with pytest.raises(ValueError, match="`less` study recorded no reject"): + assert_more_powerful(_decision_study(0.8), _interval_study(), "one side") + + +# -------------------------------------------------------------------------- +# Interval endpoints on the result object. +# -------------------------------------------------------------------------- + + +def test_one_endpoint_without_the_other_is_rejected(): + """Half an interval is not an interval.""" + with pytest.raises(ValueError, match="needs both endpoints"): + MonteCarloResult( + estimates=np.zeros(10), + standard_errors=np.ones(10), + covered=None, + rejected=None, + truth=0.0, + lowers=np.full(10, -1.0), + ) + + +def test_a_backwards_interval_is_rejected(): + """A negative width is not a width.""" + with pytest.raises(ValueError, match="upper endpoint below the lower"): + MonteCarloResult( + estimates=np.zeros(10), + standard_errors=np.ones(10), + covered=None, + rejected=None, + truth=0.0, + lowers=np.full(10, 1.0), + uppers=np.full(10, -1.0), + ) + + +def test_covered_contradicting_the_endpoints_is_rejected(): + """Two measurements of the same thing that disagree are one measurement too many.""" + with pytest.raises(ValueError, match="disagrees with the endpoints"): + MonteCarloResult( + estimates=np.zeros(10), + standard_errors=np.ones(10), + covered=np.zeros(10, dtype=bool), + rejected=None, + truth=0.0, + lowers=np.full(10, -1.0), + uppers=np.full(10, 1.0), + ) + + +def test_covered_is_filled_in_from_the_endpoints(): + """Endpoints determine coverage exactly, so a study with them has a rate.""" + study = _interval_study(reps=100, half_width=Z_95 * 0.1) + assert study.covered is not None + assert 0.85 <= study.coverage <= 1.0 + + +def test_a_study_without_endpoints_has_no_widths(): + """Zero would read as an infinitely precise estimator, not an unmeasured one.""" + with pytest.raises(ValueError, match="recorded no interval endpoints"): + _ = _study().mean_width + + # -------------------------------------------------------------------------- # The band itself. # -------------------------------------------------------------------------- @@ -278,7 +726,9 @@ def test_the_gates_still_fire_under_python_o(): import sys from simcheck import ( MonteCarloResult, assert_count_rate, assert_coverage, - assert_proportion, assert_se_calibrated, assert_unbiased, + assert_intervals_informative, assert_more_powerful, assert_narrower, + assert_power, assert_proportion, assert_se_calibrated, + assert_unbiased, ) import numpy as np @@ -307,6 +757,28 @@ def check(name, fn): check("coverage", lambda: assert_coverage(biased, 0.95)) check("se_calibrated", lambda: assert_se_calibrated(biased)) + estimates = np.linspace(0.9, 1.1, 400) + def interval(half, rejection): + flags = np.zeros(400, dtype=bool) + flags[: round(rejection * 400)] = True + return MonteCarloResult( + estimates=estimates, + standard_errors=np.full(400, 0.058), + covered=None, + rejected=flags, + truth=1.0, + lowers=estimates - half, + uppers=estimates + half, + ) + + vacuous = interval(1.0, 0.5) + check("intervals_informative", + lambda: assert_intervals_informative(vacuous, 0.95)) + check("narrower", lambda: assert_narrower(vacuous, interval(0.2, 0.5))) + check("power", lambda: assert_power(interval(1.0, 0.30), 0.80)) + check("more_powerful", + lambda: assert_more_powerful(interval(1.0, 0.50), interval(1.0, 0.49))) + print(",".join(sorted(fired))) """) result = subprocess.run( # noqa: S603 @@ -319,6 +791,10 @@ def check(name, fn): assert fired == { "count_rate", "coverage", + "intervals_informative", + "more_powerful", + "narrower", + "power", "proportion", "se_calibrated", "unbiased", diff --git a/tests/test_runner.py b/tests/test_runner.py index 77fb708..d783d90 100644 --- a/tests/test_runner.py +++ b/tests/test_runner.py @@ -15,6 +15,7 @@ from simcheck import ( Estimate, assert_coverage, + assert_intervals_informative, assert_unbiased, monte_carlo, ) @@ -94,6 +95,38 @@ def test_a_study_without_intervals_has_no_coverage_rather_than_zero_coverage(): assert_unbiased(result, "sample mean without intervals") +def test_the_runner_keeps_the_endpoints_and_not_only_the_hit(): + """Reducing each interval to a boolean throws away the width. + + Coverage cannot distinguish a calibrated interval from one so wide it could + not have failed, and the endpoints are the only evidence that would. They + used to be computed, used once and dropped on the floor. + """ + result = monte_carlo(_mean_with_interval, TRUTH, 200, seed=5) + + assert result.lowers is not None + assert result.uppers is not None + np.testing.assert_array_equal( + result.covered, (result.lowers <= TRUTH) & (result.uppers >= TRUTH) + ) + # Every interval is 1.96 standard errors either side, so the mean width is + # 2 * 1.96 * the mean reported standard error, exactly. + assert result.mean_width == pytest.approx( + 2 * 1.96 * float(np.mean(result.standard_errors)), rel=1e-12 + ) + assert_intervals_informative(result, 0.95, "1.96 interval at n=40") + + +def test_a_study_without_intervals_has_no_widths(): + """Zero width would read as an infinitely precise estimator.""" + result = monte_carlo(_mean_without_interval, TRUTH, 50, seed=0) + assert result.lowers is None + with pytest.raises(ValueError, match="recorded no interval endpoints"): + _ = result.mean_width + with pytest.raises(ValueError, match="recorded no interval endpoints"): + assert_intervals_informative(result, 0.95, "no intervals") + + def test_a_study_without_decisions_has_no_rejection_rate(): """Same argument for size and power.""" result = monte_carlo(_mean_without_interval, TRUTH, 50, seed=0) From 7743eb46826dab8d104dd29cf6e79d83d8582f22 Mon Sep 17 00:00:00 2001 From: ***** <721466+soodoku@users.noreply.github.com> Date: Sat, 8 Aug 2026 23:52:49 -0700 Subject: [PATCH 2/4] Fix three defects Codex found in the new gates All three were reproduced before being fixed and watched failing after the fix was reverted. **A width comparison could be won on one replicate.** `_mean_gap_se` returned zero when either study had fewer than two replicates, and a zero standard error reads as certainty, so `assert_narrower` certified whichever method the single draw happened to favour. It now refuses the comparison. **The plug-in Wald standard error is exactly zero at a rejection rate of 0 or 1.** One replicate rejecting against one not rejecting was a gap of 1.0 with no uncertainty, so `assert_more_powerful` certified a three-sigma difference from two observations. Agresti-Caffo gives that case 2.6 sigma and fails it, while total separation over four hundred replicates is still hundreds of sigma. **The derived SE tolerance assumed normality.** The relative noise in a sample standard deviation is `sqrt((kappa-1)/(4n))`, not `1/sqrt(2n)`, so a calibrated estimator with Student t(5) sampling error was flagged in 19 of 200 studies against the 0.3% three sigma advertises -- the "gate that fires on correct code" failure this package exists to prevent. The kurtosis is now estimated from the study and floored at 3, since the sample kurtosis is downward-biased at Monte Carlo replicate counts and must not be allowed to narrow the band. 2 of 200 now. geoinference, the only consumer relying on the default tolerance, still passes: 43 passed, 4 subtests passed. Co-Authored-By: Claude Opus 5 (1M context) --- CHANGELOG.md | 10 ++- README.md | 6 +- src/simcheck/gates.py | 119 ++++++++++++++++++++++++++++++++--- tests/test_negative.py | 138 +++++++++++++++++++++++++++++++++++++++++ 4 files changed, 260 insertions(+), 13 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 5617d46..8d7aff5 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -86,8 +86,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 was `0.15`, the one number in the package chosen by hand rather than derived, and it was wrong in both directions: `se_ratio` divides a mean of `reps` reported standard errors by a sample standard deviation of `reps` estimates, so - its Monte Carlo spread is `sqrt(cv^2/reps + 1/(2(reps-1)))`, and three of those - is 0.21 at 100 replicates and 0.05 at 2000. The fixed value was therefore tight + its Monte Carlo spread is `sqrt(cv^2/reps + (kappa-1)/(4*reps))`, and three of + those is 0.21 at 100 replicates and 0.05 at 2000 for a normal estimator. The fixed value was therefore tight enough to fail correct estimators in a fast tier and loose enough to certify a 12% error in a deep one. Passing `tolerance=` explicitly still overrides it, and `se_ratio_tolerance(result)` returns the derived band. @@ -97,6 +97,12 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 against this branch: 43 passed, 4 subtests passed. Every other consumer passes `tolerance=` explicitly or does not call the gate. + The band uses the estimator's own fourth moment, not a normal assumption: + `Var(s)/sigma^2 = (kappa-1)/(4*reps)`, with `kappa` estimated from the study and + floored at 3 so a downward-biased sample kurtosis cannot narrow it. Assuming + normality flagged a *calibrated* estimator with Student t(5) sampling error in + 19 of 200 studies; with the fourth-moment term it is 2 of 200. + ### Fixed - **`assert_se_calibrated` diagnosed a missing standard error as a constant diff --git a/README.md b/README.md index a6d30e8..af33bab 100644 --- a/README.md +++ b/README.md @@ -189,8 +189,10 @@ the package chosen by hand rather than derived — and it was wrong in both directions at once. `se_ratio` is `mean(reported se) / sd(estimates)`, and both halves are estimated from the same replicates, so it is noisy even when the estimator is perfect: the numerator's relative standard error is `cv/sqrt(reps)` -and the denominator's is `1/sqrt(2(reps-1))`. Added in quadrature and taken at -three sigma, that is 0.21 at 100 replicates and 0.05 at 2000. A fixed 0.15 was +and the denominator's is `sqrt((κ-1)/(4·reps))`, where κ is the estimator's +kurtosis — measured, not assumed, because a normal assumption flags a correct +heavy-tailed estimator. Added in quadrature and taken at three sigma, that is +0.21 at 100 replicates and 0.05 at 2000 for a normal estimator. A fixed 0.15 was therefore tight enough to fail correct estimators in a fast tier and loose enough to certify a 12% error in a deep one. diff --git a/src/simcheck/gates.py b/src/simcheck/gates.py index 007a29b..c5c6654 100644 --- a/src/simcheck/gates.py +++ b/src/simcheck/gates.py @@ -273,21 +273,94 @@ def _mean_gap_se(first: np.ndarray, second: np.ndarray) -> float: this is the conservative choice: a paired comparison passing this gate would also pass a paired one. + A study of one replicate has no estimable variance, and reporting zero for it + would say the difference is known exactly -- so one replicate would certify + whichever method the single draw happened to favour. Callers check the + replicate count before getting here. + Args: - first: Per-replicate values from one study. - second: Per-replicate values from the other. + first: Per-replicate values from one study, at least two of them. + second: Per-replicate values from the other, at least two of them. Returns: float: The standard error of ``mean(first) - mean(second)``. + + Raises: + ValueError: If either study has fewer than two replicates. """ if len(first) < 2 or len(second) < 2: - return 0.0 + raise ValueError( + f"a difference of means needs at least two replicates in each " + f"study, got {len(first)} and {len(second)}" + ) return math.sqrt( float(np.var(first, ddof=1)) / len(first) + float(np.var(second, ddof=1)) / len(second) ) +def _relative_sd_error(values: np.ndarray) -> float: + """Relative Monte Carlo standard error of a sample standard deviation. + + The delta method gives ``Var(s) / sigma^2 = (kappa - 1) / (4 * n)``, where + ``kappa = mu_4 / sigma^4`` is the kurtosis of whatever is being averaged. For + a normal sample ``kappa = 3`` and this is the familiar ``1 / (2n)``; for a + heavy-tailed one it is larger, and using the normal figure would understate + the noise in ``sampling_sd`` and flag correct estimators. A calibrated + estimator with scaled Student t(5) sampling error was flagged in about 10% of + studies by the normal-only version, against the 0.3% that three sigma + advertises. + + ``kappa`` is estimated from the sample but never allowed below 3. The sample + kurtosis is an eighth-moment quantity and is heavily downward-biased at the + replicate counts a Monte Carlo study runs at, so letting a low estimate + *narrow* the band would reintroduce the same failure through the back door, + while a high one widening it is exactly the correction that is wanted. + + Args: + values: The sample whose standard deviation's noise is wanted. + + Returns: + float: The standard error of ``s``, as a fraction of ``s``. + + Raises: + ValueError: If there are fewer than two values, which leaves no spread. + """ + count = len(values) + if count < 2: + raise ValueError(f"a standard deviation needs at least two values, got {count}") + sample = np.asarray(values, dtype=float) + spread = float(np.std(sample, ddof=1)) + kurtosis = 3.0 + if spread: + standardised = (sample - float(np.mean(sample))) / spread + kurtosis = max(float(np.mean(standardised**4)), 3.0) + return math.sqrt((kurtosis - 1.0) / (4.0 * count)) + + +def _agresti_caffo_variance(rate: float, reps: int) -> float: + """Variance of an observed rate, with one success and one failure added. + + The plug-in Wald variance ``p(1-p)/n`` is exactly zero at ``p`` of 0 or 1, + which says a rate observed to be 1.0 is known with certainty -- so a single + rejecting replicate against a single non-rejecting one would be a difference + of 1.0 with no uncertainty at all. Adding a success and a failure to each arm + (Agresti and Caffo, 2000) keeps the variance positive where it matters and is + negligible where it does not: at 400 replicates it changes the standard error + in the fourth decimal. + + Args: + rate: The observed rate. + reps: Replicates it was observed over. + + Returns: + float: The variance of the adjusted rate. + """ + adjusted_reps = reps + 2 + adjusted_rate = (rate * reps + 1.0) / adjusted_reps + return adjusted_rate * (1.0 - adjusted_rate) / adjusted_reps + + def _gap_is_unresolved(gap: float, standard_error: float, sigmas: float) -> bool: """Whether a difference is too small for the studies to resolve. @@ -430,8 +503,12 @@ def se_ratio_tolerance(result: MonteCarloResult, sigmas: float = GATE_SIGMAS) -> coefficient of variation across replicates. An estimator that reports the same standard error every time contributes nothing here. * The denominator is a sample standard deviation of ``reps`` draws, whose - relative standard error is ``1 / sqrt(2 * (reps - 1))`` -- exactly for - normal estimates and closely for anything with a finite fourth moment. + relative standard error is ``sqrt((kappa - 1) / (4 * reps))`` for an + estimator with kurtosis ``kappa``. That is ``1 / sqrt(2 * reps)`` for a + normal estimator and larger for a heavy-tailed one, and the kurtosis is + estimated from the study rather than assumed: see ``_relative_sd_error``. + Assuming normality here flagged a correct estimator with Student t(5) + sampling error in about 10% of studies. Adding them in quadrature and multiplying by ``sigmas`` gives the band. The two are in fact positively correlated for most estimators -- a replicate that @@ -467,7 +544,9 @@ def se_ratio_tolerance(result: MonteCarloResult, sigmas: float = GATE_SIGMAS) -> variation = ( float(np.std(reported, ddof=1)) / mean_reported if mean_reported else 0.0 ) - relative = math.sqrt(variation**2 / result.reps + 1.0 / (2.0 * (result.reps - 1))) + relative = math.sqrt( + variation**2 / result.reps + _relative_sd_error(result.estimates) ** 2 + ) return sigmas * relative @@ -609,9 +688,12 @@ def assert_intervals_informative( ) widths = result.widths mean_width = result.mean_width + # The ratio's numerator is a mean width and its denominator a sampling + # standard deviation, so its Monte Carlo noise is the two in quadrature -- + # the same arithmetic as `se_ratio_tolerance`, on the same reasoning. relative = math.sqrt( (float(np.var(widths, ddof=1)) / result.reps) / mean_width**2 - + 1.0 / (2.0 * (result.reps - 1)) + + _relative_sd_error(result.estimates) ** 2 if mean_width else 0.0 ) @@ -659,7 +741,10 @@ def assert_narrower( sigmas: How many Monte Carlo standard errors the gap must exceed. Raises: - ValueError: If either study recorded no interval endpoints. + ValueError: If either study recorded no interval endpoints, or has fewer + than two replicates -- one replicate has no estimable spread, so the + gate would certify whichever method the single draw happened to + favour. AssertionError: If the narrower study's intervals are not measurably narrower. """ @@ -669,6 +754,12 @@ def assert_narrower( f"{label or 'this comparison'}: the `{name}` study recorded no " "interval endpoints, so the two widths cannot be compared" ) + if study.reps < 2: + raise ValueError( + f"{label or 'this comparison'}: the `{name}` study has " + f"{study.reps} replicate, which has no spread. A width gap " + "measured on one draw is not a gap." + ) gap = wide.mean_width - narrow.mean_width standard_error = _mean_gap_se(wide.widths, narrow.widths) if not _gap_is_unresolved(gap, standard_error, sigmas): @@ -752,6 +843,15 @@ def assert_more_powerful( Comparing rejection rates under *different* alternatives compares the alternatives, not the tests. + The standard error of the difference is the Agresti-Caffo one -- one success + and one failure added to each arm before the Wald formula -- rather than the + plug-in Wald standard error, which **collapses to exactly zero at a rejection + rate of 0 or 1**. One replicate rejecting against one replicate not rejecting + gives a gap of 1.0 with a plug-in standard error of 0, so the plain Wald + version certifies a three-sigma difference from two observations. Agresti- + Caffo gives that case 2.6 sigma, so it fails, while total separation over + four hundred replicates is still hundreds of sigma and still passes. + Args: more: The study claimed to be more powerful. less: The study it is claimed to beat. @@ -771,7 +871,8 @@ def assert_more_powerful( strong, weak = more.rejection_rate, less.rejection_rate gap = strong - weak standard_error = math.sqrt( - strong * (1.0 - strong) / more.reps + weak * (1.0 - weak) / less.reps + _agresti_caffo_variance(strong, more.reps) + + _agresti_caffo_variance(weak, less.reps) ) if not _gap_is_unresolved(gap, standard_error, sigmas): return diff --git a/tests/test_negative.py b/tests/test_negative.py index 0cedb4e..fb23646 100644 --- a/tests/test_negative.py +++ b/tests/test_negative.py @@ -317,6 +317,103 @@ def test_the_same_se_error_passes_a_shallow_study_and_fails_a_deep_one(): assert_se_calibrated(_study(reps=2000, se_scale=1.12), "0.15", tolerance=0.15) +def _heavy_tailed_study(reps: int, seed: int) -> MonteCarloResult: + """A perfectly calibrated estimator whose sampling error is Student t(5). + + Scaled to unit variance, so a reported standard error of 1.0 is exactly + right and every gate should stay silent. + + Parameters + ---------- + reps + Replicates. + seed + Generator seed. + + Returns + ------- + MonteCarloResult + The constructed study. + """ + rng = np.random.default_rng(seed) + estimates = rng.standard_t(5, reps) / np.sqrt(5 / 3) + return MonteCarloResult( + estimates=estimates, + standard_errors=np.ones(reps), + covered=None, + rejected=None, + truth=0.0, + ) + + +def test_the_se_tolerance_follows_the_fourth_moment(): + """The noise in a sample standard deviation depends on the fourth moment. + + ``Var(s)/sigma^2 = (kappa - 1) / (4n)``, which collapses to the familiar + ``1/(2n)`` only when ``kappa`` is 3. Assuming normality makes the band too + narrow for a heavy-tailed estimator and flags it for being correct. + """ + reps = 400 + heavy = _heavy_tailed_study(reps, 3) + estimates = heavy.estimates + kurtosis = float( + np.mean(((estimates - estimates.mean()) / estimates.std(ddof=1)) ** 4) + ) + assert kurtosis > 3.0 + + # The reported standard errors are constant here, so the numerator of the + # ratio contributes no noise and the whole band is the spread of the + # sampling standard deviation. + assert se_ratio_tolerance(heavy) == pytest.approx( + 3.0 * math.sqrt((kurtosis - 1.0) / (4 * reps)), rel=1e-12 + ) + assert se_ratio_tolerance(heavy) > 3.0 / math.sqrt(2 * reps) + + +def test_the_se_tolerance_never_narrows_below_the_normal_case(): + """A low sample kurtosis must not tighten the band. + + The sample kurtosis is an eighth-moment quantity and is badly + downward-biased at Monte Carlo replicate counts, so letting a low estimate + narrow the band would reintroduce the very failure the fourth-moment term + exists to prevent. + """ + reps = 200 + # A two-point distribution has kurtosis 1, the smallest there is. + two_point = np.tile([-1.0, 1.0], reps // 2) + study = MonteCarloResult( + estimates=two_point, + standard_errors=np.ones(reps), + covered=None, + rejected=None, + truth=0.0, + ) + assert se_ratio_tolerance(study) == pytest.approx( + 3.0 / math.sqrt(2 * reps), rel=1e-12 + ) + + +def test_a_calibrated_heavy_tailed_estimator_is_not_flagged(): + """A gate that fires on 5% of correct code gets disabled within a week. + + Student t(5) sampling error is calibrated here by construction: the reported + standard error is exactly the true one. With the band computed as though the + estimates were normal, 19 of these 200 studies were flagged. The bound below + is the package's own stated standard for a false-positive rate, not the + number that happened to come out. + """ + flagged = 0 + for seed in range(200): + try: + assert_se_calibrated(_heavy_tailed_study(400, seed), f"seed {seed}") + except AssertionError: + flagged += 1 + assert flagged <= 10, f"{flagged} of 200 calibrated t(5) studies were flagged" + + # And the specific study that exposed this must pass. + assert_se_calibrated(_heavy_tailed_study(400, 0), "the reported counterexample") + + def test_assert_se_calibrated_says_so_when_no_standard_error_was_reported(): """A NaN standard error is an absent measurement, not a zero spread. @@ -502,6 +599,21 @@ def test_assert_narrower_fails_when_the_arguments_are_the_wrong_way_round(): assert_narrower(loose, tight, "backwards") +def test_assert_narrower_refuses_a_one_replicate_comparison(): + """One draw has no spread, so a gap measured on it is not a gap. + + The Monte Carlo standard error of the difference is zero when either study + has a single replicate, and a zero standard error reads as certainty -- so + the gate would certify whichever method the one draw happened to favour. + """ + single = _interval_study(reps=1, half_width=0.1) + pair = _interval_study(reps=2, half_width=0.2) + with pytest.raises(ValueError, match="which has no spread"): + assert_narrower(single, pair, "one replicate") + with pytest.raises(ValueError, match="which has no spread"): + assert_narrower(pair, single, "one replicate on the other side") + + def test_assert_narrower_needs_endpoints_on_both_studies(): """Comparing a width against an unmeasured one is not a comparison.""" with pytest.raises(ValueError, match="`wide` study recorded no interval"): @@ -578,6 +690,32 @@ def test_assert_more_powerful_fails_when_the_arguments_are_the_wrong_way_round() assert_more_powerful(_decision_study(0.30), _decision_study(0.80), "backwards") +def test_a_rejection_rate_of_one_is_not_known_with_certainty(): + """The plug-in Wald standard error is exactly zero at a rate of 0 or 1. + + One replicate rejecting against one not rejecting is a gap of 1.0 with a + plug-in standard error of 0, which passes any multiple of it -- a three-sigma + claim from two observations. The Agresti-Caffo standard error gives that case + 2.6 sigma, so it fails. + """ + with pytest.raises(AssertionError, match="not measurably above"): + assert_more_powerful( + _decision_study(1.0, reps=1), _decision_study(0.0, reps=1), "one each" + ) + + +def test_total_separation_over_a_real_study_still_passes(): + """The guard on the guard: the boundary fix must not blind the gate. + + All four hundred replicates rejecting against none of them is the largest + difference there is, and a correction that failed here would have traded one + defect for another. + """ + assert_more_powerful( + _decision_study(1.0, reps=400), _decision_study(0.0, reps=400), "separated" + ) + + def test_assert_more_powerful_needs_decisions_on_both_studies(): """One side without decisions is not a comparison.""" with pytest.raises(ValueError, match="`less` study recorded no reject"): From 37e350d56b146b3e333f67b084d60ac302573b76 Mon Sep 17 00:00:00 2001 From: ***** <721466+soodoku@users.noreply.github.com> Date: Sun, 9 Aug 2026 00:13:40 -0700 Subject: [PATCH 3/4] Fix two more defects the hosted Codex review found Both reproduced before the fix and watched failing after it was reverted. **Agresti-Caffo was applied to the variance but not to the gap**, which is not Agresti-Caffo and leaves the hole open whenever the two studies are different sizes: one replicate rejecting against a hundred not rejecting is a raw gap of 1.0, and it cleared three sigma against an adjusted standard error. Both halves now come from the adjusted rates, which puts that case at 2.4 sigma and the one-against-one case at 0.9, while total separation over four hundred replicates is unchanged at hundreds of sigma. **`se_ratio_tolerance` returned NaN for a study that reported no standard error.** The helper is public and the README points callers at it, so the obvious hand-written check -- `if abs(ratio - 1) > tolerance: raise` -- compared against NaN, got False, and passed silently on a study that measured nothing. It raises now. The gate itself was never affected: it rejects a non-finite ratio first. Co-Authored-By: Claude Opus 5 (1M context) --- CHANGELOG.md | 9 +++++-- src/simcheck/gates.py | 58 ++++++++++++++++++++++++++++++++++-------- tests/test_negative.py | 32 +++++++++++++++++++++++ 3 files changed, 86 insertions(+), 13 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 8d7aff5..e2d0ae8 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -61,14 +61,19 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 distribution separates conservatism from vacuity on width alone. - `assert_narrower`, for the efficiency half of an interval comparison, banded by - the Monte Carlo standard error of the difference in mean width. + the Monte Carlo standard error of the difference in mean width. A study of one + replicate is refused rather than treated as certain. - `assert_power` and `assert_more_powerful`. The package documented power as one of the four questions it answers and had no gate for it; consumers were reaching for `assert_proportion`, which is two-sided and needs a nominal you already know analytically, or hand-rolling a two-sample standard error. `assert_power` is one-sided, because power is a floor and a two-sided band - fails a test for being better than claimed. + fails a test for being better than claimed. `assert_more_powerful` compares + Agresti-Caffo adjusted rates rather than raw ones: the plug-in Wald standard + error is exactly zero at a rejection rate of 0 or 1, so one replicate rejecting + against one not rejecting would otherwise be a three-sigma finding from two + observations. - **Negative tests for every gate.** Each is exercised on input that satisfies its property, where it must stay silent, and on input that violates it, where diff --git a/src/simcheck/gates.py b/src/simcheck/gates.py index c5c6654..23a35ba 100644 --- a/src/simcheck/gates.py +++ b/src/simcheck/gates.py @@ -338,6 +338,19 @@ def _relative_sd_error(values: np.ndarray) -> float: return math.sqrt((kurtosis - 1.0) / (4.0 * count)) +def _agresti_caffo_rate(rate: float, reps: int) -> float: + """An observed rate with one success and one failure added to it. + + Args: + rate: The observed rate. + reps: Replicates it was observed over. + + Returns: + float: The adjusted rate. + """ + return (rate * reps + 1.0) / (reps + 2) + + def _agresti_caffo_variance(rate: float, reps: int) -> float: """Variance of an observed rate, with one success and one failure added. @@ -532,7 +545,12 @@ def se_ratio_tolerance(result: MonteCarloResult, sigmas: float = GATE_SIGMAS) -> Raises: ValueError: If the study has fewer than two replicates, which leaves the - spread -- and so the ratio -- undefined. + spread -- and so the ratio -- undefined; or if it reported no usable + standard error, in which case there is no ratio for the band to be a + band around. Returning NaN in that second case would be worse than + raising: a caller writing the obvious check by hand, + ``if abs(ratio - 1) > tolerance: raise``, gets ``False`` from every + comparison with NaN and so passes silently. """ if result.reps < 2: raise ValueError( @@ -541,6 +559,12 @@ def se_ratio_tolerance(result: MonteCarloResult, sigmas: float = GATE_SIGMAS) -> ) reported = np.asarray(result.standard_errors, dtype=float) mean_reported = float(np.mean(reported)) + if not np.isfinite(mean_reported): + raise ValueError( + f"this study's reported standard errors average to {mean_reported}, " + "so there is no ratio to put a band around. Have the estimator " + "report a standard error on every replicate." + ) variation = ( float(np.std(reported, ddof=1)) / mean_reported if mean_reported else 0.0 ) @@ -843,14 +867,20 @@ def assert_more_powerful( Comparing rejection rates under *different* alternatives compares the alternatives, not the tests. - The standard error of the difference is the Agresti-Caffo one -- one success - and one failure added to each arm before the Wald formula -- rather than the - plug-in Wald standard error, which **collapses to exactly zero at a rejection - rate of 0 or 1**. One replicate rejecting against one replicate not rejecting - gives a gap of 1.0 with a plug-in standard error of 0, so the plain Wald - version certifies a three-sigma difference from two observations. Agresti- - Caffo gives that case 2.6 sigma, so it fails, while total separation over - four hundred replicates is still hundreds of sigma and still passes. + The comparison is the Agresti-Caffo one: a success and a failure are added to + each arm, and **both** the gap and its standard error are computed from the + adjusted rates. The plug-in Wald standard error, ``p(1-p)/n``, is exactly zero + at a rejection rate of 0 or 1, so one replicate rejecting against one not + rejecting would be a difference of 1.0 with no uncertainty at all -- a + three-sigma claim from two observations. Adjusting only the variance and not + the gap leaves the same hole open when the two studies are different sizes: + 1 of 1 against 0 of 100 is a raw gap of 1.0, which clears three sigma against + an adjusted standard error. + + Adjusted, those two cases are 0.87 and 2.41 sigma and both fail, while total + separation over four hundred replicates is still hundreds of sigma and + passes. The adjustment is negligible wherever the answer is not in doubt: at + 400 replicates it moves a rate by a quarter of a percentage point. Args: more: The study claimed to be more powerful. @@ -869,7 +899,12 @@ def assert_more_powerful( "reject/accept decisions, so the two cannot be compared" ) strong, weak = more.rejection_rate, less.rejection_rate - gap = strong - weak + # The gap comes from the adjusted rates, not the raw ones. Pairing a raw gap + # with an adjusted variance is not the Agresti-Caffo interval and inherits + # the defect it is there to fix: one replicate rejecting against a hundred + # not rejecting is a raw gap of 1.0 with an adjusted standard error, which + # clears three sigma on the strength of a single draw. + gap = _agresti_caffo_rate(strong, more.reps) - _agresti_caffo_rate(weak, less.reps) standard_error = math.sqrt( _agresti_caffo_variance(strong, more.reps) + _agresti_caffo_variance(weak, less.reps) @@ -884,5 +919,6 @@ def assert_more_powerful( raise AssertionError( f"{label}: rejection rate {strong:.4f} over {more.reps} replicates is " f"not measurably above {weak:.4f} over {less.reps}: the gap of " - f"{gap:+.4f} is {resolved}, against a gate of {sigmas:g}" + f"{gap:+.4f} between the adjusted rates is {resolved}, against a gate of " + f"{sigmas:g}" ) diff --git a/tests/test_negative.py b/tests/test_negative.py index fb23646..96d3194 100644 --- a/tests/test_negative.py +++ b/tests/test_negative.py @@ -414,6 +414,24 @@ def test_a_calibrated_heavy_tailed_estimator_is_not_flagged(): assert_se_calibrated(_heavy_tailed_study(400, 0), "the reported counterexample") +def test_the_se_tolerance_refuses_a_study_with_no_standard_errors(): + """A NaN tolerance is worse than no tolerance at all. + + The helper is public, so a caller may write the obvious check by hand -- + ``if abs(ratio - 1) > tolerance: raise`` -- and every comparison against NaN + is False, so that check passes silently on a study that measured nothing. + """ + study = MonteCarloResult( + estimates=np.linspace(0.9, 1.1, 100), + standard_errors=np.full(100, np.nan), + covered=None, + rejected=None, + truth=1.0, + ) + with pytest.raises(ValueError, match="no ratio to put a band around"): + se_ratio_tolerance(study) + + def test_assert_se_calibrated_says_so_when_no_standard_error_was_reported(): """A NaN standard error is an absent measurement, not a zero spread. @@ -704,6 +722,20 @@ def test_a_rejection_rate_of_one_is_not_known_with_certainty(): ) +def test_one_replicate_cannot_beat_a_hundred(): + """Adjusting the variance but not the gap leaves the same hole open. + + One replicate rejecting against a hundred not rejecting is a *raw* gap of + 1.0, which clears three sigma against an Agresti-Caffo standard error. The + adjustment has to be applied to the gap as well, which puts this case at 2.4 + sigma. + """ + with pytest.raises(AssertionError, match="not measurably above"): + assert_more_powerful( + _decision_study(1.0, reps=1), _decision_study(0.0, reps=100), "1 v 100" + ) + + def test_total_separation_over_a_real_study_still_passes(): """The guard on the guard: the boundary fix must not blind the gate. From 6e8ae3b59c0c5b2871fd8d34c5d355269810b36a Mon Sep 17 00:00:00 2001 From: ***** <721466+soodoku@users.noreply.github.com> Date: Sun, 9 Aug 2026 00:20:15 -0700 Subject: [PATCH 4/4] Date the release for the day it is actually tagged The changelog entry was written before midnight rolled over. Co-Authored-By: Claude Opus 5 (1M context) --- CHANGELOG.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index e2d0ae8..77eb29a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,7 +7,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] -## [0.1.0] - 2026-08-08 +## [0.1.0] - 2026-08-09 ### Added