Sequential Probability Ratio Test (SPRT)

The sequential probability ratio test compares two simple hypotheses by tracking the likelihood ratio as observations arrive, and stopping the moment it leaves a corridor. Unlike a group sequential design it fixes neither the sample size nor the number of looks: it evaluates after every observation.

Wald's approximate bounds come from the two error rates directly — \(A = (1-\beta)/\alpha\) for accepting \(H_1\) and \(B = \beta/(1-\alpha)\) for accepting \(H_0\), with sampling continuing in between.

What the bounds deliver

import numpy as np
from scipy.stats import norm

alpha, beta_err, delta = 0.05, 0.20, 0.30
A, B = (1 - beta_err) / alpha, beta_err / (1 - alpha)
print(f"  Wald's bounds for alpha={alpha}, beta={beta_err}:  A = (1-b)/a = {A:.3f}"
      f"   B = b/(1-a) = {B:.4f}")
print(f"  in log units, stop for H1 above {np.log(A):.4f} and for H0 below {np.log(B):.4f}")
"""The test statistic is the running log likelihood ratio. Nothing else about
   the design is fixed - not the sample size, not the number of looks."""

rng = np.random.default_rng(0)
N, cap, chunk = 60_000, 1200, 4_000

def run(true_mean):
    """Simulated in chunks: the full matrix of running statistics would not fit."""
    ns, h1s, fin = [], [], []
    for _ in range(N // chunk):
        x = rng.normal(loc=true_mean, size=(chunk, cap))
        llr = delta * np.cumsum(x, axis=1) - np.arange(1, cap + 1) * delta ** 2 / 2
        up, done = llr >= np.log(A), (llr >= np.log(A)) | (llr <= np.log(B))
        n = np.where(done.any(axis=1), done.argmax(axis=1) + 1, cap)
        ns.append(n); h1s.append(up[np.arange(chunk), n - 1]); fin.append(done.any(axis=1))
    return np.concatenate(ns), np.concatenate(h1s), np.concatenate(fin)

fixed_n = int(np.ceil(((norm.ppf(1 - alpha) + norm.ppf(1 - beta_err)) / delta) ** 2))
print(f"\n  a fixed-horizon test of the same effect needs {fixed_n} observations")
print(f"\n{'truth':>22}{'chooses H1':>13}{'average n':>12}{'median n':>11}{'90th pct':>11}")
for label, mu in (("no effect (mu = 0)", 0.0), (f"the effect (mu = {delta})", delta)):
    n, h1, finished = run(mu)
    print(f"{label:>22}{h1.mean():>13.4f}{n.mean():>12.1f}{np.median(n):>11.0f}"
          f"{np.percentile(n, 90):>11.0f}")
print(f"  error rates land near the {alpha} and {beta_err} the bounds were built for,")
print(f"  and the average sample is well under the fixed test's {fixed_n}")

n, _, finished = run(0.0)
print(f"\n  but the sample size is a random variable. Under the null, capped at {cap}:")
print(f"    reached a decision {finished.mean():>8.4f} of the time")
for q in (50, 90, 99, 99.9):
    print(f"    {q:>5}th percentile of n {np.percentile(n, q):>8.0f}")
#   Wald's bounds for alpha=0.05, beta=0.2:  A = (1-b)/a = 16.000   B = b/(1-a) = 0.2105
#   in log units, stop for H1 above 2.7726 and for H0 below -1.5581
#
#   a fixed-horizon test of the same effect needs 69 observations
#
#                  truth   chooses H1   average n   median n   90th pct
#     no effect (mu = 0)       0.0440        34.3         25         73
#  the effect (mu = 0.3)       0.8310        47.9         39         91
#   error rates land near the 0.05 and 0.2 the bounds were built for,
#   and the average sample is well under the fixed test's 69
#
#   but the sample size is a random variable. Under the null, capped at 1200:
#     reached a decision   1.0000 of the time
#        50th percentile of n       24
#        90th percentile of n       72
#        99th percentile of n      145
#      99.9th percentile of n      222

With \(\alpha = 0.05\) and \(\beta = 0.20\) the bounds are 16 and 0.2105 — stop for the alternative when the data is sixteen times more likely under it, stop for the null when it is about five times less likely.

The measured type I error is 0.0440 and power is 0.8310 in this normal example. Both are conservative here, but the approximate bounds do not guarantee that each actual error rate is below its target. For H0 = Bernoulli(0.06) and H1 = Bernoulli(0.97), the same bounds 16 and 0.2105 stop after one observation: success gives LR = 16.1667 and failure gives 0.031915. The actual type I error is then 0.06, above 0.05. Calibrate the boundaries for the observation model when a strict error guarantee is required.

The efficiency is the headline. A fixed-horizon test of the same effect needs 69 observations. The SPRT averages 34.3 under the null and 47.9 under the alternative — roughly half and roughly seven tenths. This is not a free lunch extracted from nothing; it is the payoff for being allowed to stop, and for testing two specific hypotheses rather than one against everything else.

The sample size is a random variable

The last block is the constraint that governs whether you can use it. Under the null the median run ends at 24 observations, the 90th percentile at 72, the 99th at 145 and the 99.9th at 222 — nine times the median. Every run here finished within the cap, but in principle the test has no upper bound on its sample size.

That spread is a planning problem before it is a statistical one. A design whose duration is a random variable with a long tail is hard to schedule, hard to budget, and hard to explain to anyone waiting on the result. The standard response is to truncate at a maximum sample size — and truncation changes the error rates the bounds promised, so a truncated SPRT needs its boundaries recomputed rather than borrowed.

What it assumes

Two assumptions do most of the work, and both are stronger than they look.

• Two simple hypotheses. The test compares \(\theta = \theta_0\) against \(\theta = \theta_1\), both specified as points. Its optimality is for that comparison. If the true effect lies between them, the test still terminates but its error rates and expected sample size are no longer the designed ones.
• Independent, identically distributed observations arriving one at a time. The running statistic is a sum of per-observation log ratios, so correlated units, delayed outcomes, or a drifting population break it in the same way they break any likelihood-based method — quietly.

The point hypothesis requirement is the one that decides most applications. Choosing \(\theta_1\) is choosing the effect the test is tuned for, in the same way a fixed design's sample size is chosen for a target effect — the difference is that the SPRT's efficiency claim is specific to it, not merely its power.

Where the assumptions hold, the SPRT is the reason sequential methods exist at all: it is optimal in expected sample size among tests with the same error rates, for the two hypotheses it was built from. Where they do not, a group sequential design with a small number of planned looks gives up some of the efficiency in exchange for a bounded, schedulable study.

References

• SciPy documentation, scipy.stats.norm. docs.scipy.org
• NumPy documentation, numpy.random.Generator.normal. numpy.org