Time-uniform, nonparametric, nonasymptotic confidence sequences

Introduction

It has become standard practice for organizations with online presence to run large-scale randomized experiments, or “A/B tests”, to improve product performance and user experience. Such experiments are inherently sequential: visitors arrive in a stream and outcomes are typically observed quickly relative to the duration of the test. Results are often monitored continuously using inferential methods that assume a fixed sample, despite the known problem that such monitoring inflates Type I error substantially (Armitage et al., 1969; Berman et al., 2018). Furthermore, most A/B tests are run with little formal planning and fluid decision-making, compared to clinical trials or industrial quality control, the traditional applications of sequential analysis.

With only a uniform lower bound $(L_{t})$, i.e., if $U_{t}\equiv\infty$, we have a lower confidence sequence. Likewise, if $L_{t}\equiv-\infty$ we have an upper confidence sequence given by $(U_{t})$. Theorems 1, 2, and 3 and Lemma 2 are our key tools for constructing confidence sequences. All build upon the general framework for uniform exponential concentration introduced in Howard et al. (2020), which means our techniques apply in diverse settings: scalar, matrix, and Banach-space-valued observations, with possibly unbounded support; self-normalized bounds applicable to observations satisfying weak moment or symmetry conditions; and continuous-time scalar martingales. Our methods allow for flexible control of the “shape” of the confidence sequence, that is, how the sequence of intervals shrinks in width over time. As a simple example, given a sequence of i.i.d. observations $(X_{t})_{t=1}^{\infty}$ from a 1-sub-Gaussian distribution whose mean $\mu$ we would like to estimate, Theorem 1 yields the following $(1-\alpha)$-confidence sequence for $\mu$, a special case of the more general bound (10):

The $\mathcal{O}(\sqrt{t^{-1}\log\log t})$ asymptotic rate of this bound matches the lower bound implied by the law of the iterated logarithm (LIL), and nonasymptotic bounds of this form are called finite LIL bounds (Jamieson et al., 2014).

We develop confidence sequences that possess the following properties:

Nonasymptotic and nonparametric: our confidence sequences offer coverage guarantees for all sample sizes, without exact distributional assumptions or asymptotic approximations.

Unbounded sample size: our methods do not require a final sample size to be chosen ahead of time. They may be tuned for a planned sample size but always permit additional sampling.

Arbitrary stopping rules: we make no assumptions on the stopping rule used by an experimenter to decide when to end the experiment, or when to act on certain inferences.

Asymptotically zero width: the interval widths of our confidence sequences shrink towards zero at a $1/\sqrt{t}$ rate, ignoring log factors, just as with pointwise confidence intervals.

These properties give us strong guarantees and broad applicability. An experimenter may always choose to gather more samples, and may stop at any time according to any rule—the resulting inferential guarantees hold under the stated assumptions without any approximations. Of course, this flexibility comes with a cost: our intervals are wider than those that rely on asymptotics or make stronger assumptions, for example, a known stopping rule. Typical, fixed-sample confidence intervals derived from the central limit theorem do not satisfy any of (P1)-(P3), and accommodating any one property necessitates wider intervals; we illustrate this in Figure 1. It is perhaps surprising that these four properties come at a numerical cost of less than doubling the fixed-sample, asymptotic interval width—the discrete mixture bound illustrated in Figure 9 stays within a factor of two of the fixed-sample CLT bounds over five orders of magnitude in time.

The idea of a confidence sequence goes back at least to Darling and Robbins (1967a). They are called repeated confidence intervals by Jennison and Turnbull (1984, 1989) (with a focus on finite time horizons) and always-valid confidence intervals by Johari et al. (2015). They are sometimes labeled anytime confidence intervals in the machine learning literature (Jamieson and Jain, 2018).

Recent interest in confidence sequences has come from the literature on best-arm identification with fixed confidence for multi-armed bandit problems. Garivier (2013), Jamieson et al. (2014), Kaufmann et al. (2016), and Zhao et al. (2016) present methods satisfying properties (P1)-(P4) for independent, sub-Gaussian observations. Our results are sharper and more general, and our Bernstein confidence sequence scales with the true variance in nonparametric settings. Confidence sequences are a key ingredient in best-arm selection algorithms (Jamieson and Nowak, 2014) and related methods for sequential testing with multiple comparisons (Yang et al., 2017; Malek et al., 2017; Jamieson and Jain, 2018). Our results improve and generalize such methods.

Maurer and Pontil (2009) and Audibert et al. (2009) prove empirical-Bernstein bounds for fixed times or finite time horizons. Our empirical-Bernstein bound holds uniformly over infinite time. Balsubramani (2014) takes a different approach to deriving confidence sequences satisfying properties (P1)-(P4) by lower bounding a mixture martingale. This work was extended in Balsubramani and Ramdas (2016) to an empirical-Bernstein bound, the only infinite-horizon, empirical-Bernstein confidence sequence we are aware of in prior work. Our result removes a multiplicative pre-factor and yields sharper bounds. We emphasize that our proof technique is quite different from all three of these existing empirical-Bernstein bounds; see Section A.8.

The literature on self-normalized bounds makes extensive use of the method of mixtures, sometimes called pseudo-maximization (de la Peña et al., 2004, 2007; de la Peña, Klass and Lai, 2009; de la Peña, Lai and Shao, 2009; Garivier, 2013); these works introduced the idea of using a mixture to bound a quantity with a random intrinsic time $V_{t}$. These results are mostly given for fixed samples or finite time horizon, though de la Peña et al. (2004, Eq. 4.20) includes an infinite-horizon curve-crossing bound. Lai (1976b) treats confidence sequences for the parameter of an exponential family using mixture techniques similar to those of Section 3.2. Like most work on the method of mixtures, Lai’s work focused on the parametric setting (which we discuss in Section 4.4), while we focus on the application of mixture bounds to nonparametric settings.

Johari et al. (2017) adopt the mixture approach for a commercial A/B testing platform, where properties (P2) and (P3) are critical to provide an “off-the-shelf” solution for a variety of clients. Their application relies on asymptotics which lack rigorous justification. In Section 4.2 we give nonasymptotic justification for a similar confidence sequence under a finite-sample randomization inference model, and in Section 5 we demonstrate how our methods control Type I error in situations where asymptotics fail.

2 Outline

We organize our results using the sub-Gaussian, sub-gamma, sub-Bernoulli, sub-Poisson and sub-exponential settings defined in Section 2.

The stitching method gives new closed-form sub-Gaussian or sub-gamma boundaries (Theorem 1). Our sub-gamma treatment extends prior sub-Gaussian work to cover any martingale whose increments have finite moment-generating function in a neighborhood of zero; see Proposition 1. Our proof is transparent and flexible, accommodating a variety of boundary shapes, including those growing at the rate $\mathcal{O}(\sqrt{t\log\log t})$ with a focus on tight constants, though we do not recommend this bound in practice unless closed-form simplicity is vital.

Conjugate mixtures give one- and two-sided boundaries for the sub-Bernoulli, sub-Gaussian, sub-Poisson and sub-exponential cases (Section 3.2) which avoid approximations made for analytical convenience. The sub-Gaussian boundaries are unimprovable without further assumptions (Section 3.6). These boundaries include a common tuning parameter which is critical in practice and we discuss why their $\mathcal{O}(\sqrt{t\log t})$ growth rate may be preferable to the slower $\mathcal{O}(\sqrt{t\log\log t})$ rate (Section 3.5).

Discrete mixtures facilitate numerical computation of boundaries with a great deal of flexibility, at the cost of slightly more involved computations (Theorem 2). Like conjugate mixture boundaries, these boundaries avoid unnecessary approximations and are unimprovable in the sub-Gaussian case.

Finally, for sub-Gaussian processes, the inverted stitching method (Theorem 3) gives numerical upper bounds on the crossing probability of any increasing, strictly concave boundary over a limited time range. We show that any such boundary yields a uniform upper tail inequality over a finite horizon, and compute its crossing probability.

Building on this foundation, we present a a state-of-the-art empirical-Bernstein bound (Theorem 4) for any sequence of bounded observations using a new self-normalization proof technique. We illustrate our methods with two novel applications: the nonasymptotic, sequential estimation of average treatment effect in the Neyman-Rubin potential outcomes model (Section 4.2), and the derivation of uniform matrix bounds and covariance matrix confidence sequences (Corollaries 3 and 4.3). We give simulation results in Section 5. Section 6 discusses the relationship of our work to existing concepts of sequential testing. Proofs of main results are in Appendix A, with others deferred to Appendix C.

Preliminaries: linear boundaries

Note that an assumption on the upper tail of $(S_{t})$ yields a lower confidence sequence for $(\mu_{t})$; a corresponding assumption on the lower tail of $(S_{t})$ yields an upper confidence sequence for $(\mu_{t})$. In this paper we formally focus on upper tail bounds, from which lower tail bounds can be derived by examining $(-S_{t})$ in place of $(S_{t})$. In general, the left and right tails of $(S_{t})$ may behave differently and require different sets of assumptions, so that our upper and lower confidence sequences may have different forms. Regardless, we can always combine upper and lower confidence sequences using a union bound to obtain a two-sided confidence sequence (1).

When the $(X_{t})$ are independent with common mean $\mu$, the resulting confidence sequence estimates $\mu$, but the setup requires neither independence nor a common mean. In general, the estimand $\mu_{t}$ may be changing at each time $t$; Section 4.2 gives an application to causal inference in which this changing estimand is useful. In principle, $\mu_{t}$ may also be random, although none of our applications involve random $\mu_{t}$.

To construct uniform boundaries $u_{\alpha}$ satisfying inequality (3), we build upon the following general condition (Howard et al., 2020, Definition 1):

When stating that a process is sub-$\psi$, we typically omit $l_{0}$ from our terminology for simplicity. In scalar cases, we always have $l_{0}=1$, while in matrix cases $l_{0}=d$, the dimension of the (square) matrices.

Although $u$ does depend on the constant $l_{0}$ in Definition 1, for simplicity we typically omit this dependence from our notation, writing simply that $u$ is a sub-$\psi$ uniform boundary.

Five particular $\psi$ functions play important roles in our development; below, we take $1/0=\infty$ in the upper bounds on $\lambda$:

$\psi_{B,g,h}(\lambda)\coloneqq\frac{1}{gh}\log\left(\frac{ge^{h\lambda}+he^{-g\lambda}}{g+h}\right)$ on $0\leq\lambda<\infty$, the scaled CGF of a centered random variable (r.v.) supported on two points, $-g$ and $h$, for some $g,h>0$, for example a centered Bernoulli r.v. when $g+h=1$.

$\psi_{N}(\lambda)\coloneqq\lambda^{2}/2$ on $0\leq\lambda<\infty$, the CGF of a standard Gaussian r.v.

One may freely scale $\psi$ by any positive constant and divide $V_{t}$ by the same constant so that Definition 1 remains satisfied; by convention, we scale all $\psi$ functions above so that $\psi^{\prime\prime}(0_{+})=1$. When we speak of a sub-gamma process (or uniform boundary) with scale parameter $c$, we mean a sub-$\psi_{G,c}$ process (or uniform boundary), and likewise for other cases. We often write $\psi_{B}$, $\psi_{P}$, etc., dropping the range and scale parameters from our notation. As we summarize in Figure 2 and detail in Proposition 11, certain general implications hold among sub-$\psi$ boundaries. In particular, any sub-Gaussian boundary can also serve as a sub-Bernoulli boundary; any sub-Poisson boundary serves as a sub-Gaussian or sub-Bernoulli boundary; and, importantly, any sub-gamma or sub-exponential boundary can serve as a sub-$\psi$ boundary in any of the other four cases. Indeed, a sub-gamma or sub-exponential boundary applies to many cases of practical interest, as detailed below.

Suppose $\psi$ is twice-differentiable and $\smash{\psi(0)=\psi^{\prime}(0_{+})=0}$. Suppose, for each $c>0$, $u_{c}(v)$ is a sub-gamma or sub-exponential uniform boundary with crossing probability $\alpha$ for scale $c$. Then $v\mapsto u_{k_{1}}(k_{2}v)$ is a sub-$\psi$ uniform boundary for some constants $k_{1},k_{2}>0$ depending only on $\psi$.

Proposition 1 restates Howard et al. (2020, Proposition 1), which shows that any process $(S_{t})$ which is sub-$\psi$ is also sub-gamma and sub-exponential, if $\psi$ satisfies the conditions of Proposition 1. Note that these conditions are satisfied for any mean-zero random variable if the CGF exists in a neighborhood of zero, so the conditions are quite weak (Jorgensen, 1997, Theorem 2.3).

Suppose $X_{1},X_{2},\dots$ are i.i.d. draws from a $\mathcal{N}(\mu,\sigma^{2})$ distribution and we wish to sequentially estimate $\sigma^{2}$ when $\mu$ is also unknown. Let $S_{t}\coloneqq\sigma^{-2}\sum_{i=1}^{t+1}(X_{i}-\bar{X}_{t+1})^{2}-t$ for $t=1,2,\dots$, where $\bar{X}_{t}\coloneqq t^{-1}\sum_{i=1}^{t}X_{i}$ is the sample mean. This $S_{t}$ is a centered and scaled sample variance, and as in Darling and Robbins (1967a), we use the fact that $S_{t}$ is a cumulative sum of independent, centered Chi-squared random variables each with one degree of freedom (see Appendix H for details). Such a centered Chi-squared distribution has variance two and CGF equal to $2\psi_{E,2}$.

Thus $(S_{t})$ is 1-sub-exponential with variance process $V_{t}=2t$ and scale parameter $c=2$. We may uniformly bound the upper deviations of $S_{t}$ using any sub-exponential uniform boundary, for example the gamma-exponential mixture boundary of Proposition 9. Or, we can use Proposition 11 to deduce that $(S_{t})$ is sub-gamma with scale $c=2$ (and the same variance process) and use the closed-form stitched boundary of Theorem 1.

The above constructions yield lower confidence sequences for the variance. To obtain an upper confidence sequence, we use the fact that $(-S_{t})$ is 1-sub-exponential with scale parameter $\smash{c=-2}$. Now Proposition 11 implies that $(-S_{t})$ is sub-gamma with scale $c=-1$, so the stitched boundary again applies, while Proposition 11 implies that $(-S_{t})$ is also sub-Gaussian, so we may alternatively use the normal mixture boundary of Proposition 6. Since $\psi_{G,-1}$ is uniformly smaller than $\psi_{N}$, the above analysis yields tighter bounds than the sub-Gaussian approach of Darling and Robbins (1967a).

The simplest uniform boundaries are linear with positive intercept and slope. This is formalized in Howard et al. (2020), partially restated below.

For any $\lambda\in[0,\lambda_{\max})$ and $\alpha\in(0,1)$,

is a sub-$\psi$ uniform boundary with crossing probability $\alpha$.

While Lemma 1 provides a versatile building block, the $\mathcal{O}(V_{t})$ growth of $u(V_{t})$ may be undesirable. Indeed, from a concentration point of view, the typical deviations of $S_{t}$ tend to be only $\mathcal{O}(\sqrt{V_{t}})$, so the bound will rapidly become loose for large $t$. From a confidence sequence point of view, recall that the confidence radius for the mean is given by $u(V_{t})/t$. Typically, $V_{t}=\Theta(t)$ a.s. as $t\to\infty$, so the confidence radius will be asymptotically zero width if and only if $u(v)=o(v)$. In other words, we cannot achieve arbitrary estimation precision with arbitrarily large samples unless the uniform boundary is sublinear. We address this problem in Section 3, building upon Lemma 1 to construct curved sub-$\psi$ uniform boundaries.

Curved uniform boundaries

We present our four methods for computing curved uniform boundaries in Sections 3.1, 3.2, 3.3, and 3.4. In Section 3.5, we discuss how to tune boundaries, a necessity for good performance in practice, and we describe the unimprovability of sub-Gaussian mixture bounds in Section 3.6.

Our analytical “stitched” bound is useful in the sub-Gaussian case or, more generally, the sub-gamma case with scale $c$. We require three user-chosen parameters:

a scalar $\eta>1$ determines the geometric spacing of intrinsic time,

a scalar $m>0$ which gives the intrinsic time at which the uniform boundary starts to be nontrivial, and

Recalling the scale parameter $c$ for the $\psi_{G}$ function above and the constant $l_{0}$ in Definition 1, we define the stitching function $\mathcal{S}_{\alpha}$ as

The boundary shape is determined by choosing the function $h$ and setting the nominal crossing probability in the $k$th epoch to equal $\alpha/h(k)$. Then Theorem 1 gives a curved boundary which grows at a rate $\mathcal{O}\left(\sqrt{V_{t}\log h(\log_{\eta}V_{t})}\right)$ as $V_{t}\uparrow\infty$. The more slowly $h(k)$ grows as $k\uparrow\infty$, the more slowly the resulting boundary will grow as $V_{t}\uparrow\infty$. A simple choice is exponential growth, $h(k)=\eta^{sk}/(1-\eta^{-s})$ for some $s>1$, yielding $\mathcal{S}_{\alpha}(v)=\mathcal{O}(\sqrt{v\log v})$. A more interesting example is $h(k)=(k+1)^{s}\zeta(s)$ for some $s>1$, where $\zeta(s)$ is the Riemann zeta function. Then, when $l_{0}=1$, Theorem 1 yields the polynomial stitched boundary: for $c\geq 0$,

where the second term is neglected in the sub-Gaussian case since $c=0$. This is a “finite LIL bound”, so-called because $\mathcal{S}_{\alpha}(v)\sim\sqrt{sk_{1}^{2}v\log\log v}$, matching the form of the law of the iterated logarithm (Stout, 1970). We can bring $sk_{1}^{2}$ arbitrarily close to $2$ by choosing $\eta$ and $s$ sufficiently close to one, at the cost of inflating the additive term $\log(\zeta(s)/(\log^{s}\eta))$. Briefly, increasing $\eta$ increases the size of each epoch in the aforementioned peeling argument, which reduces the looseness of the union bound over epochs. But the larger we make the epochs, the further each linear boundary deviates from the ideal curved shape at the ends of the epochs, which inflates our final boundary. The choice of $s$ involves a similar tradeoff: increasing $s$ causes us to exhaust more of our total error probability budget on earlier epochs, decreasing the constant term (which matters most for early times), at the cost of a union bound over smaller error probabilities in later epochs, which shows up as an increase in the leading constant. We discuss parameter tuning in more practical terms in Section 3.5. For example, take $\smash{\eta=2,s=1.4,m=1}$; if $S_{t}$ is a sum of independent, zero-mean, 1-sub-Gaussian observations, we obtain

Figure 9 in Appendix G compares a sub-Gaussian stitched boundary to a numerically-computed discrete mixture bound with a mixture distribution roughly corresponding to $h(k)\propto(k+1)^{1.4}$, as described in Section A.6. This discrete mixture boundary acts as a lower bound (see Section 3.6) and shows that not too much is lost by the approximations involved in the stitching construction. Figure 10 compare the same stitched boundary to related bounds from the literature; our bound shows slightly improved constants over the best known bounds.

Although our stitching construction begins with a sub-gamma assumption, it applies to other sub-$\psi$ cases, including sub-Bernoulli, sub-Poisson and sub-exponential cases; see Figure 2 and Proposition 1. Further, our stitched bounds apply equally well in continuous-time settings to Brownian motion, continuous martingales, martingales with bounded jumps, and martingales whose jumps satisfy a Bernstein condition; see Corollary 8.

While our focus is on nonasymptotic results, Theorem 1 makes it easy to obtain the following general upper asymptotic LIL, proved in Section A.2:

Suppose $(S_{t})$ is sub-$\psi$ with variance process $(V_{t})$ and $\psi(\lambda)\sim\lambda^{2}/2$ as $\lambda\downarrow 0$. Then

2 Conjugate mixture boundaries

For any probability distribution $F$ on $[0,\lambda_{\max})$ and $\alpha\in(0,1)$,

is a sub-$\psi$ uniform boundary with crossing probability $\alpha$, so long as the supermartingale $(L_{t})$ of Definition 1 is product measurable when the underlying probability space is augmented with the independent random variable $\lambda$.

one-, two-sided normal mixture boundaries (sub-Gaussian case);

one-, two-sided beta-binomial mixture boundaries (sub-Bernoulli case);

one-sided gamma-Poisson mixture boundary (sub-Poisson case); and

one-sided gamma-exponential mixture boundary (sub-exponential case).

The two-sided normal mixture boundary has a closed form expression:

The one-sided normal mixture boundary has a similar, closed-form upper bound, making these especially convenient. It is clear from (14) that the normal mixture boundary grows as $\mathcal{O}(\sqrt{v\log v})$ asymptotically, and this rate is shared by all of our conjugate mixture boundaries. Indeed, Proposition 2 below, proved in Section A.4, shows that such a rate holds for any mixture boundary as given by (13) whenever the mixing distribution is continuous with positive density at and around the origin, a property which holds for all mixture distributions used in our conjugate mixture boundaries, subject to regularity conditions on $\psi$ which hold for the CGF of any nontrivial, mean-zero r.v. and specifically for the five $\psi$ functions in Section 2.

Assume $(i)$ $\psi$ is nondecreasing, $\psi(0)=\psi^{\prime}(0_{+})=0$, $\psi^{\prime\prime}(0_{+})=c>0$, and $\psi$ has three continuous derivatives on a neighborhood including the origin; and $(ii)$ $F$ has density $f$ (w.r.t. Lebesgue) which is continuous and positive on a neighborhood including the origin. Then

Note that $f$ need not place mass on all of $[0,\lambda_{\max})$, only near the origin, for the asymptotic rate to hold. Proposition 2 shows how the asymptotic behavior of any such mixture bound depends only on the behavior of $\psi$ and $f$ near the origin, a result reminiscent of the central limit theorem. Analogous, related results for the sub-Gaussian special case using $\psi(\lambda)=\lambda^{2}/2$ can be found in Robbins and Siegmund (1970, Section 4) and Lai (1976a, Theorem 2), in some cases under weaker assumptions on $F$.

In contrast to previous derivations of conjugate mixture boundaries in the literature, all of our conjugate mixture boundaries include a common tuning parameter $\rho>0$ which controls the sample size for which the boundary is optimized. Such tuning is critical in practice, as we explain in Section 3.5, but has been ignored in much prior work. Additionally, with the exception of the sub-Gaussian case, most prior work on the method of mixtures has focused on parametric settings. We instead emphasize the applicability of these bounds to nonparametric settings. For example, when the observations are bounded, one may construct a confidence sequence making use of empirical-Bernstein estimates (Theorem 4) based on our gamma-exponential mixture (Proposition 9). See Appendix J for other conditions in which mixture bounds yield nonparametric uniform boundaries.

3 Numerical bounds using discrete mixtures

In applications, one may not need an explicit closed-form expression so long as the bound can be easily computed numerically. Our discrete mixture method is an efficient technique for numerical computation of curved boundaries for processes satisfying Definition 1. It permits arbitrary mixture densities, thus producing boundaries growing at the rate $\mathcal{O}(\sqrt{v\log\log v})$. Recall that the shape of the stitched bound was determined by the user-specified function $h$. For the discrete mixture bound, one instead specifies a probability density $f$ over finite support $(0,\overline{\lambda}]$ for some $\overline{\lambda}\in(0,\lambda_{\max})$. We first discretize $f$ using a series of support points $\lambda_{k}$, geometrically spaced according to successive powers of some $\eta>1$, and an associated set of weights $w_{k}$:

is a sub-$\psi$ uniform boundary with crossing probability $\alpha$.

We suppress the dependence of $\textup{DM}_{\alpha}$ on $f$, $l_{0}$, $\overline{\lambda}$ and $\eta$ for notational simplicity. Though Theorem 2 is a straightforward consequence of the method of mixtures, our choice of discretization (16) makes it effective, broadly applicable, and easy to implement. See Section A.5 for the proof of this result. Figure 9 includes an example bound, demonstrating a slight advantage over stitching. Section A.6 describes a connection between the stitching and discrete mixture methods, including a correspondence between the alpha-spending function $h$ and the mixture density $f$. Finally, we note that the method can be applied even when $f$ is not monotone; one must simply choose the discretization (16) more carefully, using known properties of $f$.

4 Inverted stitching for arbitrary boundaries

In the method of mixtures, we choose a mixing distribution $F$ and the machinery yields a boundary $\mathcal{M}_{\alpha}$. Likewise, in the stitching construction of Theorem 1, we choose an error decay function $h$ and obtain a boundary $\mathcal{S}_{\alpha}$. Here, we invert the procedure: we choose a boundary function $g(v)$ and numerically compute an upper bound on its $S_{t}$-upcrossing probability using a stitching-like construction.

is a sub-Gaussian uniform boundary with crossing probability at most

The proof is in Section A.7. For simplicity we restrict to the sub-Gaussian case; examination of the proof will show that the method applies in other sub-$\psi$ cases as well, since we simply apply Lemma 1 to appropriately chosen lines, but more involved numerical calculations will be necessary, as the closed-form (19) no longer applies. A similar idea was considered by Darling and Robbins (1968), using a mixture integral approximation instead of an epoch-based construction to derive closed-form bounds. Theorem 3 requires numerical summation but yields tighter bounds with fewer assumptions. As an example, Theorem 3 with $\eta=2.99$ shows that

This boundary is illustrated in Figure 9.

5 Tuning boundaries in practice

All uniform boundaries involve a tradeoff of tightness at different intrinsic times: making a bound tighter for some range of times requires making it looser at other times. Roughly speaking, the choice of a uniform boundary involves choosing both what time the bound should be optimized for (e.g., should the bound be tightest around 100 observations or around 100,000 observations?) as well as how quickly the bound degrades as we move away from the optimized-for time (e.g., if we optimize for 100 samples, will the bound be twice as wide when we reach 1,000 samples, or will it stay within a factor of two until we reach 1,000,000 samples?). A boundary which decays more slowly will necessarily not be as tight around the optimized-for time. In brief, linear boundaries decay the most quickly, conjugate mixture boundaries decay substantially more slowly, and polynomial stitched boundaries decay even more slowly; we feel that mixture boundaries strike a good balance in practice.

Here, we explain how to optimize uniform boundaries for a particular time and discuss the above tradeoff in more detail. Let $W_{-1}(x)$ be the lower branch of the Lambert $W$ function, the most negative real-valued solution in $z$ to $ze^{z}=x$. Consider the unitless process $S_{t}/\sqrt{V_{t}}$, and the corresponding uniform boundary $v\mapsto u(v)/\sqrt{v}$. Since all of our uniform boundaries $u(v)$ have positive intercept at $v=0$, and all grow at least at the rate $\sqrt{v\log\log v}$ as $v\to\infty$, the normalized boundary $u(v)/\sqrt{v}$ diverges as $v\to 0$ and $v\to\infty$. For the two-sided normal mixture (14), there is a unique time $m$ at which $u(v)/\sqrt{v}$ is minimized; $m$ is proportional to tuning parameter $\rho$ as follows:

Let $u(v)$ be the two-sided normal mixture boundary (14) with parameter $\rho>0$.

For fixed $\rho>0$, the function $v\mapsto u(v)/\sqrt{v}$ is uniquely minimized at $v=m$ with $m$ given by

For fixed $m>0$, the choice of $\rho$ which minimizes the boundary value $u(m)$ is also determined by (21).

The above result is proved in Section C.1; it is a matter of elementary calculus, but addresses a question that has received little attention in the literature. Figure 4 includes the normalized versions of two normal mixture boundaries optimized for different times, $m=300$ and $m=$ 5,000. Optimizing for the range of values of $V_{t}$ most relevant in a particular application will yield the tightest confidence sequences. However, as the figure shows, one need not have a very precise range of times, so long as one uses a conservatively low value for $m$, because $u(v)/\sqrt{v}$ grows slowly after time $m$. Indeed, for the normal mixture boundary with $\alpha=0.05$ and $l_{0}=1$, we have $u(m)/\sqrt{m}\approx 3.0$ and $u(100m)/\sqrt{100m}\approx 3.6$, so that the penalty for being off by two orders of magnitude is modest.

The one-sided normal mixture boundary of Appendix Proposition 6 with crossing probability $\alpha$ is nearly identical to the two-sided normal mixture boundary with crossing probability $2\alpha$, so one may choose $\rho$ as in Proposition 3 with $\alpha$ doubled. For the gamma-exponential mixture and other non-sub-Gaussian uniform boundaries, Proposition 3 provides a good approximation in practice. Figure 4 includes gamma-exponential mixture boundaries with the same $\rho$ values as each corresponding normal mixture boundary. Though the normalized gamma-exponential mixture boundary with $m=300$ clearly reaches its minimum at $v>m$, this choice of $\rho$ seems reasonable. Discrete mixtures can be similarly tuned by adjusting the precision of the mixing distribution, but require additional considerations (Appendix E).

Comparing the sub-Gaussian stitched boundary, discrete mixture boundary, and normal mixture boundary optimized for $m=300$ in Figure 4 illustrates another important point for practice: although the normal mixture bound grows more quickly than the others as $v\to\infty$, it remains smaller over about three orders of magnitude. This makes it preferable for many real-world applications, as the longest feasible duration of an experiment is rarely more than two orders of magnitude larger than the earliest possible stopping time. For example, many online experiments run for at least one week to account for weekly seasonality effects, and very few such experiments last longer than 100 weeks. As both the normal mixture and the discrete mixture are unimprovable in general (Section 3.6), the difference is attributable to the choice of mixture, or alternatively, to the fact that the normal mixture trades tightness around the optimized-for time in exchange for looseness at much later times. The lesson is that the $\mathcal{O}(v\log\log v)$ rate, while asymptotically optimal in certain settings and useful for theory and some applications, may not be preferable in all real-world scenarios.

6 Unimprovability of uniform boundaries

Definition 2 of a sub-$\psi$ boundary $u$ involves only an upper bound on the $u$-crossing probability of any sub-$\psi$ process $(S_{t})$. One may reasonably ask for corresponding lower bounds on the $u$-crossing probability to quantify how tight this boundary is. In the ideal case, we might desire a boundary $u$ such that the true $u$-crossing probability of some process $(S_{t})$ is equal to the upper bound. In nonparametric settings, we cannot achieve this goal for every sub-$\psi$ process. However, we might still ask that there exists some sub-$\psi$ process for which the true $u$-crossing probability is arbitrarily close to the upper bound, so that the upper bound on crossing probability is unimprovable in general. That is, we might ask that the inequality on the supremum in Definition 2 holds with equality.

For any exact, 1-sub-Gaussian mixture boundary $\mathcal{M}_{\alpha}$,

We prove Proposition 4 in Section C.2. In general, for each $\alpha$ there is an infinite variety of boundaries that are unimprovable in the above sense, differing in when they are loose and tight. These different boundaries will yield confidence sequences which are loose or tight at different sample sizes, or, equivalently, are efficient for detecting different effect sizes. Such a boundary cannot be tightened everywhere without increasing the crossing probability.

Applications

After presenting an empirical-Bernstein confidence sequence for bounded observations, we apply our uniform boundaries to causal effect estimation and matrix martingales. We also consider estimation for a general, one-parameter exponential family.

Suppose $X_{t}\in[a,b]$ a.s. for all $t$. Let $(\widehat{X}_{t})$ be any $[a,b]$-valued predictable sequence, and let $u$ be any sub-exponential uniform boundary with crossing probability $\alpha$ for scale $c=b-a$. Then

This is an empirical-Bernstein bound because it uses the sum of observed squared deviations to estimate the true variance, much like a classical $t$-test. Hence the confidence radius scales with the true standard deviation for sufficiently large samples, regardless of the support diameter $b-a$, and with no prior knowledge of the true variance. Note also that this bound does not require that observations share a common mean.

For an explicit example, assume $X_{i}\in$ and define the empirical variance as $\widehat{V}_{t}\coloneqq\sum_{i=1}^{t}(X_{i}-\bar{X}_{i-1})^{2}$. Invoking Theorem 4 with the polynomial stitched bound (10) using $c=1$, $\eta=2$, $m=1$, and $h(k)\propto k^{1.4}$, we have the following 95%-confidence sequence for $\mu_{t}$:

When a closed form is not required, the gamma-exponential mixture (Proposition 9) may yield tighter bounds than stitching; simulations in Section 5 demonstrate the use of Theorem 4 with this mixture.

2 Estimating ATE in the Neyman-Rubin model

At each step $t$, having treated and observed units $1,\dots,t$, we wish to draw inference about the estimand $\textup{ATE}_{t}\coloneqq t^{-1}\sum_{i=1}^{t}[Y_{i}(1)-Y_{i}(0)]$. In particular, we seek a confidence sequence for $(\textup{ATE}_{t})_{t=1}^{\infty}$. To construct our estimator, we may utilize any predictions $\widehat{Y}_{t}(0)$ and $\widehat{Y}_{t}(1)$ for each unit’s potential outcomes; these random variables must be $\mathcal{F}_{t-1}$-measurable, for each $t$. We then employ the inverse probability weighting estimator

which is (conditionally) unbiased for the individual treatment effect $Y_{t}(1)-Y_{t}(0)$. As with Theorem 4, better predictions will lead to shorter confidence intervals, but the coverage guarantee holds for any choice of predictions, and a reasonable choice would be the average of past observed outcomes. See Aronow and Middleton (2013) for a similar strategy for fixed-sample estimation.

We assume bounded potential outcomes; for simplicity we assume $Y_{t}(k)\in$ for all $t\geq 1,k=0,1$, and we assume predictions are likewise bounded. We further assume that treatment probabilities are uniformly bounded away from zero and one. Then, an empirical-Bernstein confidence sequence for $\textup{ATE}_{t}$ follows from Theorem 4, where we use $\widehat{X}_{t}=\widehat{Y}_{t}(1)-\widehat{Y}_{t}(0)$ so that

Suppose $P_{t}\in[p_{\min},1-p_{\min}]$ a.s., $Y_{t}(k)\in$ and $\widehat{Y}_{t}(k)\in$ for all $t\geq 1,k=0,1$. Let $u$ be any sub-exponential uniform boundary with scale $2/p_{\min}$ and crossing probability $\alpha$. Then

3 Matrix iterated logarithm bounds

The result follows using the polynomial stitched boundary after invoking Fact 1(c) and Lemma 2 of Howard et al. (2020) (cf. (Tropp, 2011)), which show that $(S_{t})$ is sub-gamma with variance process $(V_{t})$, scale $c=b/3$, and $l_{0}=d$. Beyond bounded increments, the same bound holds for any sub-gamma process. As evidenced by Proposition 1, this is a very general condition.

Taking $\eta$ and $s$ arbitrarily close to one and using the final result of Theorem 1, we obtain the following asymptotic matrix upper LIL, proved in Section A.9. Here we denote the martingale increments by $\Delta Y_{t}\coloneqq Y_{t}-Y_{t-1}$.

We use a sub-Poisson uniform boundary to obtain a uniform analogue:

For example, using the polynomial stitched bound with scale $c=2b/3$ and $m=b\lVert\Sigma\rVert_{\text{op}}$, Corollary 5 gives a $(1-\alpha)$-confidence sequence for $\Sigma$ with operator norm radius $\mathcal{O}(\sqrt{t^{-1}\log\log t})$. This bound has the closed form

In other words, with high probability, we have for all $t\geq 1$ that

Compared to the fixed-sample result (30), we obtain uniform control by adding a factor of $\log\log t$. We are not aware of other results like these for sequential covariance matrix estimation. Figure 6 illustrates the confidence sequence of Corollary 5 on simulated data using a discrete mixture boundary with the mixture density $f^{\text{LIL}}_{s}$ defined in (85).

4 One-parameter exponential families

Simulations

InThe repository https://github.com/gostevehoward/cspaper contains code to reproduce all simulations and plots in this paper. Uniform boundaries themselves are implemented in R and Python packages at https://github.com/gostevehoward/confseq. Figure 7 we illustrate the error control of some of our confidence sequences for estimating the mean of an i.i.d. sequence of observations $(X_{i})$ with bounded support $[a,b]$. We compare four strategies:

The Hoeffding strategy exploits the fact that bounded observations are sub-Gaussian (Hoeffding, 1963; cf. Howard et al., 2020, Lemma 3(c)). We use a two-sided normal mixture boundary (14) with variance process $V_{t}=(b-a)^{2}t/4$.

The beta-binomial strategy uses the stronger condition that bounded observations are sub-Bernoulli (Hoeffding, 1963; cf. Howard et al., 2020, Fact 1(b)), accounting for the true mean as well as the boundedness, but possibly failing to take account of the true variance. For hypothesized true mean $\mu$, this strategy uses the beta-binomial mixture boundary given in Proposition 7, with parameters $g(\mu)=\mu-a$ and $h(\mu)=b-\mu$, and variance process $V_{t}(\mu)=g(\mu)h(\mu)t$. The confidence set for the mean is $\{\mu\in[a,b]:-f_{g(\mu),h(\mu)}(V_{t}(\mu))\leq\sum_{i=1}^{t}X_{i}-t\mu\leq f_{h(\mu),g(\mu)}(V_{t}(mu))\}$. This is more efficiently computed using the mixture supermartingale $m(S_{t},V_{t})$ of (57), as $\{\mu\in[a,b]:m(\sum_{i=1}^{t}X_{i}-t\mu,V_{t}(\mu))<1/\alpha\}$.

The pointwise Bernoulli strategy uses the same sub-Bernoulli condition as the beta-binomial strategy, but relies on a fixed-sample Cramér-Chernoff bound which is valid pointwise but not uniformly over time. Specifically, we reject mean $\mu$ if $V_{t}\psi_{B}^{\star}(S_{t}/V_{t})\geq\log\alpha^{-1}$, where $S_{t}$ is the sum of centered observations as usual, $V_{t}=(\mu-a)(b-\mu)t$, and we set $\smash{g=\mu-a,h=b-\mu}$ in $\psi_{B}$, with $\psi_{B}^{\star}$ its Legendre-Fenchel transform.

The empirical-Bernstein strategy uses an empirical estimate of variance, thus achieving a confidence width scaling with the true variance in all three cases. Here we use Theorem 4 with a gamma-exponential mixture boundary (Proposition 9). For predictions, we use the mean of past observations: $\widehat{X}_{t}=(t-1)^{-1}\sum_{i=1}^{t-1}X_{i}$.

The naive self-normalized (“Naive SN”) strategy plugs the empirical variance estimate, the sum of squared prediction errors from Theorem 4, into the two-sided normal mixture (14). It ignores the facts that the observations are not sub-Gaussian with respect to their true variance and that the variance is estimated. This strategy is similar to that of Johari et al. (2017) and does not guarantee coverage. Though it will sometimes control false positives, coverage rates can easily be inflated for asymmetric, heavy-tailed distributions, as we illustrate.

We present three cases of bounded distributions. The first case is the easiest, with $\operatorname{Ber}(0.5)$ observations. Here the sub-Gaussian variance parameter based on the boundedness of the observations is equal to the true variance, so the Hoeffding strategy performs well. The empirical-Bernstein strategy is only a little wider, and all four successfully control false positives. The story changes with the more difficult $\operatorname{Ber}(0.01)$ distribution, however. The Hoeffding boundary is far too wide, since it fails to make use of information about the true variance. The beta-binomial bound uses information about variance provided by the first moment to achieve the correct scaling. The naive self-normalized strategy, on the other hand, yields confidence intervals that are too small and fail to control false positive rate. The empirical-Bernstein strategy, though only slightly wider than the naive bound for large sample sizes, gives just enough extra width to control the false positive rate and is nearly as narrow as the beta-binomial bound. The final, three-point distribution takes values $-1.408$ and $1$ with probability $0.495$ each, and takes value $20$ with probability $0.01$. Here the beta-binomial strategy yields confidence intervals that are too wide. In this most difficult case, only the empirical-Bernstein strategy yields tight intervals while controlling false positive rates.

Implications for sequential hypothesis testing

We have organized our presentation around confidence sequences and closely related uniform concentration bounds due to our belief that they offer a useful “user interface” for sequential inference. However, our methods also yield always-valid $p$-values (Johari et al., 2015) for sequential tests. Indeed, a slew of related definitions from the literature are equivalent or “dual” to one another. Here we briefly discuss these connections. The following result, proved in Section C.4, gives equivalent formulations of common definitions in sequential testing.

Let $(A_{t})_{t=1}^{\infty}$ be an adapted sequence of events in some filtered probability space and let $A_{\infty}\coloneqq\limsup_{t\to\infty}A_{t}$. The following are equivalent:

Pros and cons of the running intersection

However, the intersected intervals $\widetilde{\textup{CI}}_{t}$ may become empty at some point. This is particularly likely if the underlying parameter is drifting over time, contrary to the assumption of stationarity or identically-distributed observations, and such a drift would be the likely interpretation of this event in practice. In this non-stationary case, the non-intersected sequence is the more sensible one to use. The solution of Johari et al. (2017) is to “reset” the experiment, discarding data accumulated up to that point, on the rationale that such an event indicates that previous data are no longer relevant to estimation of the current parameter of interest. However, this means that our confidence sequence can go from a very high precision estimate at some time $t$ to knowing almost nothing at time $t+1$, which is difficult for an experimenter to interpret and could lead to misleading inference just before the reset. Jennison and Turnbull (1989) make a case for the non-intersected intervals on slightly different grounds, arguing that estimation at time $t$ ought to be a function of the sufficient statistic at that time. Shifting to the potential outcomes model in Section 4.2 neatly avoids this issue: because the estimand is changing at each time, the non-intersected intervals are the only reasonable choice for estimating $\textup{ATE}_{t}$ and no conceptual difficulty remains.

Summary and future work

We have discussed four techniques for deriving curved uniform boundaries, each improving upon past work, with careful attention paid to constants and to practical issues. By building upon the general framework of Howard et al. (2020), we have emphasized the nonparametric applicability of our boundaries. A leading example of the utility of this approach is the general empirical-Bernstein bound, with an application to sequential causal inference, and we have also shown how our framework immediately yields novel results for matrix martingales.

We introduced the method of mixtures and the epoch-based analyses in Section 1.1. Two other methods of extending the SPRT deserve mention, though they are distinct from our approaches. First, the approach of Robbins and Siegmund (1972, 1974) examines $\prod_{i}f_{\hat{\lambda}_{i-1}}(X_{i})/f_{0}(X_{i})$ where $\hat{\lambda}_{i-1}$ is a “nonanticipating” estimate based on $X_{1},\dots,X_{i-1}$. This is similar to a generalized likelihood ratio but modified to retain the martingale property (cf. Wald (Wald, 1947, section 10.5), (Lorden and Pollak, 2005)). Second, the sequential generalized likelihood ratio approach examines $\sup_{\lambda}\prod_{i}f_{\lambda}(X_{i})/f_{0}(X_{i})$, which is not a martingale under the null (Siegmund and Gregory, 1980; Lai, 1997; Kulldorff et al., 2011).

The concept of test (super)martingales expounded by Shafer et al. (2011) is related to our methods for conducting inference based on Ville’s inequality applied to nonnegative supermartingales. Their main example is the Beta mixture for i.i.d. Bernoulli observations, an example which originated with Ville (1939) and discussed by Robbins (1970) and Lai (1976b). A recent “safe testing” framework of Grünwald et al. (2019) is also tightly related. In terms of these frameworks, our work can be viewed as constructing “safe confidence intervals” (and thus safe tests) using nonparametric test supermartingales.

A very different approach is that of group sequential methods (Pocock, 1977; O’Brien and Fleming, 1979; Lan and DeMets, 1983; Jennison and Turnbull, 2000). These methods rely on either exact discrete distributions or asymptotics to assume exact normality of group increments, either of which permits computation of sequential boundaries via numerical integration. The resulting confidence sequences are tighter than ours, but lack nonasymptotic guarantees or closed-form results and do not support continuous monitoring.

A related problem is that of terminal confidence intervals, in which one assumes a rigid stopping rule and wishes to construct a confidence interval upon termination. Siegmund (1978) gave an analytical treatment of the problem; numerical methods are also available for group sequential tests (Jennison and Turnbull, 2000, section 8.5). However, the idea of a rigid stopping rule is often restrictive.

2 Future work

We discuss in Appendix I how our work may be extended to martingales in smooth Banach spaces and real-valued, continuous-time martingales. It may be fruitful to explore applications in those areas.

Our consideration of optimality has been limited to the discussion in Section 3.6. It would be valuable to further explore various optimality properties for nonasymptotic uniform bounds. For example, it is standard in sequential testing to compute the expected sample size to reject a null under parametric alternatives. Though we target less restrictive assumptions, it may be instructive to compute bounds in special cases. Second, a natural counterpoint to our uniform concentration bounds would be a set of uniform anticoncentration bounds. This would yield a nonasymptotic extension of the “lim inf” half of the classical LIL. Balsubramani (2014, Theorem 3) gives one such interesting result. Last, in practice, one will rarely require updated inference after every observation, and may be content to take observations in groups. Further, one may be satisfied with a finite time horizon Garivier and Leonardi (2011). This is the domain in which group-sequential methods shine, but SPRT-based methods can be made competitive by estimating the “overshoot” of the stopped supermartingale (Lai and Siegmund, 1977, 1979; Siegmund, 1985; Whitehead and Stratton, 1983). It would be interesting to understand whether such improvements work out in nonparametric settings.

Acknowledgments

Howard thanks ONR Grant N00014-15-1-2367. Sekhon thanks ONR grants N00014-17-1-2176 and N00014-15-1-2367. Ramdas thanks NSF grant DMS1916320. We thank Boyan Duan and Ian Waudby-Smith as well as the referees/AE for useful suggestions.

References

Appendix A Proofs of main results

In this section we give proofs of our main results along with selected discussion of and intuition for proof techniques.

The idea behind Theorem 1 is to divide intrinsic time into geometrically spaced epochs, $\eta^{k}\leq V_{t}<\eta^{k+1}$ for some $\eta>1$. We construct a linear boundary within each epoch using Lemma 1 and take a union bound over crossing events of the different boundaries. The resulting, piecewise-linear boundary may then be upper bounded by a smooth, concave function. Figure 3 illustrates the construction.

As discussed in Section 3.1, the function $h$ determines the nominal crossing probability $\alpha/h(k)$ allocated to the $k$th epoch, and we have mentioned the choices $h(k)=\eta^{sk}/(1-\eta^{-s})$ and $h(k)=(k+1)^{s}\zeta(s)$. One may substitute a series converging yet more slowly; for example, $h(k)\propto(k+2)\log^{s}(k+2)$ for $s>1$ yields

matching related analysis in Darling and Robbins (1967b), Robbins and Siegmund (1969), Robbins (1970), and Balsubramani (2014). In practice, the bound (35) appears to behave like bound (10) with worse constants. However, the fact that the stitching approach can recover key theoretical results like these gives some indication of its power.

We prove the result in the case $m=1$ for simplicity. The general result may be obtained by considering $S_{t}/\sqrt{m}$ in place of $S_{t}$, $V_{t}/m$ in place of $V_{t}$, and $c/\sqrt{m}$ in place of $c$. See Appendix F for details.

We first compute $\psi_{G}^{-1}(u)$ by taking the positive solution to the quadratic equation given by $\psi_{G}(\lambda)=u$, yielding

where we have used the identity $\sqrt{1+x}-1=\frac{x}{\sqrt{1+x}+1}$. Let

$K(u)$ will appear below. Now we start from the line-crossing inequality of Lemma 1: reparametrizing $r=\log\alpha^{-1}$, we have for any $r>0,\lambda>0$

We divide intrinsic time into epochs $\eta^{k}\leq V_{t}<\eta^{k+1}$ for each $k=0,1,\dots$, and we will construct a linear boundary over each epoch by carefully choosing values for $\lambda_{k}$ and $r_{k}$ and using the probability bound (38). We choose $\lambda_{k}$ so that the “standardized” boundary takes equal values at both endpoints of the epoch: $g_{\lambda_{k},r_{k}}(\eta^{k})/\eta^{k/2}=g_{\lambda_{k},r_{k}}(\eta^{k+1})/\eta^{(k+1)/2}$. This equation is solved by $\lambda_{k}=\psi_{G}^{-1}(r_{k}/\eta^{k+1/2})$, which yields, after some algebra,

Our goal, after choosing $r_{k}$ below, is to upper bound this expression by a function of $v$ alone, independent of $k$. Noting that the term in square brackets in (39) reaches its maximum over the $k$th epoch at the endpoints, $v=\eta^{k}$ and $v=\eta^{k+1}$, and substituting the expression (37) for $K(u)$, we have

The inequality $\eta^{k+1/2}\geq v/\sqrt{\eta}$ yields

for all $\eta^{k}\leq v<\eta^{k+1}$. This final expression no longer depends on $k$, showing that the final boundary $\mathcal{S}_{\alpha}(v)$ majorizes the corresponding linear boundary $g_{\lambda_{k},r_{k}}(v)$ over each epoch $\eta^{k}\leq v<\eta^{k+1}$ for $k=0,1,\dots$. Hence

But the first linear boundary $g_{\lambda_{0},t_{0}}(v)$ passes through $\mathcal{S}_{\alpha}(1)$ and has positive slope, which implies

Now taking a union bound over the probability bounds given by (38) for $k=0,1,\dots$, we have

Combining (46) with (45) proves that $v\mapsto\mathcal{S}_{\alpha}(1\vee v)$ is a sub-gamma uniform boundary with crossing probability $\alpha$.

For the second statement (9), we simply restrict the union bound to epochs $k\geq\lfloor\log_{\eta}V_{t}\rfloor$, which restricts the sum in (46) accordingly. ∎

We have given a stitched bound which is constant for $v<m$, but inspection of the proof shows that one may improve the bound to be linear with positive slope on $v<m$, by extending the linear bound over the first epoch to cover all $v>0$. This seems of limited utility for theoretical work, and we recommend other bounds over the stitched bound for practice, so we do not pursue this point further.

The idea of taking a union bound over geometrically spaced epochs is standard in the proof of the classical law of the iterated logarithm (Durrett, 2017, Theorem 8.5.1). The idea has been extended to finite-time bounds by Darling and Robbins (1967b), Jamieson et al. (2014), Kaufmann et al. (2016), and Zhao et al. (2016), usually when the observations are independent and sub-Gaussian; the technique is sometimes called “peeling”. Of course, Theorem 1 generalizes these constructions much beyond the independent sub-Gaussian case, but it also achieves tighter constants for the sub-Gaussian setting. Here, we briefly discuss how the improved constants arise.

Both Jamieson et al. (2014) and Zhao et al. (2016) construct a constant boundary rather than a linear increasing boundary over each epoch. They apply Doob’s maximal inequality for submartingales (Durrett, 2017, Theorem 4.4.2), as in Hoeffding (1963, eq. 2.17), to obtain boundaries similar to that of Freedman (1975). As illustrated in Howard et al. (2020, Figure 2), the linear bounds from Lemma 1 are stronger than corresponding Freedman-style bounds, and the additional flexibility yields tighter constants.

Both Darling and Robbins (1967b) and Kaufmann et al. (2016) use linear boundaries within each epoch analogous to those of Lemma 1. Both methods share a great deal in common with ours, and Darling and Robbins give consideration to general cumulant-generating functions. Recall from Lemma 1 that such linear boundaries may be chosen to optimize for some fixed time $V_{t}=m$. Our method chooses the linear boundary within each epoch to be optimal at the geometric center of the epoch, i.e., at $V_{t}=\eta^{k+{1/2}}$, so that at both epoch endpoints the boundary will be equally “loose”, that is, equal multiples of $\sqrt{V_{t}}$. Darling and Robbins choose the boundaries to be tangent at the start of the epoch, hence their boundary is looser than ours at the end of the epoch. Kaufmann et al. choose the boundary as we do, but appear to incur more looseness in the subsequent inequalities used to construct a smooth upper bound.

A.2 Proof of Corollary 1

Fix any $\epsilon>0$ and choose $a>0$ small enough that $\psi(\lambda)\leq(1+\epsilon)\lambda^{2}/2$ for all $\lambda\in(0,a)$. Using the fact that $\psi_{G,c}(\lambda)\geq\lambda^{2}/2$ for $c\geq 0$, we have $\psi(\lambda)\leq(1+\epsilon)\psi_{G,1/a}(\lambda)$ for all $\lambda\in(0,a)$, so that $(S_{t})$ is sub-gamma with scale $c=1/a$ and variance process $((1+\epsilon)V_{t})$. Now Theorem 1 shows that

where we may choose $u(v)\sim\sqrt{2(1+\epsilon)v\log\log v}$ (see (10) and discussion thereafter), so that $u((1+\epsilon)v)\sim\sqrt{2(1+\epsilon)^{2}v\log\log v}$. It follows that

As $\epsilon>0$ was arbitrary, we are done. ∎

A.3 Conjugate mixture proofs

We claim that $(M_{t})$ is a supermartingale with respect to $(\mathcal{F}_{t})$ on this augmented space. Indeed,

Now Definition 1 and Ville’s maximal inequality for nonnegative supermartingales (Durrett, 2017, exercise 4.8.2) yield

In the sub-Gaussian case, the following boundary is well-known (Robbins, 1970, example 2).

Suppose both $(S_{t})$ and $(-S_{t})$ are sub-Gaussian with variance process $(V_{t})$. Fix $\alpha\in(0,1)$ and $\rho>0$, and define

We have included the bound in Figures 9 and 4; although its $\mathcal{O}(\sqrt{V_{t}\log V_{t}})$ rate of growth is worse than the finite LIL discrete mixture bound, it can achieve tighter control over about three orders of magnitude of intrinsic time. This makes the normal mixture preferable in many practical situations when a sub-Gaussian assumption applies. When only a one-sided sub-Gaussian assumption holds, the normal mixture still yields a sub-Gaussian uniform boundary.

For any $\alpha\in(0,1)$ and $\rho>0$, the boundary

is a sub-Gaussian uniform boundary with crossing probability $\alpha$. Furthermore, we have the following closed-form upper bound:

The boundary $\textup{NM}_{\alpha}$ is easily evaluated to high precision by numerical root-finding, and the closed-form approximation is excellent: numerical calculations indicate that $\widetilde{\textup{NM}}_{0.025}(v)/\textup{NM}_{0.025}(v)<1.007$ uniformly when $\rho=1$, for example.

To obtain the explicit upper bound $\widetilde{\textup{NM}}_{\alpha}$ in (53) from the exact boundary (52), we use the inequality $1-\Phi(x)\leq e^{-x^{2}/2}$ for $x>0$, which follows from a standard Cramér-Chernoff bound. This implies

We set the RHS equal to $l_{0}/\alpha$ and solve to conclude

so long as $\textup{NM}_{\alpha}(v)>0$. But we are guaranteed that $\textup{NM}_{\alpha}(v)>0$, because the LHS of the inequality in (52) is increasing in $s$ on $s\geq 0$ and no larger than one when $s=0$, while the RHS $l_{0}/\alpha\geq 1$.

The fact that $\textup{NM}_{\alpha}$ is a sub-Gaussian uniform boundary follows directly from Lemma 2, and therefore $\widetilde{\textup{NM}}_{\alpha}$ is as well. ∎

Suppose $(S_{t})$ is sub-Bernoulli with variance process $(V_{t})$ and range parameters $g,h$, while $(-S_{t})$ is sub-Bernoulli with variance process $(V_{t})$ and range parameters $h,g$. Fix any $\rho>gh$, let $r=\rho-gh$, and define

As with the normal mixture, we have a one-sided variant as well.

Fix any $g,h>0$, $\alpha\in(0,1)$, and $\rho>gh$. Let $r=\rho-gh$ and define

Then $f_{g,h}$ is a sub-Bernoulli uniform boundary with crossing probability $\alpha$ and range parameters $g,h$.

In the sub-Bernoulli case, we first rewrite the exponential process $\exp\left\{\lambda S_{t}-\psi_{B}(\lambda)V_{t}\right\}$ in terms of the transformed parameter $p=[1+(h/g)e^{-\lambda}]^{-1}$. This is motivated by the transform from the canonical parameter to the mean parameter of a Bernoulli family, but keep in mind that we make no parametric assumption here, these are merely analytical manipulations. Then a truncated Beta distribution on $p\in[g/(g+h),1]$ yields the one-sided beta-binomial uniform boundary, while an untruncated mixture yields the two-sided boundary.

For simplicity of notation, we will assume here that the problem has been scaled so that $g+h=1$, e.g., by replacing $X_{t}$ with $X_{t}/(g+h)$. Using the sub-Bernoulli $\psi$ function $\psi_{B}(\lambda)=\frac{1}{gh}\log\left(ge^{h\lambda}+he^{-g\lambda}\right)$, the exponential integrand in our mixture is

after substituting the one-to-one transformation

followed by some algebra. We wish to integrate against a Beta mixture density on $p$ with parameters $r/h$ and $r/g$, which has mean $p=g$, corresponding to $\lambda=0$. For Proposition 8, we must also truncate to $\lambda\geq 0$, i.e., to $p\geq g$. The appropriately normalized mixture integral is then

The proof of Proposition 7 is nearly identical, but we integrate over the full Beta mixture rather than truncating.

To verify that our choice of $r$ ensures that $\lambda$ has approximate precision $\rho$ under the full (not truncated) mixture distribution, we use the delta method to calculate the approximate variance of $\lambda$ for large $r$ based on the variance of $p$ under the full Beta mixture:

Setting this equal to $1/\rho$ yields $r=\rho-gh$ as desired. ∎

Then $\textup{GE}_{\alpha}$ is a sub-exponential uniform boundary with crossing probability $\alpha$ for scale $c$.

The gamma-exponential mixture is the result of evaluating the mixture integral in (13) with mixing density

This is a gamma distribution with shape $\rho/c^{2}$ and scale $\rho/c$ applied to the transformed parameter $u=c^{-1}-\lambda$, truncated to the support $[0,c^{-1}]$. The distribution has mean zero and variance equal to $1/\rho$, making it comparable to the normal mixture distribution used above. As $\rho\to\infty$, the gamma mixture distribution converges to a normal distribution and concentrates about $\lambda=0$, the regime in which $\psi_{E}(\lambda)\sim\psi_{N}(\lambda)$, which gives some intuition for why the gamma-exponential mixture recovers the normal mixture when $\rho\gg c^{2}$.

Then the fact that $\textup{GM}_{\alpha}$ is a sub-exponential uniform boundary follows as a special case of Lemma 2.

Proving (67) is an exercise in calculus. Substituting the definition of $\psi_{E}$ and removing common terms, it suffices to show that

After change of variables $u=\left(\frac{cs+v+\rho}{c}\right)(c^{-1}-\lambda)$, the right-hand side is equal to

Now the definition of the regularized lower incomplete gamma function and a bit of algebra finishes the argument. ∎

Then $\textup{GP}_{\alpha}$ is a sub-Poisson uniform boundary with crossing probability $\alpha$ for scale $c$.

The proof follows the same contours as that of Proposition 8. Using the sub-Poisson $\psi$ function $\psi_{P}(\lambda)=c^{-2}(e^{c\lambda}-c\lambda-1)$, the exponential integrand in our mixture is

after substituting the one-to-one transformation $\theta=\theta(\lambda)\coloneqq e^{c\lambda}$, so that $\lambda=c^{-1}\log\theta$. We integrate against a gamma mixing distribution on $\theta$ with shape and scale parameters both equal to $\beta\coloneqq\rho/c^{2}$, truncated to $\theta\geq 1$, so that $\lambda\geq 0$:

This yields the closed-form mixture (72). To verify that our choice of $\beta$ ensures that $\lambda$ has approximate precision $\rho$ under the full (not truncated) mixture distribution, we use the delta method to calculate the approximate variance of $\lambda$ for large $\beta$ based on the variance of $\theta$ under the full gamma mixture:

A.4 Proof of Proposition 2

Under the conditions of Proposition 2, we have

Note that $m(s,v)$ is nondecreasing in $s$ and nonincreasing in $v$ (since $\psi\geq 0$ by our assumptions on $\psi$).

Choose $\delta\in(0,\lambda_{\max})$ so that $\psi$ has three continuous derivatives and $f$ is continuous and positive on $[0,\delta)$; such a value of $\delta$ must exist by conditions (i) and (ii). Before proving Proposition 2, we state several lemmas.

Under the conditions of Proposition 2, for any $b\in(0,\psi^{\prime}(\delta))$, we have $m(bv,v)<\infty$ and $m(bv,v)\to\infty$ as $v\to\infty$.

Now Laplace’s asymptotic approximation (Widder, 1942, Chapter VII.2, Theorem 2b) yields

where $C>0$ is a constant not depending on $v$. (The condition $b<\psi^{\prime}(\delta)$ ensures that the maximizer of $\lambda b-\psi(\lambda)$ lies within $[0,\delta)$.) Since the LHS of (78) lower bounds $m(bv,v)$ while the RHS diverges as $v\to\infty$, we must have $m(bv,v)\to\infty$ as $v\to\infty$. ∎

Under the conditions of Proposition 2, $m(\mathcal{M}_{\alpha}(v),v)=l_{0}/\alpha$ for all $v$ sufficiently large.

Let $\mathcal{C}(v)\coloneqq[0,\psi^{\prime}(\delta)v)$ for $v>0$. Lemma 4 shows that $m(s,v)<\infty$ for all $s\in\mathcal{C}(v)$. Since $m(s,v)$ is nondecreasing in $s$, we may apply dominated convergence to find that $s\mapsto m(s,v)$ is continuous at all $s\in\mathcal{C}(v)$. Condition (i) of Proposition 2 implies $\psi\geq 0$, so that $m(0,v)\leq 1\leq l_{0}/\alpha$ for all $v$. Finally, Lemma 4 shows that $\sup_{s\in\mathcal{C}(v)}m(s,v)\to\infty$ as $v\to\infty$. Hence, for $v$ sufficiently large, there exists $s\in\mathcal{C}(v)$ such that $m(s,v)>l_{0}/\alpha$.

We have argued that, for any sufficiently large $v$, $m(0,v)\leq l_{0}/\alpha<m(\bar{s},v)<\infty$ for some $\bar{s}<\psi^{\prime}(\delta)v$, and $m(\cdot,v)$ is continuous on $[0,\bar{s}]$. The conclusion follows from the definition (13) of $\mathcal{M}_{\alpha}$. ∎

Under the conditions of Proposition 2, $\lim_{v\to\infty}\mathcal{M}_{\alpha}(v)=\infty$.

Suppose for the sake of contradiction that there exists $a>0$ such that $\mathcal{M}_{\alpha}(v)\leq a$ for all $v$. Then, since $m(s,v)$ is nondecreasing in $s$, $m(\mathcal{M}_{\alpha}(v),v)\leq m(a,v)$ for all $v$. But for sufficiently large $v$, we can write $a=bv$ for some $b<\psi^{\prime}(\delta)$, so that Lemma 4 implies $m(a,v)<\infty$ for sufficiently large $v$. Since condition (i) of Proposition 2 implies $\psi\geq 0$, we have $m(s,v)$ is decreasing in $v$, and dominated convergence yields $m(a,v)\to 0$ as $v\to\infty$. But this implies $m(\mathcal{M}_{\alpha}(v),v)\to 0$, contradicting Lemma 5.

We have shown that $\limsup_{v\to\infty}\mathcal{M}_{\alpha}(v)=\infty$. But since $m(s,v)$ is nondecreasing in $s$ and nonincreasing in $v$, Lemma 5 implies $\mathcal{M}_{\alpha}(v)$ must be nondecreasing in $v$. It follows that $\lim_{v\to\infty}\mathcal{M}_{\alpha}(v)=\infty$. ∎

Under the conditions of Proposition 2, $\mathcal{M}_{\alpha}(v)=o(v)$.

Suppose for the sake of contradiction that $\mathcal{M}_{\alpha}(v)\geq bv$ for all $v$ sufficiently large, for some $0<b<\psi^{\prime}(\delta)$. Then (again using the fact that $m(s,v)$ is nondecreasing in $s$) $\lim_{v\to\infty}m(\mathcal{M}_{\alpha}(v),v)\geq\lim_{v\to\infty}m(bv,v)=\infty$ by Lemma 4, contradicting Lemma 5. ∎

We invoke Theorem 4 of Fulks (1951), setting Fulks’ $h$ equal to our $v$, Fulks’ $k$ equal to our $\mathcal{M}_{\alpha}(v)$, Fulks’ $\phi$ equal to our $\psi$, Fulks’ $\psi$ equal to the identity function, Fulks’ $f$ equal to our $f$, and Fulks’ $b$ equal to our $\lambda_{\max}$. Fulks’ assumptions (A1)-(A4) now read as follows.

requires $\psi(0)=\psi^{\prime}(0_{+})=0$, $\psi^{\prime\prime}(0_{+})>0$, $\psi$ has three continuous derivatives in a neighborhood of the origin, and $\psi$ is positive and nondecreasing on $(0,\lambda_{\max})$.

requires conditions on the identity function which are trivially satisfied.

requires $f$ to be integrable and to be continuous and positive at the origin.

requires $\mathcal{M}_{\alpha}(v)\to\infty$ as $v\to\infty$ and $\mathcal{M}_{\alpha}(v)=o(v)$.

(A1) and (A3) are satisfied by conditions (i) and (ii) of Proposition 2. (A4) is satisfied by Lemmas 6 and 7. For Fulks’ Theorem 4, it remains to verify that $\sqrt{v}=o(\mathcal{M}_{\alpha}(v))$. But if this were not true, then we could apply Theorem 1 or Theorem 2 of Fulks (1951) to conclude that $m(\mathcal{M}_{\alpha}(v),v)\to 0$ as $v\to\infty$, contradicting Lemma 5. So Fulks’ Theorem 4 yields

Using Lemma 5 to set $m(\mathcal{M}_{\alpha}(v),v)=l_{0}/\alpha$, we may write

which can be rearranged into the desired conclusion. ∎

We have proved the result for one-sided bounds, but a nearly-identical argument applies to two-sided bounds such as Proposition 7.

A.5 Proof of Theorem 2

Recall the discrete mixture support points and weights,

Figure 8 illustrates the construction. To see heuristically why the exponentially-spaced grid $\lambda_{k}=\mathcal{O}(\eta^{-k})$ makes sense, observe that the integrand $\exp\left\{\lambda s-\lambda^{2}v/2\right\}$ is a scaled normal density in $\lambda$ with mean $s/v$ and standard deviation $1/\sqrt{v}$. In the regime relevant to our curved boundaries, $s$ is of order $\sqrt{v}$, ignoring logarithmic factors. Hence the integrand at time $v$ has both center and spread of order $1/\sqrt{v}$, so as $v\to\infty$, the relevant scale of the integrand shrinks. With the grid $\lambda_{k}=\mathcal{O}(\eta^{-k})$ we have $\lambda_{k}-\lambda_{k+1}=\mathcal{O}(\lambda_{k})$, ensuring that the resolution of the grid around the peak of the integrand matches the scale of the integrand as $v\to\infty$.

The discrete mixture bound is a valid mixture boundary in its own right, based on a discrete mixing distribution, but we may wish to know how well it approximates the continuous-mixture boundary from which it is derived. To illustrate the accuracy of the discrete mixture construction, we compare it to the one-sided normal mixture bound, Proposition 6. By using the same half-normal mixing density in Theorem 2 and setting $\eta=1.05$, $\overline{\lambda}=100$, we may evaluate a corresponding discrete mixture bound $\textup{DM}_{\alpha}$. With $\rho=14.3$, $\alpha=0.05$ and $l_{0}=1$, numerical calculations indicate that $\textup{DM}_{\alpha}(v)/\textup{NM}_{\alpha}(v)\leq 1.004$ for $1\leq v\leq 10^{6}$, suggesting that Theorem 2 gives an excellent conservative approximation to the corresponding continuous mixture boundary over a large practical range. Of course, when a closed form is available as in Proposition 6, one should use it in practice. But an exact closed form integral is rarely available as it is in Proposition 6, and substantial looseness often accompanies closed-form approximations which provably maintain crossing probability guarantees. In such cases, unless a closed form is required, Theorem 2 is preferable. See figure 10 for an example; in this figure, the bounds of Balsubramani (2014) and Darling and Robbins (1968) involve closed-form mixture integral approximations.

A.6 Stitching as a discrete mixture approximation

Suppose we wish to analytically approximate the discrete mixture boundary $\textup{DM}_{\alpha}$ of Theorem 2 in the sub-Gaussian case $\psi=\psi_{N}$. Clearly the sum is lower bounded by the maximum summand, which gives

The last expression is the pointwise minimum of a collection of linear boundaries of the form presented in Lemma 1, each chosen with a different $\lambda_{k}$, and with nominal crossing rates $w_{k}\alpha$ so that a union bound over crossing events yields total crossing probability $\sum_{k}w_{k}\alpha\leq\alpha$. This is very similar to the stitching construction, with a slightly different choice of the sequence $\lambda_{k}$.

By equating $w_{k}$ from Theorem 2 with $1/h(k)$ from Theorem 1, this observation allows us to view a stitched bound with function $h(k)$ as an approximation to a mixture bound with mixture density $f(\lambda)=\Theta(1/\lambda h(\log\lambda^{-1}))$ as $\lambda\downarrow 0$. For exponential stitching, this yields $f(\lambda)=\Theta(1)$—densities approaching a nonzero constant as $\lambda\downarrow 0$, including the half-normal distribution, correspond to exponential stitched boundaries growing at a rate $\sqrt{V_{t}\log V_{t}}$. For polynomial stitching, we have the corresponding mixture density