API

Monte Carlo tests for statistical estimators.

The tests most statistical code ships assert that it runs. The tests it needs assert the classical properties: if the assumptions hold, is the estimator unbiased, do its intervals cover at the nominal rate, is the test’s size right under the null, and does it have power under an alternative.

simcheck supplies the machinery for those four questions, with one rule running through it: the tolerance comes from the replicate count, never from a number chosen by hand. coverage > 0.85 for a nominal 95% interval is meaningless – far too loose at ten thousand replicates, tight enough to fail spuriously at fifty. Every gate here derives its band from the number of replicates, so the same assertion adapts to the tier it runs in and reports how far outside the band it fell.

Two failure modes this is built to prevent, both observed in the wild:

An aggregate threshold absorbing a systematic failure. A matrix test asserting success_rate >= 0.7 passed for releases while one input pattern in eight raised an exception for every configuration – 12.5% sits comfortably inside a 30% allowance. Assert the property, not a rate that has room to hide things.

An assertion helper that cannot fail. Every gate here has a negative test: an input that violates the property, and a check that the gate raises on it. A 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 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 assert_more_powerful() exists.

Examples: >>> import numpy as np >>> from simcheck import MonteCarloResult, assert_coverage, assert_unbiased >>> rng = np.random.default_rng(0) >>> truth, n, reps = 2.0, 40, 400 >>> draws = rng.normal(truth, 1.0, size=(reps, n)) >>> means = draws.mean(axis=1) >>> errors = draws.std(axis=1, ddof=1) / np.sqrt(n) >>> result = MonteCarloResult( … estimates=means, … standard_errors=errors, … covered=np.abs(means - truth) <= 1.96 * errors, … rejected=np.abs(means) > 1.96 * errors, … truth=truth, … ) >>> assert_unbiased(result, “sample mean”) >>> assert_coverage(result, 0.95, “sample mean”)

class simcheck.Estimate(value, standard_error=nan, lower=None, upper=None, rejected=None)

What one replicate of an estimator produced.

Parameters:

• value (float) – The point estimate.

• standard_error (float) – The standard error the estimator reported. Leave as NaN when it reports none; assert_se_calibrated() will then have nothing to check and will say so.

• lower (float | None) – Lower end of the interval, if the estimator produced one.

• upper (float | None) – Upper end of the interval.

• rejected (bool | None) – Whether the replicate rejected the null, if the estimator performs a test.

property has_interval: bool

Whether this replicate reported both interval endpoints.

class simcheck.MonteCarloResult(estimates, standard_errors, covered, rejected, truth, lowers=None, uppers=None)

Sampling behaviour of an estimator at one point.

Parameters:

• estimates (NDArray[float64]) – The estimate from each replicate.

• standard_errors (NDArray[float64]) – The standard error the estimator reported on each replicate.

• covered (NDArray[bool] | None) – Whether each replicate’s interval contained the truth, or None when the estimator reported no intervals. None rather than an array of False: the latter reads as coverage of zero, which is a broken estimator, not an absent measurement.

• rejected (NDArray[bool] | None) – Whether each replicate rejected the null, or None when the estimator performs no test.

• truth (float) – The true value of the quantity being estimated.

• lowers (NDArray[float64] | None) – 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 (NDArray[float64] | None) – 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, if one endpoint array is given without the other, if any interval runs backwards, or if covered contradicts the endpoints.

property reps: int

Number of replicates.

property bias: float

Mean deviation of the estimate from the truth.

property sampling_sd: float

The estimator’s actual spread across replicates.

property reported_se: float

The standard error the estimator claims, on average.

property se_ratio: float

Claimed standard error over actual spread. One means calibrated.

Below one, the estimator is overconfident and its intervals will under-cover; above one it is conservative.

property mc_se: float

Monte Carlo standard error of the mean estimate.

This is the precision of the study, not of the estimator: it is what makes bias interpretable, because a bias smaller than this is indistinguishable from simulation noise.

property bias_t: float

Bias in units of its own Monte Carlo standard error.

A t statistic for the null that the estimator is unbiased. Zero when the estimator is deterministic, in which case there is no sampling variation to test against.

property coverage: float

Fraction of replicates whose interval contained the truth.

Returns:

The coverage rate.

Return type:

float

Raises:

ValueError – If the study recorded no intervals. Reporting zero here would be indistinguishable from an estimator whose intervals never cover, which is the opposite conclusion.

property widths: NDArray[float64]

Width of each replicate’s interval.

Returns:

One width per replicate.

Return type:

numpy.ndarray

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.

property mean_width: float

Mean interval width, in the units of the estimand.

Reading this on a study that recorded no endpoints raises ValueError through widths, rather than reporting a width of zero.

Returns:

The mean width.

Return type:

float

property median_width: 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 widths.

Returns:

The median width.

Return type:

float

property rejection_rate: float

Fraction of replicates that rejected the null.

Under a true null this estimates the test’s size; under an alternative, its power.

Returns:

The rejection rate.

Return type:

float

Raises:

ValueError – If the study recorded no reject/accept decisions.

summary(label='')

A one-line report, for printing a table of results.

Parameters:

label (str) – Row label.

Returns:

One formatted line.

Return type:

str

simcheck.assert_count_rate(successes, reps, nominal, label='', sigmas=3.0)

Fail if a count of successes is inconsistent with the claimed rate.

Parameters:

• successes (int) – Number of successes, in [0, reps].

• reps (int) – Number of replicates.

• nominal (float) – The claimed rate.

• label (str) – Included in the failure message.

• sigmas (float) – Slack, in binomial standard errors.

Raises:

• ValueError – If successes is negative or exceeds reps.

• AssertionError – If the implied rate falls outside the band.

Return type:

None

simcheck.assert_coverage(result, nominal=0.95, label='', sigmas=3.0)

Fail if interval coverage is inconsistent with the nominal level.

Parameters:

• result (MonteCarloResult) – A completed Monte Carlo study.

• nominal (float) – The level the intervals claim.

• label (str) – Included in the failure message.

• sigmas (float) – Slack, in binomial standard errors.

Raises:

• AssertionError – If coverage falls outside the band.

• ValueError – If the study recorded no intervals, so there is no coverage rate to test.

Return type:

None

simcheck.assert_intervals_informative(result, nominal=0.95, label='', max_ratio=None, sigmas=3.0)

Fail if the intervals are so wide that their coverage means nothing.

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 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.

Parameters:

• result (MonteCarloResult) – A completed Monte Carlo study that recorded interval endpoints.

• nominal (float) – The level the intervals claim.

• label (str) – Included in the failure message.

• max_ratio (float | None) – Override for the derived width multiple.

• sigmas (float) – 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.

Return type:

None

simcheck.assert_more_powerful(more, less, label='', sigmas=3.0)

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.

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.

Parameters:

• more (MonteCarloResult) – The study claimed to be more powerful.

• less (MonteCarloResult) – The study it is claimed to beat.

• label (str) – Included in the failure message.

• sigmas (float) – 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.

Return type:

None

simcheck.assert_narrower(narrow, wide, label='', sigmas=3.0)

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 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.

Parameters:

• narrow (MonteCarloResult) – The study claimed to produce the narrower intervals.

• wide (MonteCarloResult) – The study it is claimed to beat.

• label (str) – Included in the failure message.

• sigmas (float) – How many Monte Carlo standard errors the gap must exceed.

Raises:

• 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.

Return type:

None

simcheck.assert_power(result, minimum, label='', sigmas=3.0)

Fail if a test rejects less often than claimed under an alternative.

One-sided, unlike 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 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.

Parameters:

• result (MonteCarloResult) – A completed Monte Carlo study run under the alternative.

• minimum (float) – The power being claimed, in [0, 1].

• label (str) – Included in the failure message.

• sigmas (float) – 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.

Return type:

None

simcheck.assert_proportion(observed, reps, nominal, label='', sigmas=3.0)

Fail if an observed rate is inconsistent with the claimed one.

Parameters:

• observed (float) – The observed rate, in [0, 1]. Pass a count to assert_count_rate() instead.

• reps (int) – Number of replicates the rate was computed over.

• nominal (float) – The claimed rate.

• label (str) – Included in the failure message.

• sigmas (float) – Slack, in binomial standard errors.

Raises:

• ValueError – If observed is not in [0, 1], which usually means a count was passed where a rate was expected.

• AssertionError – If the observed rate falls outside the band.

Return type:

None

simcheck.assert_se_calibrated(result, label='', tolerance=None, sigmas=3.0)

Fail if the reported standard error misstates the estimator’s spread.

Coverage can look correct while the reported standard error is wrong, if two 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 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.

Parameters:

• result (MonteCarloResult) – A completed Monte Carlo study.

• label (str) – Included in the failure message.

• tolerance (float | None) – Largest permitted relative deviation of se_ratio from one. Derived from reps when omitted.

• sigmas (float) – 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, if the estimator did not vary at all across replicates, or if it reported no standard error to check.

Return type:

None

simcheck.assert_unbiased(result, label='', sigmas=3.0)

Fail if the estimator’s mean is distinguishable from the truth.

The comparison is against the Monte Carlo standard error of the mean, so a bias too small for the study to resolve does not fail, and the study can be made more demanding simply by running more replicates.

Parameters:

• result (MonteCarloResult) – A completed Monte Carlo study.

• label (str) – Included in the failure message.

• sigmas (float) – How many Monte Carlo standard errors of slack to allow.

Raises:

AssertionError – If the bias t statistic exceeds the gate.

Return type:

None

simcheck.binomial_band(nominal, reps, sigmas=3.0)

The interval a well-calibrated rate should land in.

Parameters:

• nominal (float) – The rate the estimator claims, for instance 0.95 for coverage.

• reps (int) – Number of replicates.

• sigmas (float) – How many binomial standard errors of slack to allow.

Returns:

Lower and upper bounds, clipped to [0, 1].

Return type:

tuple of float

Raises:

ValueError – If nominal is outside [0, 1], reps is not positive, or sigmas is negative.

Examples: >>> low, high = binomial_band(0.95, 400) >>> round(low, 4), round(high, 4) (0.9173, 0.9827)

simcheck.deep_tier()

Whether the deep tier has been requested.

Returns:

True when SIMCHECK_DEEP is set to something truthy.

Return type:

bool

simcheck.monte_carlo(replicate, truth, reps, *, seed=0)

Run replicate many times and collect its sampling behaviour.

Parameters:

• replicate (Callable[[np.random.Generator], Estimate]) – Called once per replicate with its own generator. It should simulate a dataset, fit the estimator, and return what came out. It must not close over shared mutable state, or the replicates stop being independent.

• truth (float) – The true value of the quantity being estimated, which the caller knows because it generated the data.

• reps (int) – Number of replicates.

• seed (int) – Seed for the replicate stream. Replicate i is a function of (seed, i) alone, so it can be reproduced on its own and raising reps adds replicates rather than changing the existing ones.

Returns:

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.

Return type:

MonteCarloResult

Raises:

ValueError – If reps is not positive, or if the estimator reported an interval on some replicates but not others – which means it is not doing the same thing every time, and pooling the results would silently mix two estimators.

Examples

>>> import numpy as np
>>> from simcheck import Estimate, assert_unbiased, monte_carlo
>>> def draw(rng):
...     x = rng.normal(2.0, 1.0, 50)
...     se = x.std(ddof=1) / np.sqrt(50)
...     return Estimate(x.mean(), se, x.mean() - 2 * se, x.mean() + 2 * se)
>>> result = monte_carlo(draw, truth=2.0, reps=200, seed=0)
>>> result.reps
200
>>> assert_unbiased(result, "sample mean")

simcheck.reps_for(deep=None)

Replicate count for the current tier.

Parameters:

deep (bool | None) – Force a tier. When omitted, read the environment.

Returns:

Number of replicates to run.

Return type:

int

simcheck.se_ratio_tolerance(result, sigmas=3.0)

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 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 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.

Parameters:

• result (MonteCarloResult) – A completed Monte Carlo study.

• sigmas (float) – How many Monte Carlo standard errors of slack to allow.

Returns:

The largest deviation of se_ratio from one that this study cannot distinguish from noise.

Return type:

float

Raises:

ValueError – If the study has fewer than two replicates, which leaves the 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.

simcheck.vacuous_width_ratio(nominal, reps)

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. 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 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.

Parameters:

• nominal (float) – The level the intervals claim, for instance 0.95.

• reps (int) – Number of replicates in the study.

Returns:

The width multiple at which the study loses the ability to observe the interval failing.

Return type:

float

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

simcheck.width_ratio(result, nominal=0.95)

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.

Parameters:

• result (MonteCarloResult) – A completed Monte Carlo study that recorded interval endpoints.

• nominal (float) – The level the intervals claim.

Returns:

The ratio.

Return type:

float

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.