[wald_friedman] sprt_single_run can loop forever: NaN log-likelihood ratio from an endpoint beta draw

[mmcky]

sprt_single_run in lectures/wald_friedman.md can enter a loop it cannot leave. The four cells the lecture ships all terminate, but they do so by luck of their seed ranges, not by construction — raising N, changing a seed, or adding another U-shaped f0/f1 pair can hang a notebook build permanently, with no error and no output.

Found while running a benchmark triage evaluation on this lecture; it cost that run two ~10-minute all-core hangs before the cause was identified. The acceleration question is separate and is discussed in QuantEcon/skills — this issue is only about the non-termination, which is worth fixing regardless of that outcome.

The mechanism

The kernel at cell @507 has no iteration cap:

@njit
def sprt_single_run(a0, b0, a1, b1, logA, logB, true_f0, seed):
    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
        log_L += np.log(p(z, a1, b1)) - np.log(p(z, a0, b0))
        if log_L >= logA:
            return n, False
        elif log_L <= logB:
            return n, True

The density at cell @298 is p(x, a, b) = r * x**(a-1) * (1-x)**(b-1), which diverges at the endpoints whenever a < 1 or b < 1. np.random.beta(0.5, 0.4) does return exactly 1.0 — a U-shaped Beta puts its mass at the endpoints and the float64 sampler rounds there. At that draw both p(z, a1, b1) and p(z, a0, b0) are inf, so the increment is log(inf) - log(inf), which is nan, and log_L is nan from that step onward.

Every comparison against nan is False. log_L >= logA is False, log_L <= logB is False, so neither return branch can ever fire and the while True spins forever. The process does not crash, does not warn, and does not stop.

Evidence

Reachable with the lecture's own params_3 = Beta(0.5, 0.4) vs Beta(0.4, 0.5), the third pair at cell @695. Scanning that parameter set with an iteration-capped copy of the kernel:

Quantity	Value
Non-terminating paths in seeds 1..200,000	1 (first at seed 12837)
Rate	~1 path in 200,000
Step at which log_L first goes non-finite, seed 12837	n = 58, with z = 1.0
p(1.0, 0.4, 0.5) and p(1.0, 0.5, 0.4)	both inf

The four shipped cells were each replayed at their exact (seed0, N) and all four terminate:

Cell	Parameters	seed0	N	Non-terminating
@507 headline	Beta(2,5) vs Beta(5,2)	1	20,000	0
@695 params_1	Beta(2,8) vs Beta(8,2)	42	5,000	0
@695 params_2	Beta(4,5) vs Beta(5,4)	42	5,000	0
@695 params_3	Beta(0.5,0.4) vs Beta(0.4,0.5)	42	5,000	0

So the published lecture builds. But params_3 at seed=1, N=20000 — a change no more exotic than reusing the headline cell's seeding — walks straight into seed 12837 and hangs. That is how the triage run hit it, twice.

There is also a milder version of the same hazard for a > 1, b > 1 parameters: there p goes to zero at the endpoints, so a draw of exactly 0.0 or 1.0 gives log(0) - log(0), which is -inf - (-inf), also nan. It is far less likely because the density vanishes there rather than diverging, but the code path is identical.

Suggested fix

The lecture's own sibling kernels already guard, which makes the headline one the odd case out:

• markov_sprt_single_run at cell @1104 caps at max_n = 10000
• var_sprt_single_run at cell @1260 caps at max_T = 500
• sprt_single_run at cell @507 is unbounded

Adding a max_n bound to sprt_single_run would match the lecture's existing convention and bound the damage. Guarding the increment itself would be better, since a cap alone silently returns a fabricated stopping time when the real problem is a nan. Something along the lines of skipping or resampling a draw whose log-density is not finite, or clipping z away from the open interval's endpoints before evaluating p.

A related option, which would remove the endpoint blow-up rather than catch it, is to compute the log-density directly with gammaln instead of forming r * x**(a-1) * (1-x)**(b-1) and taking its log. That is also better numerically for the small shape parameters this lecture uses. It is a larger edit and changes what the reader sees, so it is worth deciding separately from the immediate guard.

Happy to open a PR for whichever shape is preferred — a minimal cap, or the guard plus the gammaln reformulation.

Reproducing

Requires numba, since the RNG stream is numba's:

import numpy as np
from numba import njit

# p and the kernel as defined in the lecture, plus an iteration cap
@njit
def capped(a0, b0, a1, b1, logA, logB, true_f0, seed, cap=200000):
    log_L, n = 0.0, 0
    np.random.seed(seed)
    while n < cap:
        z = np.random.beta(a0, b0) if true_f0 else np.random.beta(a1, b1)
        n += 1
        log_L += np.log(p(z, a1, b1)) - np.log(p(z, a0, b0))
        if log_L >= logA:
            return n, 0
        elif log_L <= logB:
            return n, 1
    return n, -1          # -1 => would have hung

A, B = (1 - 0.10) / 0.05, 0.10 / (1 - 0.05)
print(capped(0.5, 0.4, 0.4, 0.5, np.log(A), np.log(B), True, 12837))
# -> (200000, -1)   with log_L nan from step 58

Drop the cap and the same call does not return.

[mmcky]

Context on how this was found: it surfaced during a benchmark-plugin triage evaluation of this lecture. The plugin-side findings from that same run, including the missing subprocess timeout that turned this hang into two ~10-minute stalls, are in QuantEcon/skills#14.