Statistics and Monte Carlo

Do the error bars mean what they say?

An inference code that reports 68% intervals should contain the true value in 68% of repeated experiments. coverage_check repeats the experiment: it draws a truth, simulates data, asks your code for an interval and counts how often the truth is inside. An exact binomial test decides; a correct code fails at most a fraction alpha (default 0.001) of the time.

"""Do the error bars mean what they say? 68% intervals should contain the truth 68% of the time."""
import math

from physverify.studies import coverage_check


def simulate(mu, rng):
    return rng.normal(mu, 2.0, size=25)  # 25 measurements with noise 2.0


def interval(sigma_assumed):
    def iv(data, level):
        half = 0.994457883 * sigma_assumed / math.sqrt(len(data))  # z for a 68% two-sided interval
        m = data.mean()
        return m - half, m + half
    return iv


def prior(rng):  # the truth is drawn from the prior the analysis assumes
    return rng.uniform(-10, 10)


print("right noise:    ", coverage_check(simulate, interval(2.0), prior, level=0.68).message)
print("noise too small:", coverage_check(simulate, interval(1.0), prior, level=0.68).message)

Output:

right noise:     268 of 400 intervals contain the truth (0.670, stated 0.68); binomial p = 6.7e-01, false-alarm rate 0.001
noise too small: 121 of 400 intervals contain the truth (0.302, stated 0.68); binomial p = 1.2e-53, false-alarm rate 0.001

The second interval assumes half the real noise, so it is far too narrow and covers the truth only 30% of the time.

Draw the truth from the analysis prior

For a Bayesian analysis, the truth must be drawn from the same prior the analysis assumes. Drawing it from a narrower range gave a spurious coverage of 0.865 for a nominal 0.68 interval in a real study. That is why the prior sampler is a required argument.

Monte Carlo samplers

physverify.stochastic checks samplers of the 2D Ising model against exact enumeration, detailed balance and a reference replay with the same random numbers. All 21 seeded bugs in its benchmark were caught with no false alarms. See Ising checks.