26. A Problem that Stumped Milton Friedman#

(and that Abraham Wald solved by inventing sequential analysis)

26.1. Overview#

This is the first of three lectures about a statistical decision problem that a US Navy Captain presented to Milton Friedman and W. Allen Wallis during World War II when they were analysts at the U.S. Government’s Statistical Research Group at Columbia University.

This problem led Abraham Wald [Wald, 1947] to formulate sequential analysis, an approach to statistical decision problems that is intimately related to dynamic programming.

In the spirit of this earlier lecture, the present lecture and its sequel approach the problem from two distinct points of view, one frequentist, the other Bayesian.

In this lecture, we describe Wald’s formulation of the problem from the perspective of a frequentist statistician working within the Neyman-Pearson tradition who thinks about testing hypotheses and consequently uses laws of large numbers to investigate limiting properties of particular statistics under a given hypothesis, i.e., a vector of parameters that pins down a particular member of a manifold of statistical models that interest the statistician.

In this related lecture we’ll discuss another formulation that adopts the perspective of a Bayesian statistician who views parameters as random variables that are jointly distributed with observable variables that he is concerned about.

Because we are taking a frequentist perspective that is concerned about relative frequencies conditioned on alternative parameter values, i.e., alternative hypotheses, key ideas in this lecture

  • Type I and type II statistical errors

    • a type I error occurs when you reject a null hypothesis that is true

    • a type II error occurs when you accept a null hypothesis that is false

  • The power of a frequentist statistical test

  • The size of a frequentist statistical test

  • The critical region of a statistical test

  • A most powerful test

  • The role of a Law of Large Numbers (LLN) in interpreting power and size of a frequentist statistical test

  • Abraham Wald’s sequential probability ratio test

We’ll begin with some imports:

import numpy as np
import matplotlib.pyplot as plt
from numba import njit, prange, vectorize, jit
from numba.experimental import jitclass
from math import gamma
from scipy.special import betaln, digamma
from scipy.stats import beta
from collections import namedtuple
import pandas as pd
import scipy as sp

This lecture uses ideas studied in the lecture on likelihood ratio processes and the lecture on Bayesian learning.

26.2. Source of the Problem#

On pages 137-139 of his 1998 book Two Lucky People with Rose Friedman [Friedman and Friedman, 1998], Milton Friedman described a problem presented to him and Allen Wallis during World War II, when they worked at the US Government’s Statistical Research Group at Columbia University.

Note

See pages 25 and 26 of Allen Wallis’s 1980 article [Wallis, 1980] about the Statistical Research Group at Columbia University during World War II for his account of the episode and for important contributions that Harold Hotelling made to formulating the problem. Also see chapter 5 of Jennifer Burns’ book about Milton Friedman [Burns, 2023].

Let’s listen to Milton Friedman tell us what happened

In order to understand the story, it is necessary to have an idea of a simple statistical problem, and of the standard procedure for dealing with it. The actual problem out of which sequential analysis grew will serve. The Navy has two alternative designs (say A and B) for a projectile. It wants to determine which is superior. To do so it undertakes a series of paired firings. On each round, it assigns the value 1 or 0 to A accordingly as its performance is superior or inferior to that of B and conversely 0 or 1 to B. The Navy asks the statistician how to conduct the test and how to analyze the results.

The standard statistical answer was to specify a number of firings (say 1,000) and a pair of percentages (e.g., 53% and 47%) and tell the client that if A receives a 1 in more than 53% of the firings, it can be regarded as superior; if it receives a 1 in fewer than 47%, B can be regarded as superior; if the percentage is between 47% and 53%, neither can be so regarded.

When Allen Wallis was discussing such a problem with (Navy) Captain Garret L. Schuyler, the captain objected that such a test, to quote from Allen’s account, may prove wasteful. If a wise and seasoned ordnance officer like Schuyler were on the premises, he would see after the first few thousand or even few hundred [rounds] that the experiment need not be completed either because the new method is obviously inferior or because it is obviously superior beyond what was hoped for \(\ldots\).

Friedman and Wallis worked on the problem for a while but didn’t completely solve it.

Realizing that, they told Abraham Wald about the problem.

That set Wald on a path that led him to create Sequential Analysis [Wald, 1947].

26.3. Neyman-Pearson formulation#

It is useful to begin by describing the theory underlying the test that the U.S. Navy told Captain Garret L. Schuyler to use.

Captain Schuyler’s doubts motivated him to tell Milton Friedman and Allen Wallis his conjecture that superior practical procedures existed.

Evidently, the Navy had told Captain Schuyler to use what was then a state-of-the-art Neyman-Pearson hypothesis test.

We’ll rely on Abraham Wald’s [Wald, 1947] elegant summary of Neyman-Pearson theory.

Watch for these features of the setup:

  • the assumption of a fixed sample size \(n\)

  • the application of laws of large numbers, conditioned on alternative probability models, to interpret probabilities \(\alpha\) and \(\beta\) of the type I and type II errors defined in the Neyman-Pearson theory

In chapter 1 of Sequential Analysis [Wald, 1947] Abraham Wald summarizes the Neyman-Pearson approach to hypothesis testing.

Wald frames the problem as making a decision about a probability distribution that is partially known.

(You have to assume that something is already known in order to state a well-posed problem – usually, something means a lot)

By limiting what is unknown, Wald uses the following simple structure to illustrate the main ideas:

  • A decision-maker wants to decide which of two distributions \(f_0\), \(f_1\) govern an IID random variable \(z\).

  • The null hypothesis \(H_0\) is the statement that \(f_0\) governs the data.

  • The alternative hypothesis \(H_1\) is the statement that \(f_1\) governs the data.

  • The problem is to devise and analyze a test of hypothesis \(H_0\) against the alternative hypothesis \(H_1\) on the basis of a sample of a fixed number \(n\) independent observations \(z_1, z_2, \ldots, z_n\) of the random variable \(z\).

To quote Abraham Wald,

A test procedure leading to the acceptance or rejection of the [null] hypothesis in question is simply a rule specifying, for each possible sample of size \(n\), whether the [null] hypothesis should be accepted or rejected on the basis of the sample. This may also be expressed as follows: A test procedure is simply a subdivision of the totality of all possible samples of size \(n\) into two mutually exclusive parts, say part 1 and part 2, together with the application of the rule that the [null] hypothesis be accepted if the observed sample is contained in part 2. Part 1 is also called the critical region. Since part 2 is the totality of all samples of size \(n\) which are not included in part 1, part 2 is uniquely determined by part 1. Thus, choosing a test procedure is equivalent to determining a critical region.

Let’s listen to Wald longer:

As a basis for choosing among critical regions the following considerations have been advanced by Neyman and Pearson: In accepting or rejecting \(H_0\) we may commit errors of two kinds. We commit an error of the first kind if we reject \(H_0\) when it is true; we commit an error of the second kind if we accept \(H_0\) when \(H_1\) is true. After a particular critical region \(W\) has been chosen, the probability of committing an error of the first kind, as well as the probability of committing an error of the second kind is uniquely determined. The probability of committing an error of the first kind is equal to the probability, determined by the assumption that \(H_0\) is true, that the observed sample will be included in the critical region \(W\). The probability of committing an error of the second kind is equal to the probability, determined on the assumption that \(H_1\) is true, that the observed sample will fall outside the critical region \(W\). For any given critical region \(W\) we shall denote the probability of an error of the first kind by \(\alpha\) and the probability of an error of the second kind by \(\beta\).

Let’s listen carefully to how Wald applies law of large numbers to interpret \(\alpha\) and \(\beta\):

The probabilities \(\alpha\) and \(\beta\) have the following important practical interpretation: Suppose that we draw a large number of samples of size \(n\). Let \(M\) be the number of such samples drawn. Suppose that for each of these \(M\) samples we reject \(H_0\) if the sample is included in \(W\) and accept \(H_0\) if the sample lies outside \(W\). In this way we make \(M\) statements of rejection or acceptance. Some of these statements will in general be wrong. If \(H_0\) is true and if \(M\) is large, the probability is nearly \(1\) (i.e., it is practically certain) that the proportion of wrong statements (i.e., the number of wrong statements divided by \(M\)) will be approximately \(\alpha\). If \(H_1\) is true, the probability is nearly \(1\) that the proportion of wrong statements will be approximately \(\beta\). Thus, we can say that in the long run [ here Wald applies law of large numbers by driving \(M \rightarrow \infty\) (our comment, not Wald’s) ] the proportion of wrong statements will be \(\alpha\) if \(H_0\) is true and \(\beta\) if \(H_1\) is true.

The quantity \(\alpha\) is called the size of the critical region, and the quantity \(1-\beta\) is called the power of the critical region.

Wald notes that

one critical region \(W\) is more desirable than another if it has smaller values of \(\alpha\) and \(\beta\). Although either \(\alpha\) or \(\beta\) can be made arbitrarily small by a proper choice of the critical region \(W\), it is impossible to make both \(\alpha\) and \(\beta\) arbitrarily small for a fixed value of \(n\), i.e., a fixed sample size.

Wald summarizes Neyman and Pearson’s setup as follows:

Neyman and Pearson show that a region consisting of all samples \((z_1, z_2, \ldots, z_n)\) which satisfy the inequality

\[ \frac{ f_1(z_1) \cdots f_1(z_n)}{f_0(z_1) \cdots f_0(z_n)} \geq k \]

is a most powerful critical region for testing the hypothesis \(H_0\) against the alternative hypothesis \(H_1\). The term \(k\) on the right side is a constant chosen so that the region will have the required size \(\alpha\).

Wald goes on to discuss Neyman and Pearson’s concept of uniformly most powerful test.

Our \(H_0\) and \(H_1\) each specify a single distribution, so the most powerful test above is all that we need.

Here is how Wald introduces the notion of a sequential test

A rule is given for making one of the following three decisions at any stage of the experiment (at the \(m\) th trial for each integral value of \(m\)): (1) to accept the hypothesis \(H\), (2) to reject the hypothesis \(H\), (3) to continue the experiment by making an additional observation. Thus, such a test procedure is carried out sequentially. On the basis of the first observation, one of the aforementioned decision is made. If the first or second decision is made, the process is terminated. If the third decision is made, a second trial is performed. Again, on the basis of the first two observations, one of the three decision is made. If the third decision is made, a third trial is performed, and so on. The process is continued until either the first or the second decisions is made. The number \(n\) of observations required by such a test procedure is a random variable, since the value of \(n\) depends on the outcome of the observations.

26.4. Wald’s sequential formulation#

By way of contrast to Neyman and Pearson’s formulation of the problem, in Wald’s formulation

  • The sample size \(n\) is not fixed but rather a random variable.

  • Two parameters \(A\) and \(B\) that are related to but distinct from Neyman and Pearson’s \(\alpha\) and \(\beta\); \(A\) and \(B\) characterize cut-off rules that Wald uses to determine the random variable \(n\) as a function of random outcomes.

Here is how Wald sets up the problem.

A decision-maker can observe a sequence of draws of a random variable \(z\).

He (or she) wants to know which of two probability distributions \(f_0\) or \(f_1\) governs \(z\).

We use beta distributions as examples.

We will also work with the Kullback–Leibler divergence introduced in Statistical Divergence Measures.

For two beta densities, Kullback–Leibler divergence has a closed form.

@vectorize
def p(x, a, b):
    """Beta distribution density function."""
    r = gamma(a + b) / (gamma(a) * gamma(b))
    return r * x** (a-1) * (1 - x) ** (b-1)

def create_beta_density(a, b):
    """Create a beta density function with specified parameters."""
    return jit(lambda x: p(x, a, b))

def compute_KL(a1, b1, a2, b2):
    """KL divergence KL(f, g) for f = Beta(a1, b1) and g = Beta(a2, b2)."""
    return (betaln(a2, b2) - betaln(a1, b1)
            + (a1 - a2) * digamma(a1) + (b1 - b2) * digamma(b1)
            + (a2 - a1 + b2 - b1) * digamma(a1 + b1))

The next figure shows two beta distributions

f0 = create_beta_density(1, 1)
f1 = create_beta_density(9, 9)
grid = np.linspace(0, 1, 50)

fig, ax = plt.subplots()
ax.plot(grid, f0(grid), lw=2, label="$f_0$")
ax.plot(grid, f1(grid), lw=2, label="$f_1$")
ax.legend()
ax.set(xlabel="$z$ values", ylabel="density")
plt.tight_layout()
plt.show()
_images/e07a6db9bddbf33e2dcc3aab6d247ad12aab8b3fe947f417dcb0d831fc9836c0.png

Conditional on knowing that successive observations are drawn from distribution \(f_0\), the sequence of random variables is independently and identically distributed (IID).

Conditional on knowing that successive observations are drawn from distribution \(f_1\), the sequence of random variables is also independently and identically distributed (IID).

But the observer does not know which of the two distributions generated the sequence.

For reasons explained in Exchangeability and Bayesian Updating, this means that the observer thinks that the sequence is not IID.

Consequently, the observer has something to learn, namely, whether the observations are drawn from \(f_0\) or from \(f_1\).

The decision maker wants to decide which of the two distributions is generating outcomes.

26.4.1. Type I and type II errors#

If we regard \(f=f_0\) as a null hypothesis and \(f=f_1\) as an alternative hypothesis, then

  • a type I error is an incorrect rejection of a true null hypothesis (a “false positive”)

  • a type II error is a failure to reject a false null hypothesis (a “false negative”)

To repeat ourselves

  • \(\alpha\) is the probability of a type I error

  • \(\beta\) is the probability of a type II error

26.4.2. Choices#

After observing \(z_k, z_{k-1}, \ldots, z_1\), the decision-maker chooses among three distinct actions:

  • He decides that \(f = f_0\) and draws no more \(z\)’s

  • He decides that \(f = f_1\) and draws no more \(z\)’s

  • He postpones deciding and instead chooses to draw \(z_{k+1}\)

Wald defines

  • \(p_{0m} = f_0(z_1) \cdots f_0(z_m)\)

  • \(p_{1m} = f_1(z_1) \cdots f_1(z_m)\)

  • \(L_{m} = \frac{p_{1m}}{p_{0m}}\)

for \(m \geq 1\) and sets \(L_0 = 1\).

Here \(\{L_m\}_{m=0}^\infty\) is a likelihood ratio process.

Wald’s sequential decision rule is parameterized by real numbers \(0 < B < 1 < A\).

For a given pair \(A, B\), the decision rule is

\[\begin{split} \begin{aligned} \textrm { accept } f=f_1 \textrm{ if } L_m \geq A \\ \textrm { accept } f=f_0 \textrm{ if } L_m \leq B \\ \textrm { draw another } z \textrm{ if } B < L_m < A \end{aligned} \end{split}\]

Note

The sequel A Bayesian Formulation of Friedman and Wald’s Problem uses \(A\) and \(B\) for thresholds on the posterior probability \(\pi_k\) that \(f = f_1\), and \(L_0\) and \(L_1\) for losses from wrong decisions.

There we write \(\Lambda_k\) for the likelihood ratio and \(A_\Lambda\), \(B_\Lambda\) for the thresholds on it that correspond to Wald’s \(A\) and \(B\) here.

The following figure illustrates aspects of Wald’s procedure.

_images/wald_dec_rule.png

26.6. Simulations#

We experiment with different distributions \(f_0\) and \(f_1\) to examine how Wald’s test performs under various conditions.

Our goal in conducting these simulations is to understand trade-offs between decision speed and accuracy associated with Wald’s sequential probability ratio test.

Specifically, we will watch how:

  • The decision thresholds \(A\) and \(B\) (or equivalently the target error rates \(\alpha\) and \(\beta\)) affect the average stopping time

  • The discrepancy between distributions \(f_0\) and \(f_1\) affects average stopping times

We will focus on the case where \(f_0\) and \(f_1\) are beta distributions since it is easy to control the overlapping regions of the two densities by adjusting their shape parameters.

First, we define a namedtuple to store all the parameters we need for our simulation studies.

We also compute Wald’s recommended thresholds \(A\) and \(B\) based on the target type I and type II errors \(\alpha\) and \(\beta\)

SPRTParams = namedtuple('SPRTParams', 
                ['α', 'β',  # Target type I and type II errors
                'a0', 'b0', # Shape parameters for f_0
                'a1', 'b1', # Shape parameters for f_1
                'N',        # Number of simulations
                'seed'])

@njit
def compute_wald_thresholds(α, β):
    """Compute Wald's recommended thresholds."""
    A = (1 - β) / α
    B = β / (1 - α)
    return A, B, np.log(A), np.log(B)

Now we can run the simulation following Wald’s recommendation.

We’ll compare the log-likelihood ratio to logarithms of the thresholds \(\log(A)\) and \(\log(B)\).

The following algorithm underlies our simulations.

  1. Compute thresholds \(A = \frac{1-\beta}{\alpha}\), \(B = \frac{\beta}{1-\alpha}\) and work with \(\log A\), \(\log B\).

  2. Given true distribution (either \(f_0\) or \(f_1\)):

    • Initialize log-likelihood ratio \(\log L_0 = 0\)

    • Repeat:

      • Draw observation \(z\) from the true distribution

      • Update: \(\log L_{n+1} \leftarrow \log L_n + (\log f_1(z) - \log f_0(z))\)

      • If \(\log L_{n+1} \geq \log A\): stop, reject \(H_0\)

      • If \(\log L_{n+1} \leq \log B\): stop, accept \(H_0\)

  3. Repeat step 2 for \(N\) replications with \(N/2\) replications for each distribution, compute the empirical type I error \(\hat{\alpha}\) and type II error \(\hat{\beta}\) with

\[ \hat{\alpha} = \frac{\text{$\#$ of times reject } H_0 \text{ when } f_0 \text{ is true}}{\text{$\#$ of replications with } f_0 \text{ true}} \]
\[ \hat{\beta} = \frac{\text{$\#$ of times accept } H_0 \text{ when } f_1 \text{ is true}}{\text{$\#$ of replications with } f_1 \text{ true}} \]
@njit
def sprt_single_run(a0, b0, a1, b1, logA, logB, true_f0, seed):
    """Run a single SPRT until a decision is reached."""
    log_L = 0.0
    n = 0
    np.random.seed(seed)
    
    while True:
        z = np.random.beta(a0, b0) if true_f0 else np.random.beta(a1, b1)
        n += 1
        
        # Update log-likelihood ratio
        log_L += np.log(p(z, a1, b1)) - np.log(p(z, a0, b0))
        
        # Check stopping conditions
        if log_L >= logA:
            return n, False  # Reject H0
        elif log_L <= logB:
            return n, True   # Accept H0

@njit(parallel=True)
def run_sprt_simulation(a0, b0, a1, b1, α, β, N, seed):
    """SPRT simulation."""
    A, B, logA, logB = compute_wald_thresholds(α, β)
    
    stopping_times = np.zeros(N, dtype=np.int64)
    decisions_h0 = np.zeros(N, dtype=np.bool_)
    truth_h0 = np.zeros(N, dtype=np.bool_)
    
    for i in prange(N):
        true_f0 = (i % 2 == 0)
        truth_h0[i] = true_f0
        
        n, accept_f0 = sprt_single_run(
                        a0, b0, a1, b1, 
                        logA, logB, 
                        true_f0, seed + i)
        stopping_times[i] = n
        decisions_h0[i] = accept_f0
    
    return stopping_times, decisions_h0, truth_h0

def run_sprt(params):
    """Run SPRT simulations with given parameters."""
    stopping_times, decisions_h0, truth_h0 = run_sprt_simulation(
        params.a0, params.b0, params.a1, params.b1, 
        params.α, params.β, params.N, params.seed
    )
    
    # Calculate error rates
    truth_h0_bool = truth_h0.astype(bool)
    decisions_h0_bool = decisions_h0.astype(bool)
    
    type_I = np.sum(truth_h0_bool & ~decisions_h0_bool) \
            / np.sum(truth_h0_bool)
    type_II = np.sum(~truth_h0_bool & decisions_h0_bool) \
            / np.sum(~truth_h0_bool)
    
    return {
        'stopping_times': stopping_times,
        'decisions_h0': decisions_h0_bool,
        'truth_h0': truth_h0_bool,
        'type_I': type_I,
        'type_II': type_II
    }

# Run simulation
params = SPRTParams(α=0.05, β=0.10, a0=2, b0=5, a1=5, b1=2, N=20000, seed=1)
results = run_sprt(params)

print(f"Average stopping time: {results['stopping_times'].mean():.2f}")
print(f"Empirical type I  error: {results['type_I']:.3f} (target = {params.α})")
print(f"Empirical type II error: {results['type_II']:.3f} (target = {params.β})")

Average stopping time: 1.59
Empirical type I  error: 0.012 (target = 0.05)
Empirical type II error: 0.022 (target = 0.1)

As anticipated in the passage above in which Wald discussed the quality of \(a(\alpha, \beta), b(\alpha, \beta)\) given in approximation (26.1), we find that the algorithm actually gives lower type I and type II error rates than the target values.

The empirical error rates satisfy Wald’s bounds (26.2), which here are \(\alpha' \leq 0.05/0.90 \approx 0.056\), \(\beta' \leq 0.10/0.95 \approx 0.105\), and \(\alpha' + \beta' \leq 0.15\).

The bounds are far from tight because \(f_0 = \text{Beta}(2,5)\) and \(f_1 = \text{Beta}(5,2)\) are far apart.

A single draw often carries \(\log L_m\) far past \(\log A\) or \(\log B\), and that overshoot pushes \(\alpha'\) and \(\beta'\) below \(\alpha\) and \(\beta\).

Note

For recent work on the quality of approximation (26.1), see, e.g., [Fischer and Ramdas, 2024].

The following code creates a few graphs that illustrate the results of our simulation.

Hide code cell source

 
def plot_sprt_results(results, params, title=""):
    """Plot SPRT results."""
    fig, axes = plt.subplots(1, 3, figsize=(20, 6))
    
    # Distribution plots
    z_grid = np.linspace(0, 1, 200)
    f0 = create_beta_density(params.a0, params.b0)
    f1 = create_beta_density(params.a1, params.b1)
    
    axes[0].plot(z_grid, f0(z_grid), 'b-', lw=2, 
                 label=f'$f_0 = \\text{{Beta}}({params.a0},{params.b0})$')
    axes[0].plot(z_grid, f1(z_grid), 'r-', lw=2, 
                 label=f'$f_1 = \\text{{Beta}}({params.a1},{params.b1})$')
    axes[0].fill_between(z_grid, 0, 
                        np.minimum(f0(z_grid), f1(z_grid)), 
                        alpha=0.3, color='purple', label='overlap')
    if title:
        axes[0].set_title(title, fontsize=20)
    axes[0].set_xlabel('z', fontsize=16)
    axes[0].set_ylabel('density', fontsize=16)
    axes[0].legend(fontsize=14)
    
    # Stopping times
    max_n = min(results['stopping_times'].max(), 101)
    bins = np.arange(1, max_n + 2) - 0.5
    axes[1].hist(results['stopping_times'], bins=bins, 
                 color="steelblue", alpha=0.8, edgecolor="black")
    axes[1].set_title(f'stopping times (μ={results["stopping_times"].mean():.1f})', 
                      fontsize=16)
    axes[1].set_xlabel('n', fontsize=16)
    axes[1].set_ylabel('frequency', fontsize=16)
    axes[1].set_xlim(0, 100)
    
    # Confusion matrix
    plot_confusion_matrix(results, axes[2])
    
    plt.tight_layout()
    plt.show()

def plot_confusion_matrix(results, ax):
    """Plot confusion matrix for SPRT results."""
    f0_correct = np.sum(results['truth_h0'] & results['decisions_h0'])
    f0_incorrect = np.sum(results['truth_h0'] & (~results['decisions_h0']))
    f1_correct = np.sum((~results['truth_h0']) & (~results['decisions_h0']))
    f1_incorrect = np.sum((~results['truth_h0']) & results['decisions_h0'])
    
    confusion_data = np.array([[f0_correct, f0_incorrect], 
                              [f1_incorrect, f1_correct]])
    row_totals = confusion_data.sum(axis=1, keepdims=True)
    
    im = ax.imshow(confusion_data, cmap='Blues', aspect='equal')
    ax.set_title(f'errors: I={results["type_I"]:.3f} II={results["type_II"]:.3f}', 
                 fontsize=16)
    ax.set_xticks([0, 1])
    ax.set_xticklabels(['accept $H_0$', 'reject $H_0$'], fontsize=14)
    ax.set_yticks([0, 1])
    ax.set_yticklabels(['true $f_0$', 'true $f_1$'], fontsize=14)
    
    for i in range(2):
        for j in range(2):
            percent = confusion_data[i, j] / row_totals[i, 0] \
                        if row_totals[i, 0] > 0 else 0
            color = 'white' if confusion_data[i, j] > confusion_data.max() * 0.5 \
                    else 'black'
            ax.text(j, i, f'{confusion_data[i, j]}\n({percent:.1%})',
                   ha="center", va="center", color=color, fontweight='bold', 
                   fontsize=14)

Let’s plot the results of our simulation

plot_sprt_results(results, params)
_images/6aed83aad1256b4671b776f987eb08149e7c1dce6ffb86a9d77a7d614dff6219.png

In this example, the stopping time never exceeds 10, and the mean stopping time is 1.59.

We can construct a \(2 \times 2\) “confusion matrix” whose diagonal elements count the number of times that Wald’s decision rule correctly accepts and rejects the null hypothesis.

print("Confusion Matrix data:")
print(f"Type I error: {results['type_I']:.3f}")
print(f"Type II error: {results['type_II']:.3f}")

Confusion Matrix data:
Type I error: 0.012
Type II error: 0.022

Next we use our code to study three different \(f_0, f_1\) pairs having different discrepancies between distributions.

We plot the same three graphs we used above for each pair of distributions

params_1 = SPRTParams(α=0.05, β=0.10, a0=2, b0=8, a1=8, b1=2, N=5000, seed=42)
results_1 = run_sprt(params_1)

params_2 = SPRTParams(α=0.05, β=0.10, a0=4, b0=5, a1=5, b1=4, N=5000, seed=42)
results_2 = run_sprt(params_2)

params_3 = SPRTParams(α=0.05, β=0.10, a0=0.5, b0=0.4, a1=0.4, 
                      b1=0.5, N=5000, seed=42)
results_3 = run_sprt(params_3)

plot_sprt_results(results_1, params_1)
_images/313651808e40965fa972f41993551a21956fa1cd849aa472fd1e8dc9abeea5fc.png 
plot_sprt_results(results_2, params_2)
_images/e458dbcdf981246020a9e99f70854e9ada769a25c14ebc298fcb804cf6fd85b3.png 
plot_sprt_results(results_3, params_3)
_images/498497e010d3e4900792b51565ebc584cfadd76353818a1d6824604619957da1.png

Notice that the stopping times are less when the two distributions are farther apart.

This makes sense.

When two distributions are “far apart”, it should not take too long to decide which one is generating the data.

When two distributions are “close”, it should take longer to decide which one is generating the data.

The Kullback–Leibler divergence described in Statistical Divergence Measures measures how far apart two distributions are.

Wald linked it to the expected number of draws that his test requires.

Under \(H_0\), \(\log L_N = \sum_{i=1}^N \log \frac{f_1(z_i)}{f_0(z_i)}\) is a sum of a random number \(N\) of IID terms, each with mean \(E_0 \log \frac{f_1(z)}{f_0(z)} = -KL(f_0, f_1)\).

Under \(H_1\), each term has mean \(E_1 \log \frac{f_1(z)}{f_0(z)} = KL(f_1, f_0)\).

Because \(N\) is a stopping time with finite mean, Wald’s identity gives

(26.3)#\[ E_0 [\log L_N] = - E_0[N] \, KL(f_0, f_1), \qquad E_1 [\log L_N] = E_1[N] \, KL(f_1, f_0) \]

This section of Likelihood Ratio Processes describes the link between KL divergence and the expected log likelihood ratio for a fixed number of draws.

Wald ignored overshoot: he approximated \(\log L_N\) by \(\log A\) when the test rejects \(H_0\) and by \(\log B\) when it accepts \(H_0\).

Setting \(P_0\{L_N \geq A\} \approx \alpha\) and \(P_1\{L_N \leq B\} \approx \beta\) in (26.3) then gives Wald’s approximations to the average sample numbers

(26.4)#\[ E_0[N] \approx \frac{(1-\alpha)\log(1/B) - \alpha \log A}{KL(f_0, f_1)}, \qquad E_1[N] \approx \frac{(1-\beta)\log A - \beta \log(1/B)}{KL(f_1, f_0)} \]

KL divergence is not symmetric, so (26.4) divides by \(KL(f_0, f_1)\) under \(H_0\) and by \(KL(f_1, f_0)\) under \(H_1\).

In the experiments below, \(f_0 = \text{Beta}(a, b)\) and \(f_1 = \text{Beta}(b, a)\), so \(f_1(z) = f_0(1-z)\) and \(KL(f_0, f_1) = KL(f_1, f_0)\).

Half of our replications draw from \(f_0\) and half from \(f_1\), so (26.4) makes the mean stopping time proportional to \(1/KL(f_0, f_1)\).

We compute \(KL(f_0, f_1)\) and the mean stopping time for 100 such pairs and compare the mean stopping time with (26.4).

def generate_β_pairs(N=100, T=10.0, d_min=0.5, d_max=9.5):
    ds = np.linspace(d_min, d_max, N)
    a0 = (T - ds) / 2
    b0 = (T + ds) / 2
    return list(zip(a0, b0, b0, a0))

def wald_asn(α, β, KL_01, KL_10):
    """Wald's approximations to E_0[N] and E_1[N]."""
    A, B, logA, logB = compute_wald_thresholds(α, β)
    E0_N = ((1 - α) * (-logB) - α * logA) / KL_01
    E1_N = ((1 - β) * logA - β * (-logB)) / KL_10
    return E0_N, E1_N

param_comb = generate_β_pairs()

KL_divs = []
mean_stopping_times = []
asn_approx = []
param_list = []

for a0, b0, a1, b1 in param_comb:
    KL_01 = compute_KL(a0, b0, a1, b1)
    KL_10 = compute_KL(a1, b1, a0, b0)
    
    params_d = SPRTParams(α=0.05, β=0.10, a0=a0, b0=b0, 
                          a1=a1, b1=b1, N=5000, seed=42)
    results_d = run_sprt(params_d)
    E0_N, E1_N = wald_asn(params_d.α, params_d.β, KL_01, KL_10)
    
    KL_divs.append(KL_01)
    mean_stopping_times.append(results_d['stopping_times'].mean())
    asn_approx.append(0.5 * (E0_N + E1_N))
    param_list.append((a0, b0, a1, b1))

KL_divs = np.array(KL_divs)
mean_stopping_times = np.array(mean_stopping_times)
asn_approx = np.array(asn_approx)

fig, ax = plt.subplots()
ax.scatter(KL_divs, mean_stopping_times, s=40, alpha=0.7, 
           label='simulated mean stopping time')
ax.plot(KL_divs, asn_approx, color='C1', lw=2, 
        label="Wald's approximation")
ax.axhline(1, color='black', linestyle='--', lw=1, label='$N = 1$')
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_xlabel('$KL(f_0, f_1)$', fontsize=14)
ax.set_ylabel('mean stopping time', fontsize=14)
ax.legend()
plt.tight_layout()
plt.show()

ratio = mean_stopping_times / asn_approx
print(f"corr(mean stopping time, 1/KL) = "
      f"{np.corrcoef(mean_stopping_times, 1 / KL_divs)[0, 1]:.4f}")
print(f"simulated / approximate mean stopping time: "
      f"{ratio.min():.2f} (closest pair) to {ratio.max():.1f} (farthest pair)")
_images/895d6a89a5cd8592901cc5dd00b662a0ed313789440860dfd56828edf8606723.png 
corr(mean stopping time, 1/KL) = 0.9999
simulated / approximate mean stopping time: 1.14 (closest pair) to 28.1 (farthest pair)

On log scales, Wald’s approximation is a straight line with slope \(-1\).

Across the 100 pairs, the correlation between the mean stopping time and \(1/KL(f_0, f_1)\) is 0.9999.

The mean stopping time is close to a linear function of \(1/KL(f_0, f_1)\).

The simulated mean stopping times lie above Wald’s approximation for two reasons.

Overshoot makes \(|\log L_N|\) larger than \(\log A\) or \(\log(1/B)\), so by (26.3) it makes \(E_0[N]\) and \(E_1[N]\) larger than (26.4).

The test always takes at least one draw, so the mean stopping time cannot fall below 1, while (26.4) goes to zero as \(KL(f_0, f_1)\) grows.

Exercise 26.4 asks you to study these two sources of discrepancy.

Below are sampled examples from the experiments we have above

def plot_beta_distributions_grid(param_list, KL_divs, mean_stopping_times, 
                                selected_indices=None):
    """Plot grid of beta distributions with KL divergence and stopping times."""
    if selected_indices is None:
        selected_indices = [0, len(param_list)//6, len(param_list)//3, 
                          len(param_list)//2, 2*len(param_list)//3, -1]
    
    fig, axes = plt.subplots(2, 3, figsize=(15, 8))
    z_grid = np.linspace(0, 1, 200)
    
    for i, idx in enumerate(selected_indices):
        row, col = i // 3, i % 3
        a0, b0, a1, b1 = param_list[idx]
        
        f0 = create_beta_density(a0, b0)
        f1 = create_beta_density(a1, b1)
        
        axes[row, col].plot(z_grid, f0(z_grid), 'b-', lw=2, label='$f_0$')
        axes[row, col].plot(z_grid, f1(z_grid), 'r-', lw=2, label='$f_1$')
        axes[row, col].fill_between(z_grid, 0, 
                                  np.minimum(f0(z_grid), f1(z_grid)), 
                                  alpha=0.3, color='purple')
        
        axes[row, col].set_title(f'KL divergence: {KL_divs[idx]:.3f}'
                               f'\nMean time: {mean_stopping_times[idx]:.1f}', 
                               fontsize=12)
        axes[row, col].set_xlabel('z', fontsize=10)
        if i == 0:
            axes[row, col].set_ylabel('density', fontsize=10)
            axes[row, col].legend(fontsize=10)

    plt.tight_layout()
    plt.show()

plot_beta_distributions_grid(param_list, KL_divs, mean_stopping_times)
_images/4ab33396683ce338f4679e32ca2e49ccd6ec09392e66cb41407deb0c891f0a8c.png

Again, we find that the stopping time is shorter when the distributions are more separated, as measured by KL divergence.

Let’s visualize individual likelihood ratio processes to see how they evolve toward the decision boundaries.

def plot_likelihood_paths(params, n_highlight=10, n_background=200):
    """visualize likelihood ratio paths."""
    A, B, logA, logB = compute_wald_thresholds(params.α, params.β)
    f0, f1 = map(lambda ab: create_beta_density(*ab),
             [(params.a0, params.b0), 
              (params.a1, params.b1)])
    
    fig, axes = plt.subplots(1, 2, figsize=(14, 7))
    
    for dist_idx, (true_f0, ax, title) in enumerate([
        (True, axes[0], 'true distribution: $f_0$'),
        (False, axes[1], 'true distribution: $f_1$')
    ]):
        rng = np.random.default_rng(seed=42 + dist_idx)
        paths_data = []
        
        # Generate paths
        for path in range(n_background + n_highlight):
            log_L_path, log_L, n = [0.0], 0.0, 0
            
            while True:
                z = rng.beta(params.a0, params.b0) if true_f0 \
                    else rng.beta(params.a1, params.b1)
                n += 1
                log_L += np.log(f1(z)) - np.log(f0(z))
                log_L_path.append(log_L)
                
                if log_L >= logA or log_L <= logB:
                    paths_data.append((log_L_path, n, log_L >= logA))
                    break
        
        # Plot background paths
        for path, _, decision in paths_data[:n_background]:
            ax.plot(range(len(path)), path, color='C1' if decision else 'C0', 
                   alpha=0.2, linewidth=0.5)
        
        # Plot highlighted paths with labels
        for i, (path, _, decision) in enumerate(paths_data[n_background:]):
            ax.plot(range(len(path)), path, color='C1' if decision else 'C0', 
                   alpha=0.8, linewidth=1.5,
                   label='reject $H_0$' if decision and i == 0 else (
                         'accept $H_0$' if not decision and i == 0 else ''))
        
        # Add threshold lines and formatting
        ax.axhline(y=logA, color='C1', linestyle='--', linewidth=2, 
                  label=f'$\\log A = {logA:.2f}$')
        ax.axhline(y=logB, color='C0', linestyle='--', linewidth=2, 
                  label=f'$\\log B = {logB:.2f}$')
        ax.axhline(y=0, color='black', linestyle='-', alpha=0.5, linewidth=1)
        
        ax.set_xlabel(r'$n$') 
        ax.set_ylabel(r'$\log(L_n)$')
        ax.set_title(title, fontsize=20)
        ax.legend(fontsize=18, loc='center right')
        
        y_margin = max(abs(logA), abs(logB)) * 0.2
        ax.set_ylim(logB - y_margin, logA + y_margin)
    
    plt.tight_layout()
    plt.show()

plot_likelihood_paths(params_3, n_highlight=10, n_background=100)
_images/dd6aca7b6a7daf50e118209c4d91996052b88547fd8dd320c70311a78d3bf6a9.png

Next, let’s adjust the decision thresholds \(A\) and \(B\) and examine how the mean stopping time and the type I and type II error rates change.

In the code below, we adjust Wald’s rule by adjusting the thresholds \(A\) and \(B\) using factors \(A_f\) and \(B_f\).

@njit(parallel=True)  
def run_adjusted_thresholds(a0, b0, a1, b1, α, β, N, seed, A_f, B_f):
    """SPRT simulation with adjusted thresholds."""
    
    # Calculate original thresholds  
    A_original = (1 - β) / α
    B_original = β / (1 - α)
    
    # Apply adjustment factors
    A_adj = A_original * A_f
    B_adj = B_original * B_f
    logA = np.log(A_adj)
    logB = np.log(B_adj)
    
    # Pre-allocate arrays
    stopping_times = np.zeros(N, dtype=np.int64)
    decisions_h0 = np.zeros(N, dtype=np.bool_)
    truth_h0 = np.zeros(N, dtype=np.bool_)
    
    # Run simulations in parallel
    for i in prange(N):
        true_f0 = (i % 2 == 0)
        truth_h0[i] = true_f0
        
        n, accept_f0 = sprt_single_run(a0, b0, a1, b1, 
                        logA, logB, true_f0, seed + i)
        stopping_times[i] = n
        decisions_h0[i] = accept_f0
    
    return stopping_times, decisions_h0, truth_h0, A_adj, B_adj

def run_adjusted(params, A_f=1.0, B_f=1.0):
    """Wrapper to run SPRT with adjusted A and B thresholds."""
    
    stopping_times, decisions_h0, truth_h0, A_adj, B_adj = run_adjusted_thresholds(
        params.a0, params.b0, params.a1, params.b1, 
        params.α, params.β, params.N, params.seed, A_f, B_f
    )
    truth_h0_bool = truth_h0.astype(bool)
    decisions_h0_bool = decisions_h0.astype(bool)
    
    # Calculate error rates
    type_I = np.sum(truth_h0_bool 
                    & ~decisions_h0_bool) / np.sum(truth_h0_bool)
    type_II = np.sum(~truth_h0_bool 
                    & decisions_h0_bool) / np.sum(~truth_h0_bool)
    
    return {
        'stopping_times': stopping_times,
        'type_I': type_I,
        'type_II': type_II,
        'A_used': A_adj,
        'B_used': B_adj
    }

adjustments = [
    (5.0, 0.5), 
    (1.0, 1.0),    
    (0.3, 3.0),    
    (0.2, 5.0),    
    (0.15, 7.0),   
]

results_table = []
for A_f, B_f in adjustments:
    result = run_adjusted(params_2, A_f, B_f)
    results_table.append([
        A_f, B_f, 
        f"{result['stopping_times'].mean():.1f}",
        f"{result['type_I']:.3f}",
        f"{result['type_II']:.3f}"
    ])

df = pd.DataFrame(results_table, 
                 columns=["A_f", "B_f", "mean stop time", 
                          "Type I error", "Type II error"])
df = df.set_index(["A_f", "B_f"])
df
mean stop time Type I error Type II error
A_f B_f
5.00 0.5 16.1 0.006 0.036
1.00 1.0 11.1 0.033 0.070
0.30 3.0 5.5 0.086 0.195
0.20 5.0 3.4 0.120 0.304
0.15 7.0 2.2 0.146 0.410

Let’s pause and think about the table more carefully by referring back to (26.1).

Recall that \(A = \frac{1-\beta}{\alpha}\) and \(B = \frac{\beta}{1-\alpha}\).

When we multiply \(A\) by a factor less than 1 (making \(A\) smaller), we are effectively making it easier to reject the null hypothesis \(H_0\).

This increases the probability of Type I errors.

When we multiply \(B\) by a factor greater than 1 (making \(B\) larger), we are making it easier to accept the null hypothesis \(H_0\).

This increases the probability of Type II errors.

The table confirms this intuition: as \(A\) decreases and \(B\) increases from Wald’s approximate values (26.1), both Type I and Type II error rates increase, while the mean stopping time decreases.

26.8. Exercises#

In the first two exercises below, please try to rewrite the entire SPRT suite in this lecture.

Exercise 26.1

In the first exercise, we apply the sequential probability ratio test to distinguish two models generated by 3-state Markov chains

(For a review on likelihood ratio processes for Markov chains, see this section.)

Consider distinguishing between two 3-state Markov chain models using Wald’s sequential probability ratio test.

You have competing hypotheses about the transition probabilities:

  • \(H_0\): The chain follows transition matrix \(P^{(0)}\)

  • \(H_1\): The chain follows transition matrix \(P^{(1)}\)

Given transition matrices:

\[\begin{split} P^{(0)} = \begin{bmatrix} 0.7 & 0.2 & 0.1 \\ 0.3 & 0.5 & 0.2 \\ 0.1 & 0.3 & 0.6 \end{bmatrix}, \quad P^{(1)} = \begin{bmatrix} 0.5 & 0.3 & 0.2 \\ 0.2 & 0.6 & 0.2 \\ 0.2 & 0.2 & 0.6 \end{bmatrix} \end{split}\]

For a sequence of observations \((x_0, x_1, \ldots, x_t)\), the likelihood ratio is:

\[ \Lambda_t = \frac{\pi_{x_0}^{(1)}}{\pi_{x_0}^{(0)}} \prod_{s=1}^t \frac{P_{x_{s-1},x_s}^{(1)}}{P_{x_{s-1},x_s}^{(0)}} \]

where \(\pi^{(i)}\) is the stationary distribution under hypothesis \(i\).

Tasks:

  1. Implement the likelihood ratio computation for Markov chains

  2. Implement Wald’s sequential test with Type I error \(\alpha = 0.05\) and Type II error \(\beta = 0.10\)

  3. Run 1000 simulations under each hypothesis and compute empirical error rates

  4. Analyze the distribution of stopping times

The test stops when:

  • \(\Lambda_t \geq A = \frac{1-\beta}{\alpha} = 18\): Reject \(H_0\)

  • \(\Lambda_t \leq B = \frac{\beta}{1-\alpha} = 0.105\): Accept \(H_0\)

Exercise 26.2

In this exercise, apply Wald’s sequential test to distinguish between two VAR(1) models with different dynamics and noise structures.

For a review of the likelihood ratio process with VAR models, see Likelihood Processes For VAR Models.

Given VAR models under each hypothesis:

  • \(H_0\): \(x_{t+1} = A^{(0)} x_t + C^{(0)} w_{t+1}\)

  • \(H_1\): \(x_{t+1} = A^{(1)} x_t + C^{(1)} w_{t+1}\)

where \(w_t \sim \mathcal{N}(0, I)\) and:

\[\begin{split} A^{(0)} = \begin{bmatrix} 0.8 & 0.1 \\ 0.2 & 0.7 \end{bmatrix}, \quad C^{(0)} = \begin{bmatrix} 0.3 & 0.1 \\ 0.1 & 0.3 \end{bmatrix} \end{split}\]
\[\begin{split} A^{(1)} = \begin{bmatrix} 0.6 & 0.2 \\ 0.3 & 0.5 \end{bmatrix}, \quad C^{(1)} = \begin{bmatrix} 0.4 & 0 \\ 0 & 0.4 \end{bmatrix} \end{split}\]

Tasks:

  1. Implement the VAR likelihood ratio using the functions from the VAR lecture

  2. Implement Wald’s sequential test with \(\alpha = 0.05\) and \(\beta = 0.10\)

  3. Analyze performance under both hypotheses and with model misspecification, when a third VAR(1) model generates the data

  4. Compare with the Markov chain case in terms of stopping times and accuracy

Exercise 26.3

Wald and Wolfowitz [1948] proved that among all tests, sequential or not, whose type I and type II error probabilities do not exceed those of an SPRT, the SPRT minimizes both \(E_0[N]\) and \(E_1[N]\).

This exercise compares the SPRT with the fixed-sample Neyman–Pearson test that the Navy told Captain Schuyler to use.

Use results_2, which we computed above for \(f_0 = \text{Beta}(4, 5)\) and \(f_1 = \text{Beta}(5, 4)\).

  1. Compute the empirical error probabilities \(\hat \alpha\) and \(\hat \beta\) of the SPRT and its mean stopping times under \(f_0\) and under \(f_1\).

  2. For each \(n = 1, \ldots, 40\), construct the most powerful test of size \(\hat \alpha\) based on a fixed sample of \(n\) draws, and compute its type II error probability \(\beta(n)\).

  3. Find the smallest \(n\) with \(\beta(n) \leq \hat \beta\) and compare it with the SPRT’s mean stopping times.

Hint: for these two densities, \(\log \left( f_1(z)/f_0(z) \right) = \log \left( z/(1-z) \right)\).

Exercise 26.4

Wald’s approximation (26.4) ignores overshoot and the fact that \(N \geq 1\).

  1. Modify sprt_single_run so that it also returns \(\log L_N\).

  2. For every fifth beta pair generated by generate_β_pairs, simulate 50,000 replications under \(f_0\) and 50,000 under \(f_1\) and estimate \(E_0[N]\), \(E_1[N]\), \(E_0[\log L_N]\), and \(E_1[\log L_N]\).

  3. Verify Wald’s identity (26.3).

  4. Plot the ratios of simulated \(E_0[N]\) and \(E_1[N]\) to Wald’s approximations (26.4) against \(KL(f_0, f_1)\) and explain the pattern.