Suppose a sequential probability ratio test (SPRT) with $$\alpha=P(\text{Type I error})=P(\text{Type II error})=\beta = 10\%$$ is used to decide between $$\ H_0: p=1/2\quad\text{vs.}\quad H_1: p=b\,$$ when in fact the data are$$(X_1,X_2,...)\sim\text{i.i.d. Bernoulli($p=1/2$)}.$$ Let $EN(b)$ denote the expected number of data the SPRT requires to reach a decision (i.e. the expected stopping time).

Simulations strongly suggest that $EN(b)$ has a countable infinity of "jump discontinuities" (most being extremely small) occurring at those $b\in(0,1)$ such that $2^{s+t} b^s (1-b)^t\in\{1/9,9\}$ for some $s,t\in\mathbb{Z^+}.$ (This comes from setting the likelihood ratio $$2^n b^{\sum_{i=1}^nX_i}(1-b)^{n-\sum_{i=1}^nX_i}$$ equal to the "threshold" values used by the SPRT, namely ${\beta\over 1-\alpha}=1/9,{1-\beta\over \alpha}=9$.) Some of these discontinuities seem highly counter-intuitive due to the direction of the jump.

Question 1: How can it happen that $EN(b)$ has a negative jump for some $b<1/2$? Or a positive jump for some $b>1/2$? (Examples include $b=b^*=0.05054567564...$ and $b=1-b^*$, which are among the roots of $2^8b^1(1-b)^7=9$ and $2^8b^7(1-b)^1=9$, respectively.) Is there an intuitive explanation?

NB: Intuitively, making $b$ closer to $1/2$ should make it harder for the SPRT to decide between the two, so this should cause $EN(b)$ to increase. Thus, for $\epsilon>0$ we should have $EN(b) < EN(b+\epsilon)$ when $b<1/2$, and $EN(b) > EN(b+\epsilon)$ when $b>1/2$. But at $b=b^*<1/2$ and $b=1-b^*>1/2$ just the opposite happens!

Question 2:The locations of the jumps are given by roots of the stated polynomial equations, but is there some way to determine the sign of the jump?

The following three pictures, from left-to-right, zoom in on the counter-intuitive negative jump at $b^*=0.05054567564...$ (the red arrows point to $b^*$ and the very light grey dotted lines are $3\sigma$ error bounds on the estimated $EN$):

overview of EN(b), zooming in on a counter-intuitive negative jump

The following Python code clearly identifies the negative jump at $b^*$ ($10^8$ trials may take few minutes, and $10^7$ may be sufficient):

import math
import random

def trial(b):
    factor1, factor0 = 2*b, 2*(1-b)
    Rlo, R, Rhi = 1/9, 1, 9 
    N = 0
    while True: 
        N += 1
        X = 1 if (random.random() <= .5) else 0 
        R *= (factor1 if X else factor0)
        if not (Rlo <= R <= Rhi): 
            return N 
def MonteCarlo(n, b):
    """simulate n trials to estimate the expected outcome and std.dev. of the est. 
    ('Welford's online algorithm')"""
    M = SSD = 0
    for i in range(1,n+1):
        outcome = trial(b)
        DEV = outcome - M
        M += DEV/i
        SSD += DEV*(outcome - M)
    return M, math.sqrt(SSD/(n-1)/n)  # M = sample mean, sqrt(SSD/(n-1)/n) = est.std.dev. of M

# show the negative jump at b = 0.05054567564...
for b in [.0505455, .0505456, .0505457, .0505458]:
    (avgN, sd_avgN) = MonteCarlo(10**8, b) 
    print(f"est. EN(b={b}) +- 3*std.dev. = {avgN:.3f} +- {3*sd_avgN:.3f}")

          > est. EN(b=0.0505455) +- 3*std.dev. = 2.936 +- 0.001
          > est. EN(b=0.0505456) +- 3*std.dev. = 2.936 +- 0.001
          > est. EN(b=0.0505457) +- 3*std.dev. = 2.907 +- 0.001
          > est. EN(b=0.0505458) +- 3*std.dev. = 2.907 +- 0.001

(Some time ago I noticed this phenomenon in a similar scenario and asked about it at MathSE, but got no replies.)


A suggestion:

Consider coupling the SPRT for two values $b_0<b_1$, so you use the same random numbers x but different $b_i$. For the downward step, you must at least sometimes have the test stop for $b_1$ before it stops for $b_0$.

  1. You can actually run the simulation that way. You don't need to record results for the 99.9+% of trials where either the $b_0$ test stops first or they stop at the same time. You can just record trials where the result is surprising. That means you can record it in more detail, which should help with understanding. You will also not need as large a sample size, because the variation in x has been taken out of consideration

  2. You can consider analytically how the two coupled tests might have different results, which may be simpler since x is the same for both of them.

Looking at 2, there seem to be at least two possibilities. First, the counterintuitive jump is explained by an increase in the probability of crossing the 'wrong' boundary. You could check that by running a one-sided version of the test that only stops at the correct boundary and verifying that it didn't have the counterintuitive steps. Second, rounding error: is it possible that multiplying together 7 copies of $b$ and 1 of $1-b$ gives a very slightly different answer depending on the order of multiplication, and this matters right at the edge of the step?

You could also check the rounding error hypothesis by using a higher precision accumulator for the likelihood ratio (eg, cumprod in R uses long double if the compiler provides it; I don't know what's available in Python)

  • $\begingroup$ Thanks for these interesting suggestions -- I'm mulling them over. (Rounding error can be excluded, I think, because the phenomenon occurs also using SageMath with arbitrary-precision reals; e.g., 1000-bit precision). $\endgroup$
    – r.e.s.
    Jul 30 at 16:50
  • $\begingroup$ Ok, so the high-precision test does rule out rounding error, which is good to known. I think coupling will still be useful. $\endgroup$ Jul 31 at 3:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.