The beta distribution is ubiquitous in statistics due to its flexibility. Among many applications, it is used in analyses of uniform order statistics [17], problems in Euclidean geometry [10], general theory of stochastic processes [24], and applied statistical inference; the last category of applications includes time estimation in project management [5], hypothesis testing [34], A/B testing in business [28], modelling in life sciences [31] and others [19, 20].

Unfortunately, the importance and simplicity do not go in pairs. The distribution of X ∼ Beta ( α , β ) with parameters α , β > 00$]]> is given by (1) P { X ⩽ ϵ } = ∫ 0 ϵ x α − 1 ( 1 − x ) β − 1 B ( α , β ) d x , 0 ⩽ ϵ ⩽ 1 , where B ( α , β ) ≜ ∫ 0 1 x α − 1 ( 1 − x ) β − 1 d x is the normalizing constant; the integral (1), also known as the incomplete beta function [7], is intractable. Thus, there is a strong demand for closed-form approximations of tail probabilities, for example, in the context of adaptive Bayesian inference [8], Bayesian nonparametric statistics [4], properties of random matrices [10, 26], and (obviously) large deviation theory [33].

The main goal of this paper is to give accurate closed-form bounds for the beta distribution in the form of sub-gamma type exponential concentration inequalities [3], that is, (2) P { X − E [ X ] < − ϵ } , P { X − E [ X ] > ϵ } ⩽ exp − ϵ 2 2 v + 2 c ϵ \epsilon \}\leqslant \exp \left(-\frac{{\epsilon ^{2}}}{2v+2c\epsilon }\right)\]]]> for any ϵ ⩾ 0 and the two parameters: the variance proxy v and the scale c, both depending on α, β. Such concentration bounds, pioneered by Bernstein [2] and popularized by Hoeffding [14], are capable of modelling both the sub-Gaussian and sub-exponential behaviors. Due to this flexibility, bounds of this sort are the working horse of approximation arguments used in modern statistics [18, 3, 22, 30].

The contributions of this work are as follows: •

Optimal Bernstein-type concentration inequality. We establish an exponential concenration bound (2) with v that matches the variance (which is optimal) and, respectively, the optimal value of c which depends on the ratio of the third and second moment. This bound is shown optimal in the regime of small deviations.

When judging the bounds in the form of (2), it is important to insist on the optimality of the sub-Gaussian behavior. For small deviations ϵ the bound (2) becomes approximately Gaussian with variance v, thus we ideally want v 2 = Var [ X ] and exponential bounds of this type are considered optimal in the literature [1]. On the other hand, bounds with v 2 > Var [ X ]\mathbf{Var}[X]$]]> essentially overshoot the variance, leading to unnecessary wide tails and increasing uncertainty in statistical inference.

Bearing this in mind, we review prior suboptimal bounds in the form of (2): •

Folklore methods give some crude bounds, for example, one can express a beta random variable in terms of gamma distributions and utilize their concentration properties; such techniques do not give optimal exponents.

The remainder of the paper is organized as follows: Section 2 presents the results, Section 3 provides the technical background, Section 4 gives proofs and Section 5 concludes the work. The implementation and numerical evaluation are shared in the notebook [27].

The central result of this work is a Bernstein-type tail bound with numerically optimal parameters. The precise statement is given in the following theorem.

The result below shows that the parameters c and v are best possible, in the regime of small deviations ϵ, for the exponential moment method.

Comparing this result with the tail from Theorem 1 as ϵ → 0 we obtain: Corollary 1

The values of constants c and v in the Bernstein-type inequality in Theorem 1 are optimal.

Our optimal Bernstein inequality is proved using the novel recursion for central moments, presented below in Theorem 3 as the contribution of independent interest. Being of order 2 it is not only efficient to evaluate numerically, but also easy to manipulate algebraically when working with closed-form formulas.

The implementation of the recursion, performed in Python’s symbolic algebra package Sympy [23], is presented in Listing 1. Scalability is ensured by memoization, which caches intermediate results to avoid repeated calls.

To demonstrate the algorithm in action, we list some first central moments of the beta distribution in Table 1.

Note that the algorithm addresses the lack of simple derivations for such formulas, as well as the lack of formulas beyond skewness and kurtosis.

Finally, we note that the recurrence implies, by induction, the following important property:

The experiment summarized in Figure 1 (the code shared in [27]) demonstrates the advantage of sub-gamma bounds obtained in Theorem 1 over the optimal sub-Gaussian bounds from prior work [21].

For skewed beta distributions, the bound from this work is more accurate than the sub-Gaussian approximation. The chosen range of parameters covers the cases where the expectation E [ X ] is a relatively small number like 0.02 or 0.2, a typical range for many practical Beta models, particularly A/B testing. Note that there is still room for improvement in the regime of larger deviations ϵ, where the bounds could potentially benefit from refining the numeric value of c. The experiment details are shared in the Colab Notebook [27].

The Gaussian hypergeometric function is defined as follows [11]: (8) 2 F 1 ( a , b ; c ; , z ) = ∑ k = 0 + ∞ ( a ) k ( b ) k ( c ) k · z k k ! , where we use the Pochhammer symbol defined as (9) ( x ) k = 1 , k = 0 , x ( x + 1 ) ⋯ ( x + k − 1 ) , k > 0 . 0.\end{array}\right.\]]]> We call two functions F of this form contiguous when their parameters differ by integers. Gauss considered 2 F 1 ( a ′ , b ′ ; c ′ ; , z ) where a ′ = a ± 1, b ′ = b ± 1, c ′ = c ± 1 and proved that between F and any two of these functions, there exists a linear relationship with coefficients linear in z. It follows [29, 15] that F and any two of its contiguous series are linearly dependent, with the coefficients being rational in parameters and z. For our purpose, we need to express F by the series with increased second argument. The explicit formula comes from [32]: Lemma 1

The following recurrence holds for the Gaussian hypergeometric function: (10) 2 F 1 ( a , b ; c ; , z ) = 2 b − c + 2 + ( a − b − 1 ) z b − c + 1 2 F 1 ( a , b + 1 ; c ; , z ) + ( b + 1 ) ( z − 1 ) b − c + 1 2 F 1 ( a , b + 2 ; c ; , z ) .

We use the machinery of hypergeometric functions to establish certain properties of the beta distribution. The first result expresses the central moment in terms of the Gaussian hypergeometric function.

The proof appears in Section 4.1.

Below, we review the canonical method of obtaining concentration inequalities from the moment generating function, following the discussion in [3].

Recall that for a given random variable Z we define its moment generating function (MGF) and, respectively, cumulant generating function as (12) ϕ Z ( t ) ≜ E exp ( t Z ) , ψ Z ( t ) ≜ log E exp ( t Z ) . In this notation, Markov’s exponential inequality becomes (13) P Z ⩾ ϵ ⩽ exp ( − t ϵ ) E exp ( t Z ) = exp ( − t ϵ + ψ Z ( t ) ) , valid for any ϵ and any nonnegative t, and nontrivial when ϕ Z ( t ) is finite. The optimal value of t is determined by the so-called Cramér transform: (14) ψ ∗ ( ϵ ) = sup t ϵ − ψ Z ( t ) : t ⩾ 0 , and leads to the Chernoff inequality (15) P Z ⩾ ϵ ⩽ exp ( − ψ Z ∗ ( ϵ ) ) .

It is well known that log ( 1 + x ) = x + x 2 2 + ⋯ for nonnegative x. Below, we present a useful refinement of this expansion:

The bound is illustrated in Figure 2 below, and we see that it matches up to the term O ( x 4 ). The proof appears in Section 4.3.

We know that the raw higher-order moments of X are given by [16] (17) E [ X d ] = ( α ) d ( α + β ) d . Combining this with the binomial theorem, we obtain (18) E [ ( X − E [ X ] ) d ] = ∑ k = 0 d ( − 1 ) d − k d k E [ X k ] α α + β d − k = ∑ k = 0 d ( − 1 ) d − k d k ( α ) k ( α + β ) k α α + β d − k . Finally, by d k = ( − 1 ) k ( − d ) k k ! and the definition of the hypergeometric function, (19) E [ ( X − E [ X ] ) d ] = ( − 1 ) d ∑ k = 0 d ( − d ) k ( α ) k k ! ( α + β ) k α α + β d − k = − α α + β d ∑ k = 0 d ( − d ) k ( α ) k k ! ( α + β ) k α + β α k = − α α + β d 2 F 1 α , − d ; α + β ; α + β α , which finishes the proof.

The goal is to prove the recursion formula (10). One way to accomplish this is to reuse the summation formulas developed in the proof of Lemma 2. Below, we give an argument that uses the properties of hypergeometric functions to better highlight their connection to the beta distribution.

To simplify the notation, we define a = α, b = − d, c = α + β, and z = α + β α . Define μ d ≜ E [ ( X − E [ X ] ) d ]. Then using Lemma 2 and Lemma 1 we obtain (20) ( − z ) d · μ d = 2 F 1 ( a , b ; c ; z ) = 2 b − c + 2 + ( a − b − 1 ) z b − c + 1 2 F 1 ( a , b + 1 ; c ; , z ) + ( b + 1 ) ( z − 1 ) b − c + 1 2 F 1 ( a , b + 2 ; c ; , z ) . In terms of α, β, d we obtain (21) 2 b − c + 2 + ( a − b − 1 ) z b − c + 1 = ( d − 1 ) · α − β α ( α + β + d − 1 ) , ( b + 1 ) ( z − 1 ) b − c + 1 = ( d − 1 ) · β α ( α + β + d − 1 ) . The computations are done in SymPy package [23], as shown in Listing 2.

Since we have (22) 2 F 1 ( a , b + 1 ; c ; , z ) = 2 F 1 ( a , − d + 1 ; c ; , z ) = μ d − 1 · ( − z ) d − 1 , 2 F 1 ( a , b + 2 ; c ; , z ) = 2 F 1 ( a , − d + 2 ; c ; , z ) = μ d − 2 · ( − z ) d − 2 , it follows that (23) z 2 μ d = − z ( d − 1 ) ( α − β ) α ( α + β + d − 1 ) · μ d − 1 + ( d − 1 ) β α ( α + β + d − 1 ) · μ d − 2 . Recalling that z = α + β α we finally obtain (24) μ d = − ( d − 1 ) ( α − β ) ( α + β ) ( α + β + d − 1 ) · μ d − 1 + ( d − 1 ) α β ( α + β ) 2 ( α + β + d − 1 ) · μ d − 2 , which finishes the proof.

Consider the function f ( x ) = log ( 1 + x ) − ( x − x 2 2 ( 1 + 2 x 3 ) ) for nonnegative x. Using Sympy we find that (25) f ′ ( x ) = − x 3 x + 1 2 x + 3 2 , as demonstrated in Listing 3. Thus, f is nondecreasing; since f ( 0 ) = 0 we obtain f ( x ) ⩽ 0 as claimed.

If X ∼ Beta ( α , β ), then for X ′ = 1 − X we obtain X − E [ X ] = E [ X ′ ] − X ′ and thus the upper tail of X equals the lower tail of X ′ , and the lower tail of X equals the upper tail of X ′ . Furthermore, from (1) it follows that X ′ ∼ Beta ( β , α ). Thus, by exchanging the roles of α and β accordingly, it suffices to prove the theorem for the upper tail.

To facilitate calculations, we introduce the centred version of X: (26) Z ≜ X − E [ X ] . With the help of normalized moments (27) m d ≜ E [ Z d ] d ! , we expand the moment generating function of Z into the series (28) ϕ ( t ) ≜ E exp ( t Z ) = 1 + ∑ d = 2 + ∞ m d t d , t ∈ R (which converges everywhere as | Z | ⩽ 1), and expand its derivative as (29) ϕ ′ ( t ) = ∑ d = 2 + ∞ d m d t d − 1 , t ∈ R (everywhere, accordingly).

The proof strategy is as follows: we start by expressing the cumulant generating function as an integral involving the moment generating function and its derivative; then we study series expansions of the functions involved, using the moment recursion; this leads to a tractable upper-bound on the integral; finally, the cumulant bound is optimized by the Cramér–Chernoff method. To achieve the promised optimal constants, it is critical to keep all intermediate estimates sharp up to the order of 4. The proof flow is illustrated on Figure 3.

Observe that the cumulant generating function ψ can be expressed in terms of the moment generating function ϕ as follows.

We now present a recursion for coefficients of the series expansions.

The relation between ϕ ′ and ϕ is given in the following auxiliary result. Claim 3.

For nonnegative t and c = β − α ( α + β ) ( α + β + 2 ) it holds that (33) ϕ ′ ( t ) ϕ ( t ) ⩽ v t 1 − c t , α ⩽ β and t < 1 c , ϕ ( t ) ⩽ exp v t 2 2 , α > β . \beta .\end{array}\]]]>

Finally, in the next claim, we present the Chernoff inequality.

Theorem 1 follows now by (15) and replacing c = β − α ( α + β ) ( α + β + 2 ) by 2 c.

Fix ϵ ⩾ 0. Denote Z = X − E [ X ], and let ϕ ( t ) = E [ exp ( t Z ) ], ψ ( t ) = log ϕ ( t ) be, respectively, the moment and cumulant generating functions of Z. By Jensen’s inequality ϕ ( t ) ⩾ exp ( t E [ Z ] ) = 1, and so ψ ( t ) ⩾ 0 and t ϵ − ψ ( t ) ⩽ 0 for t ⩽ 0. Since ψ ( 0 ) = 0 we conclude that (47) ψ ∗ ( t ) = sup t ϵ − ψ ( t ) : t ⩾ 0 = sup t ϵ − ψ ( t ) : t ∈ R , which shows that the Cramér–Chernoff exponent equals the global maximum of the function t → t ϵ − ψ ( t ) (the Legendre–Fenchel transform), achieved for t ⩾ 0.

It is known that the cumulant generating function is convex, and strictly convex for nonconstant random variables [3]. Thus, ψ ( t ) is strictly convex. As a consequence, we see that the function t → t ϵ − ψ ( t ) is strictly concave. Thus, (47) is maximized at the value t which is a unique solution to (48) ψ ′ ( t ) − ϵ = 0 .

The last equation defines t = t ( ϵ ) implicitly and seems to be not solvable with elementary functions. Nevertheless, the Lagrange inversion theorem ensures that t is analytic as a function of ϵ and gives initial terms of the series expansion. More precisely, if y = f ( x ) where f is analytic at 0 and f ( 0 ) = 0, f ′ ( 0 ) ≠ 0 then x expands into a power series in y around y = 0, with the coefficient of the term y k equal to 1 k times the coefficient of x k − 1 in the series expansion of x / f ( x ) k when k > 00$]]> and 0 when k = 0; for a reference, see for example, [9]. Applying this to y = ϵ, x = t and f = ψ ′ we obtain (49) t = ϵ v − c ϵ 2 2 v 2 + O ( ϵ 3 ) , where we do computations in Sympy as shown in Listing 6. Note that we use only up to three terms of the expansion of f, since we are interested in the first two terms for x.

Therefore, we obtain the following exponent in the Chernoff inequality: (50) ψ ∗ ( ϵ ) = t ϵ − ψ ( t ) = ϵ 2 2 v − c ϵ 3 6 v 2 + O ϵ 4 , which finishes the proof of the formula in Theorem 2. This expansion matches the expansion of the exponent in Theorem 1 up to O ( ϵ 4 ), which proves the optimality of the numeric constants.