{'fraction_non_alphanumeric': 0.0670672296200065, 'fraction_numerical': 0.023221825267944137, 'mean_word_length': 3.8649289099526065, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A closed-form approximation for the median of the beta distribution

A simple closed-form approximation for the median of the beta distribution Beta(a, b) is introduced: (a-1/3)/(a+b-2/3) for (a,b) both larger than 1 has a relative error of less than 4%, rapidly decreasing to zero as both shape parameters increase.


Introduction
Consider the the beta distribution Beta(a, b), with the density function, The mean of Beta(a, b) is readily obtained by the formula a/(a + b), but there is no general closed formula for the median. The median function, here denoted by m(a, b), is the function that satisfies, The relationship m(a, b) = 1 − m(b, a) holds. Only for the special cases a = 1 or b = 1 we may obtain an exact formula: m(a, 1) = 2 −1/a and m(1, b) = 1 − 2 −1/b . Moreover, when a = b, the median is exactly 1/2.
There has been much literature about the incomplete beta function and its inverse (see e.g. Dutka (1981) for a review). The focus in literature has been on finding accurate numerical results, but a simple and practical approximation that is easy to compute has not been found.

A new closed-form approximation for the median
Trivial bounds for the median can be derived (Payton et al., 1989), which are a consequence of the more general mode-median-mean inequality (Groeneveld and Meeden, 1977). In the case of the beta distribution with 1 < a < b, the median is bounded by the mode (a − 1)/(a + b − 2) and the mean a/(a + b): For a ≤ 1 the formula for the mode does not hold as there is no mode. If 1 < b < a, the order of the inequality is reversed. Equality holds if and only if a = b; in this case the mean, median, and mode are all equal to 1/2. This inequality shows that if the mean is kept fixed at some p, and one of the shape parameters is increased, say a, then the median is sandwiched between p(a − 1)/(a − 2p) and p, hence the median tends to p.
From the formulas for the mode and mean, it can be conjectured that the median m(a, b) could be approximated by m(a, b; d) = (a − d)/(a + b − 2d) for some d ∈ (0, 1), as this form would satisfy the above inequality while agreeing with the symmetry requirement, that is, m(a, b; d) = 1 − m(b, a; d). . The smaller of the shape parameters is fixed, i.e. for p ≤ 0.5, the median is computed for Beta(a, a(1 − p)/p) and for p > 0.5, the median is computed for Since a Beta(a, b) variate can be expressed as the ratio γ 1 /(γ 1 + γ 2 ) where γ 1 ∼ Gamma(a) and γ 2 ∼ Gamma(b) (both with unit scale), it is useful to have a look at the median of the gamma distribution. Berg and Pedersen (2006) studied the median function of the unit-scale gamma distribution median function, denoted here by M (a), for any shape parameter a > 0, and obtained M (a) = a − 1/3 + o(1), rapidly approaching a − 1/3 as a increases. It can therefore be conjectured that the distribution median may be approximated by, Figure (1) shows that this approximation indeed appears to approach the numerically computed median asymptotically for all distribution means p = a/(a + b) as the (smaller) shape parameter a → ∞. For a ≥ 1, the relative error is less than 4%, and for a ≥ 2 this is already less than 1%. Figure (2) shows the relative error over all possible distribution means p = a/(a + b), as the smallest of the two shape parameters varies from 1 to 4. This illustrates how the relative error tends uniformly to zero over all p as the shape parameters increase. The figure also shows that the formula consistently either underestimates or overestimates the median depending on whether p < 0.5 or p > 0.5.
However, the function m(a, b; d) approximates the median fairly accurately if some other d close to 1/3 (say d = 0.3) is chosen. Figure (3) displays curves of the logarithm of the absolute difference from the numerically computed median for a fixed p = 0.01, as the shape parameter a increases. The absolute difference has been scaled by p before taking the logarithm: due to this scaling, the error stays approximately constant as p decreases so the picture and its scale will not essentially change even if the error is computed for other values of p < 0.5. The figure shows that although some As the smaller of the two shape parameters increases, the tail probability tends rapidly and uniformly to 0.5.
approximations such as d = 0.3 has a lower absolute error for some a, the error of m(a, b; 1/3) tends to be lower in the long run, and moreover performs more consistently by decreasing at the same rate on the logarithmic scale. In practical applications, d = 0.333 should be a sufficiently good approximation of d = 1/3. Another measure of the accuracy is the tail probability Pr(θ ≤ m(a, b; 1/3)) of a Beta(a, b) variate θ: good approximators of the median should yield probabilities close to 1/2. Figure (4) shows that as long as the smallest of the shape parameters is at least 1, the tail probability is bound between 0.4865 and 0.5135. As the shape parameters increase, the probability tends rapidly and uniformly to 0.5.
Finally, let us have a look at a well-known paper that provides further support for the uniqueness of m(a, b; 1/3). Peizer and Pratt (1968) and Pratt (1968) provide approximations for the probability function Pr(θ ≤ x) of a Beta(a, b) variate θ. Although they do not provide a formula for the inverse, it is the probability function at the approximate median. According to Peizer and Pratt (1968), Pr(θ ≤ x) is well approximated by Φ(z(a, b; x)) where Φ is the standard normal probability function, and z is a function of the shape parameters and the quantile x. Consider m = m(a, b; d): z(a, b; m) should be close to zero and at least tend to zero fast as a and b increase. Now assume that p is fixed, a varies and b = a(1 − p)/p. The function z(a, b; m) equals, rewritten with the notation in this paper, where the function f (a, p; d) tends to zero as a increases, being exactly zero only when d = 1/2 or m = 1/2. It is evident that for the fastest convergence rate to zero, one should choose d = 1/3. This is of the order O(a −3/2 ); if d = 1/3, for example if we choose the mean p as the approximation of the median (d = 0), the rate is at most O(a −1/2 ).