Designing an experiment before there is any data · step 5 of 6
When to stop
The witness runs for 110 hours. Do you have to wait?
A hundred and ten hours at one angle is a long time to learn nothing you could not have learned in thirty. The same sequential machinery that stops a clinical trial early applies here: look at the witness on a schedule, against a boundary that spends the error rate across the looks.
The boundary is what makes looking legitimate. Peeking at a running count and stopping when it looks good is how you manufacture a discovery; an O'Brien-Fleming boundary is the price of being allowed to look at all.
import math
from axiom.design import (
LookSchedule, StoppingRule, information_fractions,
obrien_fleming, operating_characteristics,
)
look_hours = (5.0, 15.0, 30.0, 60.0, 110.0)
looks = LookSchedule(
labels=tuple(f"{h:g} h" for h in look_hours),
information=information_fractions(look_hours),
)
boundary = obrien_fleming(0.05, looks, side="upper", kind="efficacy")
rule = StoppingRule(name="witness_150", looks=looks, boundaries=(boundary,))
print(f"{'look':>8} {'information':>12} {'stop above z':>13}")
for label, frac, z in zip(looks.labels, looks.information, boundary.z):
print(f"{label:>8} {frac:12.2f} {z:13.3f}")
signal = float(
s.rate_per_steradian([150.0], s.HARD_CENTRE)[0] - s.BACKGROUND_DENSITY
) * s.OMEGA_MAX
background = s.BACKGROUND_DENSITY * s.OMEGA_MAX
full = 110 * 3600.0
print()
for factor in (1.0, 1e-3, 1e-4):
drift = signal * factor * full / math.sqrt((signal * factor + background) * full)
oc = operating_characteristics(rule, drift=drift)
print(f"a source {factor:7.0e} of the modelled one: "
f"P(stop) = {oc.crossings.stop_probability:.3f}")
look information stop above z
5 h 0.05 7.914
15 h 0.14 4.569
30 h 0.27 3.231
60 h 0.55 2.285
110 h 1.00 1.687
a source 1e+00 of the modelled one: P(stop) = 1.000
a source 1e-03 of the modelled one: P(stop) = 1.000
a source 1e-04 of the modelled one: P(stop) = 0.205
Every line above was captured by running this code at build time. The steps share one namespace, so the code on this page is the code you would type after the steps before it.