Bayesian Sequential Testing

Bayesian sequential testing evaluates evidence as it arrives, using a posterior probability or a Bayes factor rather than a p-value against a fixed \(\alpha\). At each look the choice is to stop for one hypothesis, stop for the other, or keep sampling.

It is often presented as immune to the problem that makes peeking costly in a frequentist test. That claim needs to be stated much more carefully than it usually is.

What optional stopping does and does not break

import numpy as np
from scipy.stats import norm

RUNS, MAX_N = 20_000, 2_000

def simulate(p_a, p_b, threshold, every, prior=(1, 1)):
    """Update two Beta posteriors after every user and stop the first time
       P(B > A) leaves the corridor. P(B > A) uses the normal approximation to
       the difference of two Beta posteriors, whose accuracy needs checking, especially at early looks with few successes."""
    rng = np.random.default_rng(0)
    a, b = prior
    ca = np.cumsum(rng.random((RUNS, MAX_N)) < p_a, axis=1)
    cb = np.cumsum(rng.random((RUNS, MAX_N)) < p_b, axis=1)
    looks = np.arange(every, MAX_N + 1, every)
    n = looks[None, :]
    ka, kb = ca[:, looks - 1], cb[:, looks - 1]
    var = lambda k: (a + k) * (b + n - k) / ((a + b + n) ** 2 * (a + b + n + 1))
    mean = lambda k: (a + k) / (a + b + n)
    prob = norm.cdf((mean(kb) - mean(ka)) / np.sqrt(var(ka) + var(kb)))
    crossed = (prob > threshold) | (prob < 1 - threshold)
    stopped = crossed.any(axis=1)
    first = np.where(stopped, crossed.argmax(axis=1), len(looks) - 1)
    called_b = prob[np.arange(RUNS), first] > threshold
    return float((called_b & stopped).mean()), looks[first], stopped

print(f"{'threshold':>11}{'looks taken':>18}{'calls B better':>17}{'median n':>11}")
for threshold, every in ((0.95, 2_000), (0.95, 500), (0.95, 100), (0.95, 20), (0.99, 20)):
    rate, stop, done = simulate(0.10, 0.10, threshold, every)
    label = "once, at the end" if every == MAX_N else f"every {every}"
    print(f"{threshold:>11}{label:>18}{rate:>17.4f}{np.median(stop):>11.0f}")

for prior in ((1, 1), (5, 45), (50, 450)):
    rate, stop, _ = simulate(0.10, 0.10, 0.95, 100, prior)
    print(f"    Beta{prior} on each arm ({prior[0] + prior[1]:>3} prior observations):"
          f" wrong-call rate {rate:.4f}, median n {np.median(stop):>5.0f}")

for threshold, every in ((0.95, 100), (0.99, 100)):
    rate, stop, done = simulate(0.10, 0.12, threshold, every)
    print(f"    threshold {threshold}: calls B better {rate:.4f}, median n"
          f" {np.median(stop):>5.0f}, reached a decision {done.mean():.4f}")
#   threshold       looks taken   calls B better   median n
#        0.95  once, at the end           0.0503       2000
#        0.95         every 500           0.1163       2000
#        0.95         every 100           0.2064       2000
#        0.95          every 20           0.2796       1140
#        0.99          every 20           0.0838       2000
#     Beta(1, 1) on each arm (  2 prior observations): wrong-call rate 0.2064, median n  2000
#     Beta(5, 45) on each arm ( 50 prior observations): wrong-call rate 0.1814, median n  2000
#     Beta(50, 450) on each arm (500 prior observations): wrong-call rate 0.0851, median n  2000
#     threshold 0.95: calls B better 0.7743, median n   800, reached a decision 0.8072
#     threshold 0.99: calls B better 0.5198, median n  1900, reached a decision 0.5244

Two identical variants, both converting at 10%, stopped the first time \(P(B > A)\) leaves a corridor. A single look at the end calls B better 5.03% of the time in this simulated setting using the approximate rule; the 0.95 threshold does not itself guarantee that long-run error rate. Looking every twenty users raises that to 27.96% — on data with no effect at all.

With a fixed prior, correct likelihood and ignorable observation-based stopping, an exact posterior can be updated at each look without an optional-stopping correction to that posterior. This code uses a normal approximation to the difference of Beta variables, not an exact posterior probability. Its error can matter with sparse successes or boundary rates. The reported simulation measures the approximate decision rule implemented here.

So the common formulation — “Bayesian methods let you peek” — is half right in a way that misleads. The exact posterior remains coherent under those conditions, while a numerical approximation needs its own accuracy check. What you may not do is assume that a rule built on those numbers has the long-run error rate its threshold suggests. Raising the threshold to 0.99 brings the rate to 8.38%, not to 1%.

The second block shows the other input. The same rule with a Beta(50, 450) prior on each arm — 500 prior observations centered on the truth — gives a wrong-call rate of 0.0851 against 0.2064 under a flat prior. The prior is part of the stopping rule, not context around it, because it decides how far a run of luck can push the posterior.

What you get in exchange

The last block runs the same machinery on a real 2-point lift. At a 0.95 threshold it calls B better 77.43% of the time with a median stop at 800 users; at 0.99, 51.98% with a median of 1,900. Earlier decisions, bought with more wrong ones — which is the same trade a group sequential design makes, expressed in a different vocabulary.

Three genuine advantages survive the caveat above. The output is directly interpretable — “a 95% probability that B is better” means what it sounds like, which a p-value does not. No sample size has to be fixed in advance, so a test can run until it resolves or funds run out. And the stopping rule does not enter the likelihood, under ignorable stopping and the stated model. This does not make a posterior mean or other stopped point estimate frequentist-unbiased.

Running one honestly

The practical upshot is that a Bayesian sequential design needs the same discipline as any other, applied to different objects. Fix the prior, the threshold, the look schedule and the maximum sample before the data, because all four determine the error rate. Then simulate the whole rule under the null to find out what that rate is — which is exactly what the first block above does, and it takes a few seconds.

If a calibrated false-positive rate is a requirement — a regulator, a published claim, a decision that has to be defensible — then measuring it is not optional and a threshold alone will not deliver it. If what you need is a coherent statement about which variant is better given everything observed, the posterior gives that at every look and the frequentist machinery never does.

The thresholds themselves are conventions. 0.95 carries no more authority than 0.05 does, and the number that belongs there comes from the cost of each mistake — see Bayesian decision theory and, for the rules in more detail, Bayesian stopping rules.

References

• SciPy documentation, scipy.stats.norm. docs.scipy.org