A New Confidence Interval for the Mean of a Bounded Random Variable

TL;DR

This paper introduces a novel confidence interval for the mean of a bounded random variable that guarantees coverage without normality assumptions, applicable across all sample sizes and confidence levels.

Contribution

A new method for constructing confidence intervals that guarantees coverage for bounded variables without relying on normality assumptions.

Findings

Competitive with Student's t-intervals in accuracy

Guarantees coverage for all distributions on a bounded interval

Applicable for all sample sizes and confidence levels

Abstract

We present a new method for constructing a confidence interval for the mean of a bounded random variable from samples of the random variable. We conjecture that the confidence interval has guaranteed coverage, i.e., that it contains the mean with high probability for all distributions on a bounded interval, for all samples sizes, and for all confidence levels. This new method provides confidence intervals that are competitive with those produced using Student's t-statistic, but does not rely on normality assumptions. In particular, its only requirement is that the distribution be bounded on a known finite interval.

Equations

$$\mathrm{Pr}\left(m_{\alpha }(X_1,\dots ,X_n)\ge \mu \right)\ge 1-\alpha .$$

$$m(\mathbf{z},\mathbf{u})\stackrel{\underset{\mathrm{def}}{}}{=}1-\sum _{i=1}^n{u}_i(z_{i+1}-z_i),$$

$$Q(1-\alpha ,Y)\stackrel{\underset{\mathrm{def}}{}}{=}\inf \{y\in \mathbb{R}:F_Y(y)\ge 1-\alpha \},$$

$$m_{\alpha }(\mathbf{z})\stackrel{\underset{\mathrm{def}}{}}{=}Q(1-\alpha ,m(\mathbf{z},\mathbf{U})).$$

$$\mathrm{Pr}\left(m_{\alpha }(\mathbf{Z})\ge \mu \right)\ge 1-\alpha ,$$

$$(z_1,u_1),(z_2,u_2),...,(z_n,u_n).$$

$$F_{(\mathbf{z},\mathbf{u})}(x)\stackrel{\underset{\mathrm{def}}{}}{=}\{\begin{array}{ll}0,& \text{for }x<z_1\\ u_1,& \text{for }{z}_1\le x<z_2\\ u_2,& \text{for }{z}_2\le x<z_3\\ ...,& ...\\ u_n,& \text{for }{z}_n\le x<1\\ 1,& \text{for }x\ge 1.\end{array}$$

$$m(\mathbf{z},\mathbf{u})=\sum _{i=1}^{n+1}{z}_i(u_i-u_{i-1}),$$

$$m(\mathbf{z},\mathbf{u})=1-\sum _{i=1}^n{u}_i(z_{i+1}-z_i).$$

$$m_{\alpha }(\mathbf{z})\stackrel{\underset{\mathrm{def}}{}}{=}Q(1-\alpha ,m(\mathbf{z},\mathbf{U})).$$

$$\mathcal{Z}\overset{\underset{\mathrm{def}}{}}{=}\big{\{}\mathbf{z}=(z_{1},z_{2},...,z_{n}):0\leq z_{1}\leq z_{2}...\leq z_{n}\leq 1\big{\}}.$$

$$p_k=\frac{1-\mu }{1-k}.$$

$$[k,k,\dots ,k,k]$$

$$[k,k,\dots ,k,1]$$

$$[k,k,\dots ,1,1]$$

$$⋮$$

$$[k,1,\dots ,1,1]$$

$$[1,1,\dots ,1,1].$$

$${\mathbf{z}}_{j,n}\stackrel{\underset{\mathrm{def}}{}}{=}[\underset{j}{\underbrace{k,\dots ,k}},\underset{n-j}{\underbrace{1,\dots ,1}}].$$

$$\mathtt{If}\forall i,y_i\le z_i\mathtt{then}{m}_{\alpha }(\mathbf{y})\le m_{\alpha }(\mathbf{z}).$$

$$j_{\text{min}}(H_{k,\mu},1-\alpha,n)=\min\big{\{}j\in\{0,1,\dotsc,n\}:\mathbf{z}_{j,n}\in\overline{\mathcal{Z}}_{\alpha}^{\mu}\big{\}},$$

$$\underset{H_{k,\mu }}{\mathrm{Pr}}\left(\mathbf{Z}\in {\overline{\mathcal{Z}}}_{\alpha }^{\mu }\right)=$$

$$=$$

$$m({\mathbf{z}}_{j,n},\mathbf{U})$$

$$=$$

$$\underset{p_k:j_{\text{min}}=j}{\mathrm{argmax}}\underset{H_{k,\mu }}{\mathrm{Pr}}(\mathbf{Z}\in {\overline{\mathcal{Z}}}_{\alpha }^{\mu })\stackrel{\text{(a)}}{=}$$

$$\stackrel{\text{(b)}}{=}$$

$$\displaystyle\max\Big{\{}p_{k}\in[0,1]:j_{\text{min}}=j\Big{\}}$$

$$=$$

$$\stackrel{(a)}{=}$$

$$=$$

$$\stackrel{(b)}{=}$$

$$=$$

$$\stackrel{(c)}{=}$$

$$=$$

Full text

1 Introduction

Consider one of the fundamental problems in statistics: how to use $n$ samples of a real-valued random variable to obtain a confidence interval on its mean. Methods for constructing such confidence intervals are used across all branches of science. In the natural and social sciences, the confidence interval based on Student’s $t$-statistic (Student, 1908) is one of the standard tools for quantifying uncertainty about the results of an empirical study. In theoretical work, concentration inequalities like Hoeffding’s inequality (Hoeffding, 1963) are often used to analyze properties of algorithms in machine learning, data science, and other areas. Providing methods for obtaining tighter confidence intervals from fewer samples is critical to scientific advancement, enabling stronger conclusions to be drawn from the same experimental data.

In particular, there is a practical need today for confidence intervals that hold for small sample sizes. Since the confidence interval produced using Student’s $t$-statistic, which we refer to hereafter as the Student-$t$ interval, relies on the (near) normality of the sample mean, it is recommended that sample sizes be at least 30 for it to be used, unless there is a specific reason to believe that the population distribution is approximately normal. While other confidence intervals that hold for small sample sizes exist (such as Anderson’s (1969)), they produce intervals that are so wide as to be of little use in practice. This leaves the practitioner with the following choices:

• •

use methods, such as bootstrap methods, with no performance guarantees;

• •

use methods with unrealistic assumptions, such as the Student-$t$ interval;

• •

use valid but weak methods such as Hoeffding or Anderson’s inequalities that provide little information about the mean;

• •

abandon the idea of obtaining useful confidence intervals from the data.

In this paper, we introduce a new confidence interval for bounded distributions that is much tighter than other confidence intervals that come with guarantees. We conjecture that it holds for all bounded distributions, all samples sizes, and all confidence levels. We suggest that for many applications, this is the first practical confidence interval for sample sizes less than 30. We now offer a formal statement of the problem we are addressing.

We became aware of similar previous work after the initial publication of this paper (Gaffke, 2005). We thank Aaditya Ramdas for bringing it to our attention.

1.1 Problem Statement

Let $X_{1},\dotsc,X_{n}$ be $n$ independent and identically distributed real-valued random variables. Let each $X_{i}$ take values in the interval $[0,1]$ and have expected value $\mu$. For now we focus on a high-confidence upper bound, i.e., we desire a function $m_{\alpha}$ such that, for all distributions of $X_{i}$, all sample sizes $n\geq 1$, and all confidence levels $1-\alpha\in[0,1]$:

[TABLE]

That is, with probability at least $1-\alpha$, $m_{\alpha}(X_{1},\dotsc,X_{n})$ should be greater than or equal to the mean. Critically, in this statement the random quantity is the high-confidence upper bound, not the mean. Any definition of $m_{\alpha}$ that satisfies (1) for all bounded distributions, samples sizes $n\geq 1$, and $1-\alpha\in[0,1]$, is said to have guaranteed coverage.

In this paper we present a new method for constructing confidence intervals on the mean: a new $m_{\alpha}$. We conjecture that our function $m_{\alpha}$ satisfies (1), i.e., that it has guaranteed coverage. If our conjecture holds, this is the first confidence interval with tightness comparable to the Student-$t$ interval, but with guaranteed coverage in this setting. We prove in Section 7 that it dominates several other known confidence intervals with guaranteed coverage. That is, for every possible sample $(x_{1},x_{2},...,x_{n})$, it produces a confidence interval with width less than or equal to these previous methods, often with a much smaller width. This makes our confidence interval suitable for small sample sizes where other methods are not practical.

After defining $m_{\alpha}$ in the next section, we present the following results:

• •

a proof of (1) for a class of distributions that includes Bernoulli distributions, for all samples sizes $n$, and for all confidence levels $1-\alpha$;

• •

the sketch of a proof that our intervals are always at least as tight as those provided by Anderson (1969), which, in turn, are strictly tighter than those of Hoeffding (1963);

• •

results of extensive simulations on a wide variety of distributions that are consistent with (1) for many sample sizes and confidence intervals;

• •

empirical comparisons (through Monte Carlo simulations) with previous methods, demonstrating that the confidence intervals produced by $m_{\alpha}$ are consistently tighter than or as tight as the intervals produced by existing methods.

2 A New Confidence Interval for the Mean

In this section we present our new confidence interval. We also present our conjecture that it holds for all distributions bounded on $[0,1]$, for all sample sizes, and for all confidence levels.

Let $\mathbf{X}\overset{\underset{\mathrm{def}}{}}{=}(X_{1},X_{2},\dotsc,X_{n})$. Let $\mathbf{Z}\overset{\underset{\mathrm{def}}{}}{=}(Z_{1},Z_{2},\dotsc,Z_{n})$ be the order statistics of $\mathbf{X}$, i.e., $\mathbf{Z}$ is a vector containing the sorted values of $\mathbf{X}$ such that $Z_{1}\leq Z_{2}\leq\cdots\leq Z_{n}$. Let $\mathbf{z}\overset{\underset{\mathrm{def}}{}}{=}(z_{1},z_{2},...,z_{n})$ denote a particular sample of $\mathbf{Z}$ and $\mathbf{x}\overset{\underset{\mathrm{def}}{}}{=}(x_{1},x_{2},\dotsc,x_{n})$ a sample of $\mathbf{X}$. For notational convenience, we alternate between viewing $m_{\alpha}$ as a function of $\mathbf{z}$ or $\mathbf{x}$. So, when we write $m_{\alpha}(\mathbf{z})$ subsequently, this corresponds to a definition of $m_{\alpha}(\mathbf{x})$ where $\mathbf{z}$ are the order statistics of $\mathbf{x}$.

Let $\mathbf{U}$ be the order statistics of a sample of size $n$ from the continuous uniform distribution on $[0,1]$, with $\mathbf{u}\overset{\underset{\mathrm{def}}{}}{=}(u_{1},u_{2},\dotsc,u_{n})$ being a particular sample of $\mathbf{U}$. Since $\mathbf{u}$ are order statistics, $0\leq u_{1}\leq u_{2}\leq...\leq u_{n}\leq 1$. We define a function of two ordered vectors:

[TABLE]

where $z_{n+1}\overset{\underset{\mathrm{def}}{}}{=}1$. Let $Q(1-\alpha,Y)$ be the quantile function of the scalar random variable $Y$, i.e.,

[TABLE]

where $F_{Y}(y)$ is the cumulative distribution function (CDF) of $Y$.

Consider the random quantity $m(\mathbf{z},\mathbf{U})$, which depends upon a fixed sample $\mathbf{z}$ (non-random) and also on the random variable $\mathbf{U}$. We define $m_{\alpha}(\mathbf{z})$ to be the $(1-\alpha)$-quantile of $m(\mathbf{z},\mathbf{U})$, i.e.,

[TABLE]

We conjecture that this definition of $m_{\alpha}$ satisfies (1):

Conjecture 1**.**

Let $\mathbf{X}=(X_{1},\dotsc,X_{n})$ be $n$ independent and identically distributed random variables bounded in the interval $[0,1]$, each with mean $\mu$. Let $\mathbf{Z}=(Z_{1},\dotsc,Z_{n})$ be the order statistics of $\mathbf{X}$. Then for all $\alpha\in[0,1]$:**

$\Pr\left(m_{\alpha}(\mathbf{Z})\geq\mu\right)\geq 1-\alpha,$

(5)

where $m_{\alpha}$ is defined in (4).

Several extensions of this conjecture are apparent. First, since each $X_{i}$ is bounded above by $1$, this conjecture implies that $1-m_{\alpha}(1-\mathbf{z})$ is a $(1-\alpha)$-confidence lower bound on $\mu$. Second, if our main conjecture holds, we further conjecture that the assumption that the random variables are in $[0,1]$ can be extended to $(-\infty,1]$, or $[0,\infty)$ for the high-confidence lower bound. Furthermore, the deterministic upper bound of $1$ can be loosened to only require an almost-sure upper bound of $1$. Although these extensions may be important for some applications, hereafter we focus on the basic setting introduced previously.

3 Understanding $m_{\alpha}(\mathbf{z})$

Our high-confidence bound (for brevity, hereafter we refer to it as simply our bound) is given by the function $m_{\alpha}$ defined above. In this section we introduce the following concepts, which provide intuition for $m_{\alpha}$:

• •

ordered CDF pairs,

• •

the conservative completion of a set of ordered CDF pairs,

• •

the induced mean of a set of ordered CDF pairs, via conservative completion.

3.1 Ordered CDF pairs

For any order statistic vector $\mathbf{z}$, each element of $\mathbf{z}$ can be paired with an element from a non-decreasing sequence of numbers, $u_{1},u_{2},...,u_{n}$, to form $n$ pairs:

[TABLE]

Assuming the $u$’s are all in the interval $[0,1]$ (as is the case if $\mathbf{u}$ is a sample of $\mathbf{U}$), these pairs can be viewed as points on a CDF $F$, i.e., $u_{i}=F(z_{i})$. For this reason, we refer to these $n$ pairs as ordered CDF pairs, and write $(\mathbf{z},\mathbf{u})$ to denote such a set of ordered CDF pairs. We say that a set of ordered CDF pairs is consistent with a CDF $F$ if $u_{i}=F(z_{i})$ for all $i\in\{1,2,\dotsc,n\}$. Notice that a set of ordered CDF pairs is consistent with many (usually infinitely many) different CDFs—all non-decreasing functions on the interval $[0,1]$ that pass through these $n$ points (see Figure 1).

3.2 Conservative Completion of Ordered CDF Pairs

Given a set of ordered CDF pairs, one may ask which of the (usually infinitely many) CDFs that are consistent with the ordered CDF pairs represents the distribution with the greatest mean, and is this CDF unique? This CDF is unique, and we refer to it as the conservative completion of the ordered CDF pairs. That is, the mean of the distribution characterized by the conservative completion represents an upper bound on the mean of any distribution consistent with the set of ordered CDF pairs.

The conservative completion for a set of ordered CDF pairs, $(\mathbf{z},\mathbf{u})$, is illustrated in Figure 1. It is given by the CDF:

[TABLE]

3.3 The Induced Mean, $m(\mathbf{z},\mathbf{u})$

We introduce $m(\mathbf{z},\mathbf{u})$ to represent the mean of the distribution characterized by $F_{(\mathbf{z},\mathbf{u})}(x)$. This quantity is, for distributions over $[0,1]$, equivalent to the area of the region above the CDF, as depicted in Figure 2. The geometry of this figure suggests two methods for calculating this mean. The first, derived by decomposing the blue in Figure 2 into a set of horizontal strips, is

[TABLE]

where $u_{0}\overset{\underset{\mathrm{def}}{}}{=}0$, $u_{n+1}\overset{\underset{\mathrm{def}}{}}{=}1$, and $z_{n+1}\overset{\underset{\mathrm{def}}{}}{=}1$. Another formula is given by dividing the pink region into vertical strips, and is given by

[TABLE]

Next, we consider the distribution of such means obtained by allowing $\mathbf{u}$ to vary in a particular fashion.

3.4 A Distribution of Induced Means

Recall that $\mathbf{U}$ is a random vector containing $n$ samples from the continuous uniform distribution on $[0,1]$, sorted such that $0\leq U_{1}\leq U_{2}\leq...\leq U_{n}\leq 1$. We now consider a distribution of induced means obtained by replacing the fixed $\mathbf{u}$ in $m(\mathbf{z},\mathbf{u})$ with a random vector $\mathbf{U}$ to form a new scalar random variable $m(\mathbf{z},\mathbf{U})$.

Recall the definition of the quantile function from (3). We define $m_{\alpha}(\mathbf{z})$ to be the $(1-\alpha)$-quantile of the random variable $m(\mathbf{z},\mathbf{U})$, i.e.,

[TABLE]

Thus, $m_{\alpha}(\mathbf{z})$ considers the set of all $\mathbf{u}$’s that can be used to form an ordered CDF pair with a particular sample $\mathbf{z}$ and chooses the $(1-\alpha)$-quantile of the resulting induced means. As we shall see, this turns out to be just “conservative enough” to provide a valid high-confidence bound for Bernoulli-like distributions, and appears to be looser for distributions that are not Bernoulli-like.

4 The Order Statistic Simplex and Feasible Set

In this section, we define the order statistic simplex and the notion of feasible and infeasible sets of samples of order statistics. These definitions will be used in Section 5 to prove that our bound holds for all Bernoulli distributions and also for a more general set of Bernoulli-like distributions.

4.1 A Conditional Analysis

Our general method of proof (in the next section) will rely on a conditional analysis for a specific set of distributions. In particular, we will analyze our bound specifically for

• •

a fixed sample size $n$,

• •

a subset of the distributions with a specific mean, $\mu$,

• •

a specific confidence level, $1-\alpha$ (or equivalently, failure rate $\alpha$).

If we can show that the bound, (1), holds for each tuple $(n,\mu,1-\alpha)$, then we have a complete proof for the set of distributions under consideration.

Before proceeding with this method of proof, we define a few necessary terms.

4.2 The Order Statistic Simplex

Consider the order statistic simplex in $n$ dimensions—the set of all possible order statistic vectors, $\mathbf{z}$, which forms a polytope of dimension $n$ with $n+1$ vertices, i.e., a simplex. For distributions on $[0,1]$, we define the order statistic simplex as:

[TABLE]

The order statistic simplex for $n=2$ is depicted by the blue region in Figure 3.

4.3 The Infeasible Set of $\mathbf{z}$’s

Let the sample size, $n$, be fixed. For distributions with a specific mean, $\mu$, and for a specific confidence level, $1-\alpha$, let $\mathcal{Z}_{\alpha}^{\mu}\subseteq\mathcal{Z}$ denote the set of $\mathbf{z}$’s for which $m_{\alpha}(\mathbf{z})\geq\mu$. We refer to $\mathcal{Z}_{\alpha}^{\mu}$ as the feasible set, and say that $\mathbf{z}$ satisfies the bound if $\mathbf{z}\in\mathcal{Z}_{\alpha}^{\mu}$. Let $\overline{\mathcal{Z}}_{\alpha}^{\mu}$ denote the complement of this set, the set of $\mathbf{z}$’s for which $m_{\alpha}(\mathbf{z})<\mu$. We refer to $\overline{\mathcal{Z}}_{\alpha}^{\mu}$ as the infeasible set of $\mathbf{z}$’s, and say that $\mathbf{z}$ does not satisfy the bound if it is in $\overline{\mathcal{Z}}_{\alpha}^{\mu}$. Note that, given the mean, $\mu$, of a distribution, the feasible and infeasible sets have no other dependency on the unknown distribution of $\mathbf{X}$.

These ideas are illustrated in Figure 4. For sample size $n=2$, each plot shows the feasible (blue) and infeasible (green) regions for given confidence levels $1-\alpha$ and means $\mu$. Each row shows results for the same $\mu$ and different confidence levels. Note that in some cases, such as $1-\alpha=0.6$ and $\mu=0.3$, the entire order statistic simplex is feasible (there are no green pixels).

5 Bernoulli and Half-Bernoulli Distributions

In this section, we present a proof of our conjecture for Bernoulli distributions and for a generalization of Bernoullis that we refer to as half-Bernoulli distributions. While Bernoulli distributions have point masses on both 0 and 1, half-Bernoullis can have point masses at two positions: $k$ and $1$, where $0\leq k<1$. Thus, they are a generalization of Bernoullis that allow the lower value to be any non-negative value less than 1. Let $H_{k,\mu}$ be the half-Bernoulli distribution where the point masses are at $k$ and $1$, and the mean is $\mu$. With $k$ and $\mu$ specified by $H_{k,\mu}$, the probability $p_{k}$ of sampling $k$ is given by

[TABLE]

Bernoullis (and half-Bernoullis) are prime candidates for distributions for which the bound will fail (violate (1)), since it is not uncommon to have a sample of all 0’s (or all $k$’s), despite having a relatively large mean. For example, the probability of obtaining a sample of $[0,0,0,0]$ from a Bernoulli distribution with parameter $p=0.5$ is $0.5^{4}=0.0625$. Since the probability of getting this sample is greater than $0.05$, the bound must produce a result greater than $\mu=0.5$ for this sample at the confidence level $1-\alpha=0.95$. As we demonstrate below, the bound holds for all Bernoulli and half-Bernoulli distributions.

5.1 Finding Worst Case Half-Bernoullis for $(n,\mu,1-\alpha)$

Our approach will be to find “worst case” distributions among the set of half-Bernoulli distributions. By worst case, we mean that the probability of drawing a sample $\mathbf{z}$ for which the bound fails is as high as possible. We will show that for the worst case half-Bernoulli distributions (there can be more than one of these for each tuple $(n,\mu,1-\alpha)$, the probability of drawing a sample $\mathbf{z}\in\overline{\mathcal{Z}}_{\alpha}^{\mu}$ is no more than $\alpha$. Since the bound holds for the worst case half-Bernoullis, it must hold for all half-Bernoullis.

5.1.1 Enumerating possible $\mathbf{z}$’s

Consider a half-Bernoulli distribution with probability masses at $k$ and $1$, and with mean $\mu$. For a given $n$, there are $n+1$ possible order statistics, $\mathbf{z}$:

[TABLE]

Let $\mathbf{z}_{j,n}$ be the sample $\mathbf{z}$ with $j$ out of $n$ values of $k$, i.e., for $j\in\{0,1,\dotsc,n\}$:

[TABLE]

5.1.2 Monotonicity of $m_{\alpha}(\mathbf{z})$

Notice that $m_{\alpha}$ is monotonic in the following sense. For two samples $\mathbf{y}$ and $\mathbf{z}$:

[TABLE]

It follows from (20) that $m_{\alpha}(\mathbf{z}_{i,n})\leq m_{\alpha}(\mathbf{z}_{j,n})$ whenever $i<j$.

5.1.3 An expression for $\Pr(\mathbf{Z}\in\overline{\mathcal{Z}}_{\alpha}^{\mu})$

For any half-Bernoulli distribution $H_{k,\mu}$, sample size $n$, and confidence level $1-\alpha$, let

[TABLE]

where $j_{\text{min}}(H_{k,\mu},1-\alpha,n)=n+1$ if $\mathbf{z}_{n,n}\in\mathcal{Z}_{\alpha}^{\mu}$.

For example, suppose $n=5$ and $[k,k,k,1,1]\in\overline{\mathcal{Z}}_{\alpha}^{\mu}$ but $[k,k,1,1,1]\in\mathcal{Z}_{\alpha}^{\mu}$. Then $j_{\text{min}}=3$, where here and in the following the arguments of $j_{\text{min}}$ are implicit. By the monotonicity of the bound (see (20)), all of the samples with ${\tt count}(k)\geq j_{\text{min}}$ will be in $\overline{\mathcal{Z}}_{\alpha}^{\mu}$.

For a given half-Bernoulli distribution we can now write an expression for the probability that $\mathbf{Z}$ will not satisfy the bound:

[TABLE]

where $\beta_{\text{cdf}}(x;a,b)$ is the CDF of a beta distribution with parameters $a$ and $b$. The above derivation uses the property that each $\mathbf{z}_{i,n}$ can be viewed as a sample from a binomial distribution, and in the last step we use a well-known identity that relates the sum of binomials to the CDF of a beta distribution.

5.1.4 Simplification of $m(\mathbf{z}_{j,n})$ due to the simple structure of $\mathbf{z}_{j,n}$

Before continuing with deriving the $p_{k}$ that maximizes the failure rate of the bound, we show how $m(\mathbf{z},\mathbf{U})$ simplifies for samples from half-Bernoulli distributions.

Recall that our bound is a quantile of the function $m(\mathbf{z},\mathbf{U})$ with respect to the uniform random variable $\mathbf{U}$. For samples of the form $\mathbf{z}_{j,n}$, this function reduces to a particularly simple form:

[TABLE]

That is, with the exception of the $j^{\text{th}}$ term, all of the successive differences of $\mathbf{z}_{j,n}$ are 0,111We can ignore the case where $j=n$, since the bound is trivial for $\mathbf{z}_{n,n}=(1,1,...,1)$. leaving us with a simple function of the $j^{\text{th}}$ order statistic, $U_{j}$. Later it will be useful to note that the $j^{\text{th}}$ order statistic when taking $n$ samples from a uniform distribution is beta distributed with parameters $j$ and $n-j+1$ (Casella and Berger, 2002, Example 5.4.5).

5.1.5 Choosing $p_{k}$ to maximize the failure rate

For a fixed $(n,\mu,1-\alpha)$, consider the set of distributions, $H_{k,\mu}$, with $j_{\text{min}}=j$, for some value $j$. We are interested in the $H_{k,\mu}$ that maximizes $\Pr(\mathbf{Z}\in\overline{\mathcal{Z}}_{\alpha}^{\mu})$. Since we are only considering half-Bernoulli distributions with a particular mean, $\mu$, the entire distribution is specified if $p_{k}$ is specified, and so we solve for the $p_{k}$ that maximizes $\Pr(\mathbf{Z}\in\overline{\mathcal{Z}}_{\alpha}^{\mu})$:

[TABLE]

Step (a) follows from (24). Step (b) follows since all beta CDFs are monotonic in their first argument. In other words, within the set of $H_{k,\mu}$ that have the same $\mu$ and $j_{\text{min}}$, the failure rate of the bound (the probability that $\mathbf{Z}\in\overline{\mathcal{Z}}_{\alpha}^{\mu}$) is monotonic in $p_{k}$.

Although the failure rate is monotonic in $p_{k}$, this does not mean that the worst-case is when $p_{k}=1$, since this monotonicity result is restricted to the set of half-Bernoulli distributions with $j_{\text{min}}=j$. We therefore now solve for the maximum $p_{k}$ such that $j_{\text{min}}=j$ in order to obtain the half-Bernoulli distribution with mean $\mu$ that maximizes the failure rate:

[TABLE]

Step (a) follows from the monotonicity of the bound and Section 5.1.2. Step (b) uses the result of (27). Step (c) uses the the fact that the $j^{\text{th}}$ order statistic of a uniform sample of size $n$ is beta distributed with parameters $j$ and $n-j+1$.

5.1.6 Bringing the pieces together

We have established that, for a given $n$, of all half-Bernoulli distributions with mean $\mu$ and $j_{\text{min}}=j$, the one that maximizes the failure rate of the bound has $p_{k}=\beta_{\text{cdf}}^{-1}\left(\alpha;j,n-j+1\right)$. Plugging this into (24), we have that

[TABLE]

Thus, we have seen that by maximizing the probability that $\mathbf{Z}$ causes $m_{\alpha}(\mathbf{Z})>\mu$, i.e., maximizing $\Pr(\mathbf{Z}\in\overline{\mathcal{Z}}_{\alpha}^{\mu})$, we can produce a probability of violation of at most $\alpha$. Thus, the bound holds with probability at least $1-\alpha$. Since this is true for all half-Bernoulli distributions for all values of $j_{\text{min}}$ and arbitrary tuples $(n,\mu,1-\alpha)$, it is true for all half-Bernoullis, all sample sizes and all confidence levels.

6 Computing $m_{\alpha}$

In this section, we discuss two methods for computing our bound, $m_{\alpha}(\mathbf{z})$, for a particular sample $\mathbf{z}$. The first is based upon a geometric analysis of the bound and the second uses a Monte Carlo sampling technique.

6.1 Geometric Computation of the Bound

Recall that the random variable $\mathbf{U}$ represents the order statistics of a uniform sample, and hence lies in the order statistic simplex defined in Section 4.2. Figure 5 shows the order statistic simplex for $n=2$. Note that the order statistic simplex of dimension $n$ has volume $\frac{1}{n!}$.

Let $\mathbf{s}\overset{\underset{\mathrm{def}}{}}{=}[z_{2}-z_{1},z_{3}-z_{2},...,z_{n+1}-z_{n}]$ be the spacings of the sample. Consider an example in which $\mathbf{z}=[0.3,0.8]$. Then $\mathbf{s}=[0.5,0.2]$, as shown in the figure.

Starting from the definition of our bound $m_{\alpha}(\mathbf{z})$ and expanding the definition of the quantile function, we have

[TABLE]

This final expression has a clear geometric interpretation. The set of points $\mathbf{u}$ such that $\mathbf{u}\cdot\mathbf{s}$ is greater than some value is the upper right region of the order statistic simplex, depicted by the pink region in Figure 5. The bound is defined to be the least value of $\hat{\mu}$ such that the volume of the pink region is a fraction $1-\alpha$ of the simplex volume.

This value of $\hat{\mu}$ can be found by evaluating the volume of the section of the order statistic simplex above the hyperplane $l$, a hyperplane orthogonal to the spacings vector $\mathbf{s}$. Thus, we seek the smallest value of $\hat{\mu}$ such that the volume of this section is $100(1-\alpha)\%$ of the volume of the order statistic simplex. This value of $\hat{\mu}$ is our bound.

Closed-form expressions for sections of simplexes cut by hyperplanes have been published by several authors, including Lasserre (2015). These expressions lead to efficient calculations of the bound in most cases. However, these formulas have singularities that cause problems for certain samples $\mathbf{z}$, such as samples with repeated values. Thus, we explore a more reliable Monte Carlo approach for computing our bound below.

6.2 Monte Carlo Estimate of the Bound

Since the bound is defined in terms of a quantile of a function that depends upon a uniform random variable, it is simple to develop a Monte Carlo estimate. This is provided in Algorithm 1.

Algorithm 1 can be implemented more efficiently if $m_{\alpha}$ will be estimated multiple times for the same $n$, since samples of $\mathbf{U}$ can be computed and sorted a single time. Also, notice that the number of Monte Carlo samples, $l$, does not scale poorly for distributions with rare values, since we are estimating the $(1-\alpha)$-quantile of $m(\mathbf{z},\mathbf{U})$, which is robust to outliers (e.g., $\alpha=0.5$ makes this the median, which is well known to be robust to outliers).

In practice, we find that $l=10,\!000$ tends to provide a reasonable approximation of $m_{\alpha}(\mathbf{x})$ for $\alpha=0.05$. Note that (1) may not hold when using this Monte Carlo estimate of $m_{\alpha}$ due to error in the estimate. In practice this can be remedied by increasing $l$, or by incorporating high-probability bounds on the error in the Monte Carlo estimate into the bound. Also, note that as $\alpha$ decreases, $l$ should be increased.

7 Related Work

In this section, we review other methods for computing high-confidence upper bounds on the mean of a random variable from samples (several of which we compare to in the subsequent numerical analysis section). Although some of these methods extend to more general settings (e.g., Hoeffding’s inequality does not require identically distributed samples, and Anderson’s inequality does not require a lower-bound on the random variable), here we consider only the standard setting that we have discussed in this paper, wherein the samples are i.i.d. and the random variable always takes values in the interval $[0,1]$. We divide this section into two parts: prior methods that provide guaranteed coverage, and prior methods that do not provide guaranteed coverage. We present these prior methods as functions, $m_{\alpha}^{\text{Hoeffding}}$, $m_{\alpha}^{\text{Maurer\&Pontil}}$, etc., each of which provides an alternative to $m_{\alpha}$.

7.1 Prior Methods with Guaranteed Coverage

The methods presented in this subsection have guaranteed coverage in the setting that we have described—they satisfy (1) if used in place of $m_{\alpha}$.

Using Hoeffding’s inequality (Hoeffding, 1963) to construct a high-confidence upper-bound on $\mu$ is perhaps the best known, and simplest, prior method with guaranteed coverage:

[TABLE]

where $\bar{\mathbf{x}}$ is the sample mean, i.e., $\bar{\mathbf{x}}\overset{\underset{\mathrm{def}}{}}{=}\frac{1}{n}\sum_{i=1}^{n}x_{i}$. In cases where the variance of the random variable is significantly less than one, the upper bounds provided by Maurer and Pontil’s empirical Bernstein bound (Maurer and Pontil, 2009) can be tighter than those produced by Hoeffding’s inequality. This is achieved by leveraging not just the sample mean, $\bar{\mathbf{x}}$, but also the sample variance, $\widehat{\operatorname{Var}}(\mathbf{x})\overset{\underset{\mathrm{def}}{}}{=}\frac{1}{n-1}\sum_{i=1}^{n}(x_{i}-\bar{\mathbf{x}})^{2}$:

[TABLE]

Going one step further, Anderson’s inequality provides high-confidence upper bounds on the mean by using the entire sample CDF (rather than only the sample mean and variance):

[TABLE]

where for $i\in\{1,2,\dotsc,n\}$,

[TABLE]

is a vector that Anderson derived from the Dvoretsky-Kiefer-Wolfowitz (DKW) inequality (Dvoretzky et al., 1956). Note that the form we present for Anderson’s inequality uses the tight constant for the DKW inequality found by Massart (1990), which relies on the assumption that $\alpha\leq 0.5$. This restriction is not restrictive because high-confidence bounds are typically applied with small values of $\alpha$, e.g., $\alpha=0.05$.

Following these three methods, several alternatives have been proposed, including other approaches that rely only on the sample mean (Chen, 2008), methods that extend Maurer and Pontil’s empirical Bernstein bound to provide tighter bounds for random variables with long tails (Bubeck et al., 2012; Thomas et al., 2015a), and methods that provide alternatives to Anderson’s inequality that use alternative methods of defining inclusion envelopes for a distribution’s CDF (Learned-Miller and DeStefano, 2008; Diouf and Dufour, 2005).

7.2 Prior Methods without Guaranteed Coverage

All of the methods presented in this subsection do not have guaranteed coverage in the setting that we have described, but are often used to compute high-confidence upper bounds on the mean.

Perhaps the most common method for constructing high-confidence upper bounds is based on Student’s $t$-statistic (Student, 1908):

[TABLE]

where $t_{1-\alpha,\nu}$ denotes the $100(1-\alpha)$ percentile of the Student’s $t$ distribution with $\nu$ degrees of freedom. We refer to this confidence interval as the Student-$t$ interval. If $\bar{\mathbf{X}}$ is normally distributed, then $m_{\alpha}^{\text{Student}}$ does provide guaranteed coverage. The central limit theorem implies that $\bar{\mathbf{X}}$ tends towards a normal distribution as $n$ increases, and so this method is often applied in scientific research if $n\geq 30$, even though this does not provide coverage guarantees.

Bootstrap methods tend to provide the tightest confidence intervals for the mean. However, this comes at a large cost: they do not have guaranteed coverage, even with normality assumptions. Despite concerns about their reliability, bootstrap methods remain in common use due to their tight confidence intervals and tendency to produce error rates roughly around $\alpha$ for many common distributions (Hanna et al., 2017; Thomas et al., 2015b). Since bootstrap methods are not easily expressed as closed-form alternatives to $m_{\alpha}$, we refer the reader to the work of Efron and Tibshirani (1993) for details on these approaches. The two that we focus on in our subsequent experiments are the most common, the percentile bootstrap, and one of the most sophisticated, the bias corrected and accelerated (BCa) bootstrap.

One limitation of BCa, the more sophisticated bootstrap method, is that it is not defined in some cases (e.g., if all of the samples take the same value) and can encounter numerical issues in other cases. In our implementations, whenever numerical issues are detected, the method automatically reverts to the percentile bootstrap.

8 Theoretical Analysis

In this section, we provide an analytic comparison of our bound to two prior methods that provide guaranteed coverage: Hoeffding’s inequality and Anderson’s inequality. We will show that, for any sample $\mathbf{z}$, our high-confidence upper bound is less than that resulting from Hoeffding’s inequality, and never greater than that of Anderson’s inequality. We break this result into two components. First we show that for all $\mathbf{z}$, $m_{\alpha}(\mathbf{z})\leq m_{\alpha}^{\text{Anderson}}(\mathbf{z})$. Second, we show that for all $\mathbf{z}$, $m_{\alpha}^{\text{Anderson}}(\mathbf{z})\leq m_{\alpha}^{\text{Hoeffding}}(\mathbf{z})$, which implies that $m_{\alpha}(\mathbf{z})\leq m_{\alpha}^{\text{Hoeffding}}(\mathbf{z})$, where these inequalities are strict if $\alpha\leq 0.5$.

8.1 Theoretical Comparison to Anderson’s Inequality

In this section we compare $m_{\alpha}$ to $m_{\alpha}^{\text{Anderson}}$.

Theorem 1**.**

For all possible values $\mathbf{z}$ of $\mathbf{Z}$ and all $\alpha\in[0,0.5]$,

[TABLE]

Proof.

We present a sketch of the proof. Consider the diagram in Figure 6. This figure depicts, for $n=2$, the space of possible vectors $\mathbf{u}$, which are sorted uniform samples. The point $(u_{1},u_{2})$ represents $\mathbf{u}^{\text{DKW}}$, defined in (52). $\mathcal{U}_{\text{DKW}}$ denotes the set of vectors that are element-wise greater than $\mathbf{u}^{\text{DKW}}$. It follows from the DKW inequality, with the tight constants found by Massart (1990), that the probability $\mathbf{U}$ is in $\mathcal{U}_{\text{DKW}}$ is at least $1-\alpha$. The region $\mathcal{U}_{\text{LMT}}$ is any set of $\mathbf{u}$’s that result in the lowest induced means, $m(\mathbf{z},\mathbf{u})$, while ensuring that the probability that $\mathbf{U}$ is in $\mathcal{U}_{\text{LMT}}$ is precisely $1-\alpha$. Note that any point that is not contained within the pink region must represent a vector $\mathbf{u}$ that results in an induced mean, $m(\mathbf{z},\mathbf{u})$, which is greater than the induced mean of any point in the pink region. Our bound is effectively the maximum over induced means of points in the pink region and Anderson’s bound is the maximum over points in the blue region. Since the probability that $\mathbf{U}$ is in $\mathcal{U}_{\text{DKW}}$ cannot be less than the probability that $\mathbf{U}$ is in $\mathcal{U}_{\text{LMT}}$, and $\mathcal{U}_{\text{LMT}}$ contains the $\mathbf{u}$ vectors that minimize the induced mean, our bound cannot be larger than Anderson’s. ∎

8.2 Analytic Comparison to Hoeffding’s Inequality

In this section we prove the following theorem:

Theorem 2**.**

For all possible values $\mathbf{z}$ of $\mathbf{Z}$ and all $\alpha\in[0,1]$,

[TABLE]

where the inequality is strict if $\alpha\leq 0.5$.

Proof.

We begin with $m_{\alpha}^{\text{Anderson}}(\mathbf{z})$ and present a sequence of inequalities that conclude with $m_{\alpha}^{\text{Hoeffding}}(\mathbf{z})$, where one inequality is strict if $\alpha\leq 0.5$:

[TABLE]

If $\alpha\leq 0.5$, then this final inequality is strict because, when $i=1$, we have that for any $n$, $0>i/n-\sqrt{\ln(1/\alpha)/2n}$, and so

[TABLE]

Continuing, we have:

[TABLE]

Notice that (64) provides an expression similar to Hoeffding’s inequality, but where the lower bound on the random variable (in our case, zero) is replaced by the smallest observed sample, $z_{1}$. This presents a tighter variant of Hoeffding’s inequality that holds when $\alpha\leq 0.5$ and the random variables are i.i.d. (the general form of Hoeffding’s inequality holds for random variables that are not necessarily identically distributed). ∎

It then follows from Theorem 2 that our bound is always at least as tight as Hoeffding’s inequality, and is strictly tighter if $\alpha\leq 0.5$:

Corollary 1**.**

For all possible values $\mathbf{z}$ of $\mathbf{Z}$ and all $\alpha\in[0,1]$,

[TABLE]

where the inequality is strict if $\alpha\leq 0.5$.

Proof.

This follows immediately from Theorems 1 and 2. ∎

9 Numerical Analysis

In this section we present results from a numerical analysis of our bound. These empirical results aim to answer the following research questions:

1. RQ1

For a variety of distributions, confidence levels, and number of samples, are results consistent with (1)? 2. RQ2

For a variety of distributions that resemble common use-cases, how do the confidence intervals produced by our bound compare to those of previous methods that have guaranteed coverage (i.e., those that satisfy (1))? 3. RQ3

This question is the same as RQ2, but for methods that do not have guaranteed coverage. 4. RQ4

Can our bound provide confidence intervals that are practical for scientific experiments with fewer than $30$ samples?

9.1 Numerical Studies on Guaranteed Coverage

In this subsection we study RQ1 with experiments that are consistent with the conjecture that $m_{\alpha}$, as we have defined it, has guaranteed coverage (satisfies (1)). Although they show that (1) appears to hold for a variety of settings, this does not imply that settings do not exist under which (1) does not hold.

To study RQ1, we selected a variety of different distributions (uniform, beta, and Bernoulli, each with various parameters), confidence levels $1-\alpha$, and number of samples, $n$. For each such tuple, $(\text{distribution},1-\alpha,n)$, we collected $10,\!000$ samples of $\mathbf{Z}$, computed $m_{\alpha}(\mathbf{Z})$, and checked whether $m_{\alpha}(\mathbf{Z})\geq\mu$. From these $10,\!000$ tests, we estimated the coverage—the probability that the bound holds. That is, we estimated $\Pr(m_{\alpha}(\mathbf{Z})\geq\mu)$ by dividing the number of samples of $\mathbf{Z}$ such that $m_{\alpha}(\mathbf{Z})\geq\mu$ by $10,\!000$.

Although or goal here is to study RQ1, to facilitate the interpretation of the presented results, we provide results using two prior methods that exhibit the different types of behavior that our bound might produce. This comparison also provides some insight into RQ2 (via the comparison to Heoffding’s inequality) and RQ3 (via comparison to the Student-$t$ interval). However, note that subsequent experiments further study these two research questions.

First consider Figure 7, which presents results using Hoeffding’s inequality. The top left plot shows the coverage for a variety of beta distributions, but with $n$ fixed. To interpret this plot, consider the curve $\beta(1,5)\,\,n=10$ at the $0.7$ position on the horizontal axis. The position $0.7$ on the horizontal axis indicates that we requested an upper bound that holds with probability at least $0.7$. Since, at horizontal position $0.7$, the curve $\beta(1,5)\,\,n=10$ (red curve) lies above $0.7$ (blue curve), the upper bound produced by Hoeffding’s inequality held with probability greater than $0.7$. Hence, a curve remaining above the blue line indicates that the desired confidence level was achieved. However, notice that the curve is far above the blue line—this indicates that Hoeffding’s inequality was overly conservative. It provided a high-confidence upper bound that was greater than or equal to the mean far more often than requested. This means that for this distribution the confidence interval provided by Hoeffding’s inequality is not tight, and could be improved.

The three other plots in Figure 7 are similar, but use different parameters. The top row presents results for beta distributions, while the bottom row presents results for Bernoulli distributions. The left column presents results as the parameters of the distributions are varied, while the right column presents results as the number of samples is varied. Overall Figure 7 shows the behavior that we would expect of a bound that has guaranteed coverage, but which provides loose high-confidence upper bounds.

Now consider Figure 8, which is identical to Figure 7, except that it uses the Student-$t$ interval instead of Hoeffding’s inequality. This plots shows very different behavior: the curves tend to be much lower, indicating tighter confidence intervals around the sample mean. However, the curves often cross the blue line, indicating that in these settings the Student-$t$ interval does not provide guaranteed coverage—if you ask for an $0.8$-confidence upper bound, you may only get a $0.6$-confidence upper bound. Hence Figure 8 shows the behavior that we would expect of a bound that does not have guaranteed coverage, but which provides tight “high-confidence” upper bounds.

The desired behavior of a high-confidence upper bound would blend the desirable properties of Hoeffding’s inequality and the Student-$t$ interval. In these plots, this would result in curves that always remain above the blue curve (guaranteed coverage), but are otherwise as low as possible (tight). Figure 9 presents the results of this same experiment, conducted using our bound. It achieves this desired behavior—it always remains above the blue line (consistent with guaranteed coverage), but tends to be significantly lower than Hoeffding’s inequality.

9.2 Numerical Comparison to Previous Methods

In this subsection we focus on RQ2 and RQ3 with experiments that compare the tightness of our bound to that of previous methods (both with and without guaranteed coverage). A variety of different statistics can be used to capture how tight the high-confidence upper bounds produced by a method are, including the mean upper bound, i.e., $\mathbf{E}[m_{\alpha}(\mathbf{Z})]$, and the median upper bound. Here we report the mean upper bound: we gather $1,\!000$ samples of $\mathbf{Z}$ from a distribution and compute the upper bounds produced by our bound and several previous methods and report the sample mean of the upper bounds for each method. For simplicity, here we vary the distribution and $n$ but fix $\alpha=0.05$ to obtain $95\%$-confidence upper bounds.

First compare the blue curve (our bound) to the black curves (previous methods with guaranteed coverage), noting that the horizontal axis uses a logarithmic scale. In every case, the blue curve remains strictly below the black curves, indicating that in every setting our bound produces lower values on average. Notice that frequently our bound obtains mean upper bounds that previous methods require an order of magnitude more samples to achieve, indicating that our bound is a drastic improvement in tightness and/or data efficiency.

Next, compare the blue curve (our bound) to the red curves (previous methods that do not have guaranteed coverage). The two bootstrap methods do not provide guaranteed coverage, even with normality assumptions. So, although they produce tight confidence intervals (as is evident in these plots), the high-confidence bounds that they produce cannot be relied upon.

Next consider the Student-$t$ interval: for the uniform distribution, it produces high-confidence upper bounds that are similar to those produced by our bound. When the Student-$t$ interval is computed from normally distributed data, it produces a tighter high-confidence upper-bound than our bound. However, when the sampling distribution includes right-skew (e.g., $\beta(1,10)$), the Student-$t$ interval tends to be overly optimistic—it does not have guaranteed coverage (as is evident in Figure 8 with $\beta(1,5)$). Hence, although the upper bound of the Student-$t$ interval tends to be lower than those produced by our method for $\beta(1,10)$, it does not have guaranteed coverage. On the other hand, when the sampling distribution includes left-skew (e.g., $\beta(10,1)$, the Student-$t$ interval is overly-conservative (like Hoeffding’s inequality). Hence, in Figure 10, the $\beta(10,1)$ plot indicates that our bound is tighter than the Student-$t$ interval. This is further evidence that our bound is combining the desirable properties of Hoeffding’s inequality and the Student-$t$ interval: it roughly preserves the tightness of the Student-$t$ interval, except in the cases where the Student-$t$ interval is too tight to provide guaranteed coverage (in which case our bound is sufficiently looser to provide guaranteed coverage).

9.3 Numerical Support for Practical use With $n<30$

Of the many potential uses of our bound, one stands out: it provides a valid method for constructing confidence intervals for scientific studies with fewer than $30$ samples. Even though the Student-$t$ interval does not have guaranteed coverage when the sampling distribution is not normal, the central limit theorem tells us that the sample mean tends towards a normal distribution as $n$ increases. Hence, the Student-$t$ interval becomes reasonable when $n$ is large.222Notice that even with arbitrarily large $n$, the Student-$t$ interval may not have guaranteed coverage, so here saying that the Student-$t$ interval is “reasonable” does not mean that it has guaranteed coverage. However, without knowing the sampling distribution, it is not clear how large $n$ must be for the Student-$t$ interval to be reasonable. A common rule of thumb used in current scientific research is that $n$ must be at least $30$.

This raises the question: what should one do when fewer than $30$ samples are available? Our bound provides an answer (assuming our conjecture is true), as it provides confidence intervals of comparable tightness, but with guaranteed support for any $n$ and without any normality assumptions. The only requirement is the ability to identify limits on the support of the distribution. To answer RQ3, we present an experiment that shows how our bound can be used to obtain confidence intervals based on fewer than $30$ empirical measurements.

Specifically, we used data from the United States Census from the year 2000 to obtain an estimate of the distribution of people’s ages, considering only people zero to 84 years old. We then consider the problem of obtaining a tight high-confidence upper bound on the mean age of people ages zero to 84 based on $n<30$ samples. The results of this experiment are presented in Figure 11, which is a similar form to Figure 10 (but without the logarithmic horizontal axis). The key observations from this plot are: 1) previous methods with guaranteed coverage are too loose to provide useful high-confidence bounds with so few samples, 2) the Student-$t$ interval is sufficiently tight, but it cannot be applied responsibly with such a small $n$, and 3) our bound produces high-confidence bounds that are comparable to the Student-$t$ interval (while maintaining guaranteed coverage, if our conjecture holds).

10 Acknowledgements

This work benefitted significantly from conversations with Vince Lysinski and George Bissias, as well as from discussions with Andrew McGregor, Don Towsley, Berthold Horn, Archan Ray, Justin Domke, Gary Huang, Dan Sheldon, Luc Rey-Bellet, and Markos Katsoulakis. Some of this work grew out of early efforts to improve Anderson’s bound by Benjamin Mears while at the University of Massachusetts.