On closed-form tight bounds and approximations for the median of a gamma distribution

The median of a gamma distribution, as a function of its shape parameter k, has no known representation in terms of elementary functions. In this work we use numerical simulations and asymptotic analyses to bound the median, finding bounds of the form 2−1/k(A + Bk), including an upper bound that is tight for low k and a lower bound that is tight for high k. These bounds have closed-form expressions for the constant parameters A and B, and are valid over the entire range of k > 0, staying between 48 and 55 percentile. Furthermore, an interpolation between these bounds yields closed-form expressions that more tightly bound the median, with absolute and relative margins to both upper and lower bounds approaching zero at both low and high values of k. These bound results are not supported with analytical proofs, and hence should be regarded as conjectures. Simple approximation expressions between the bounds are also found, including one in closed form that is exact at k = 1 and stays between 49.97 and 50.03 percentile.


Introduction
The known Laurent series for the median of a gamma distribution produces some upper and lower bounds, and very accurate approximations at high enough values of the shape parameter k, but poor results near k = 1 and below.The known approximate and bounding results for k close to zero do not extend to high k.We find new bounds that correspond to percentiles close to 50 across the entire range of shape parameter, and apply them to making even tighter bounds as well as better approximation formulae.An asymptotic approximation and bound at low k, as opposed to the usual focus on high k, provides the key functional form needed for this approach.

Problem formulation
The gamma distribution PDF is θ , but we'll use θ = 1 because both the mean and median simply scale with this parameter.Thus we use this PDF with just the shape parameter k, with k > 0 and x ≥ 0: The mean of this distribution, µ, is well known to be µ = k.The median ν is the value of x at which the CDF equals one-half: This equation has no easy solution, but the median is well known to be a bit below the mean, bounded by [6] k − 1 3 < ν < k and 0 < ν Bounds that are tighter in some part of the shape parameter range can be obtained from the known Laurent series partial sums, or from the low-k asymptote and bounds of Berg and Pedersen [3].
We seek upper and lower bounds that are tighter, especially in the middle part of the k range, than are previously known.Further, we seek simple approximation formulae for the median, leveraging these bounds.

Prior work
A Laurent series for ν with rational coefficients has been discovered, with deep connections to some math by Ramanujan.Choi [7] applied Ramanujan's work to this particular question, providing 4 coefficients (through the k −3 term).Berg and Pedersen [3], based on work by Marsaglia [8], extended this to 10 coefficients.Neither commented on the radius of convergence, which appears to be in the neighborhood of k = 1.So for large enough k and N the series yields excellent approximations, but for k < 1 it is useless.
Thus, where Choi [7] had 144 instead of the correct 2 3 • 23 = 184: Partial sums of the Laurent series are not generally bounds, but the first two (k and k − 1/3) are upper and lower bounds, respectively, and the sums ending with −3 and −5 powers of k are also upper and lower bounds, respectively.
Berg and Pedersen [3] also derived an asymptote for small k, which we call where γ ≈ 0.577216 is the Euler-Mascheroni constant.This asymptote is a lower bound, as we will show.Their factor 2 −1/k is the key to good approximations and bounds, and we divide by it to reduce the dynamic range in some of our later plots.For the range 0.01 < k < 100, this reduces the dynamic range we need to work with by nearly 30 orders of magnitude-or 300 orders of magnitude for k down to 0.001-but we still need log-log plots to show the bounds and approximations, and their errors, across these ranges.Berg and Pedersen [3] also provide an upper bound ν < e −1/3k k that is just above k − 1  3 , and a lower bound ν > 2 −1/k k.These previous known bounds are illustrated in Figure 1, where upper and lower bounds and their errors (or "margins") are distinguished by different line styles.
Others have shown good bounds and approximations where k is an integer, that is, for the Erlang distribution [1,10,5].These results do not extend to the low-k region.
The approach of approximating functions by interpolating between upper and lower bounds has been discussed by Barry [2], who used minimax optimization to find numeric parameters in an interpolation-between-bounds approximation to the exponential integral.We are not aware of an interpolation approach being used to find improved closed-form bounds.

Tight upper and lower bounds
The Berg lower bounds 2 −1/k e −γ and 2 −1/k k are tight within their families 2 −1/k A and 2 −1/k Bk, but not tight within the wider two-parameter family 2 −1/k (A + Bk).We show that the sum of these two is the uniquely tight upper bound in that family, and that there is a range of tight lower bounds in the family.
To improve on the lower bounds, and to motivate the family that we consider further, first solve for ν in a simple approximation to the distribution's integral, using e −x < 1 for x > 0: which is a tight lower bound, and is a good approximation for k < 0.1, but not so great at high k-and not what we consider a closed form, due to the gamma function.This expression resembles the "quantile mechanics" boundary condition for the gamma distribution from Steinbrecher and Shaw [9], and converges Figure 1: Previously published bounds (red; lower bounds solid, upper bounds dashed) for the median of a gamma distribution (black dotted), are good at high k or low k, but not both.Their margins (errors) are shown in blue (corresponding line styles).At the left, at k = 0.01, the median is near 10 −30 , and at k = 0.001 near 10 −300 .
with the Berg asymptote at low k, but we can improve Berg's result by utilizing another term.A symbolic calculus system finds for us the next Taylor series terms about k = 0 for the power of the gamma function: So we have this improved approximation, which has much less relative error than Berg's at low k, and has a high-k behavior nearly proportional to k (but is no longer a lower bound because we made it larger by ignoring a next negative term): Inspired by this asymptotic approximation, we consider members of this family of functions, with coefficients A and B, and analyze which ones are bounds: A graphical characterization of this family is most informative.Given some values of k and corresponding numerical ν(k), we can find the lines A + Bk = ν(k)2 1/k in A-B space, and plot them-see Figure 2. Regions full of lines are not bounds, and regions without lines are where bounds are found (including some of Berg's bounds); we're interested in the boundaries between these regions, where tight bounds are to be found.
The improved asymptote ν 1 is great at low k, but is neither an upper nor a lower bound, as shown in 3. We can modify it to approach k − 1 3 at high k by a few adjustments, via this asymptotic approximation that we get from a symbolic calculus system: Thus we find this approximation for high k, which is a lower bound as illustrated in Figure 2: A compromise approximation for low k mixes these two, differing from the high-k approximation in only the A coefficient, leaving a result consistent to the same order as Berg's at low k, and forming an upper bound as illustrated in 2: This mixed approximation has absolute and relative errors approaching zero at low k, and relative error approaching zero at high k; but the absolute error remains high, near log 2 − 1 3 − e γ ≈ 0.20, at high k.These approximations and their errors are illustrated in Figure 4.  Figure 3: Zooming in to ν 1 and ν L0 , note that the point with B = e −γ π 2 /12, which we got from the Taylor series of the power of the gamma function, is actually inside the shaded area, so does not represent a bound; but a point at slightly lower B = 0.45965 is on the edge, so represents a lower bound.These points give zero error at approximately k = 0.1003 and k = 0.0708, respectively (see the next figure).We do not have analytic formulations for these numeric and graphical observations.The lower bounds ν L0 , ν L1 , and ν L∞ (solid), and upper bound ν U (dashed) are shown in red over the ideal median (black heavy dots), with their absolute errors in blue, all premultiplied by 2 1/k to reduce the required plot range.The approximation ν 1 , which is not a bound, is also shown; note that its error curve changes from solid to dashed at the cusp, while the errors for ν L0 and ν L1 have log(0) cusps where the error grazes zero but does not change sign.
The k parameters at these cusps correspond to the sloped lines indicated in the previous figures.
To support the graphical/numerical observation that ν L∞ and ν U are lower and upper bounds, respectively, of the true median (ν L∞ < ν < ν U ), we examine their asymptotic behaviors in more detail.At low k, it is easy to see, using log 2 − 1 3 ≈ 0.359814 < e −γ ≈ 0.561459, and e −γ π 2 12 ≈ 0.461781 < 1, that these differences are positive, for k → 0: At high k, a symbolic calculus system gives us for ν L∞ : which we can use to construct comparisons to the Laurent series terms for ν [7,3].Again we find positive differences, with log 2 3 − log 2 2 2 ≈ −0.009177 < 8 405 , and log 2 − e −γ ≈ 0.131688 < 1  3 , for k → +∞: In addition to these high-k and low-k asymptotic results, we can show the inequalities also hold at k = 1 where 1 2 log 2 − 1 3 + 1 < log 2 < 1 2 (e −γ + 1), but otherwise we're relying on the graphical and numerical results.Since there is considerable margin in the asymptotes, and the median is well behaved (unique, monotonic, and smooth, with positive second derivative for all k > 0 [3,4]), this seems reliable enough in concluding that these are bounds (but mathematicians are invited to interpret these as conjectures to be proved).The positivity constraints would not all hold, since higher-order terms would not cancel, if A L∞ or B L∞ were any higher, or if A U or B U were any lower.In that sense, these upper and lower bounds are proved tight.
For the lower bound ν L1 that is tight at k = 1, both the value and the slope need to match the true median.The value is the median of the exponential distribution, ν(1) = log 2. The slope v ′ (k), which is somewhat more troublesome to work out, but is tractable at the special point k = 1, is: which is a mathematical expression with a definite value, but is not a closed form due to the exponential integral, so still requires a numerical approach to evaluate it.Therefore, we describe this bound with approximate numerical parameters instead of closed-form analytic expressions.
The lower bound ν L0 is in worse shape, as we have to search for the k value that gives the lowest B value with A = e −γ .So its B parameter has no concise mathematical expression, but can be computed to high precision: B ≈ 0.4596507.
The bounds and asymptotic approximations discussed here are summarized in Table 1.In subsequent sections we focus primarily on the new upper and lower bounds with closed-form coefficients, ν U and ν L∞ , as a basis for even tighter closed-form bounds.
Figure 5 shows the percentile values achieved by these bounds and approximations, compared to the ideal 50% that defines the median-ν L∞ always comes in between 48% and 50% and ν U always between 50% and 55% (percentiles are calculated in Matlab as 100*gammainc(x, k), using the normalized lower incomplete gamma function that is the CDF for our PDF).
The coefficient of k −1 for ν L∞ is negative, so the lower bound ν L∞ slightly violates the k − 1 3 lower bound in spite of having asymptotically zero absolute error.That is, for k > 3.021, it's a looser lower bound and a worse approximation than k − 1/3, even though it is the tightest lower bound of the form we're considering.On the other hand, ν U is a much tighter upper bound than k is.

Formulae for tighter bounds
Letting A and B be functions of k, rather than constants, allows tighter bounding expressions (and potentially exact expressions) for ν(k), but not enough structure.Allowing only one of them to vary, and tying the other to values used in the tight bounds above, allows a more constrained space of bounds.(dotted), upper and lower bounds from Berg and Pedersen [3] (dash-dot), and a pair of closer bounds formed by interpolation between ν U and ν L∞ using a one-parameter rational function (dashed).The bounds 2 1/k k, ν U , ν L∞ , and the interpolated bounds converge on 50th percentile at both low and high k, while the other six do not.Both the upper and lower interpolated bounds are close to ν U at low k and close to ν L∞ at high k; tighter such interpolated bounds, developed in a later section, would crowd the center of the graph.

Ideal and bounding interpolators and margins
Figure 7: The ideal interpolator g(k) (heavy dotted sigmoid) is compared with upper and lower bounds g(k); their margins g(k) − g(k) are also plotted, magnified and displaced, with the same curve styles.The curves with largest absolute margins (dashed), which correspond to the interpolated bounds shown in Figure 5, are for the first-order rational-function interpolator k k+b0 , while the curves with smaller margins (solid) are for arctan interpolators 2 π tan −1 k b .In each case, the one parameter (b 0 or b) is chosen to give a tight bound (analytically in closed form in three of the four cases).Lower bounds of g(k) make upper bounds of ν(k), and vice versa.
At high k: And at k = 1: How such approximation goals relate to bounds is not immediately clear.We construct some examples.Figure 7 shows the ideal interpolator, computed numerically, and compares it to bounding interpolators of the forms g(k) = k k+b0 and g(k) = 2 π tan −1 k b , with parameters b 0 and b chosen to yield tight upper and lower bounds, as discussed in the next two sections.The effects of these interpolator bounds on the median bounds is shown in Figure 8.

Rational-function interpolators
Consider rational functions as interpolators, of the form For N = 1, the only parameter is b 0 , so we have a one-parameter family: For N = 2 we have three parameters: and so forth.
We can easily constrain the coefficients to match the properties of the ideal interpolator.At low k: And at k = 1: For N = 1, the low asymptote is tightly approached with b 0 = 1 P0 , yielding an upper bound for the median (lower bound for the interpolator g ideal ).Or the high asymptote is tightly approached with b 0 = P ∞ , yielding a lower bound for the median (upper bound for g(k)).These bounds are illustrated in Figure 7.For b 0 between these values, the resulting interpolated function is not a bound, but is exact at one value of k, for example at k = 1 with b 0 = 1 P1 − 1.For N = 2, there are enough parameters to use any or all of the three constraints, but it's not immediately clear which sets of constraints can lead to bounds.Certainly using all three constraints does not lead to a bound, but to an interesting approximation.As we did with the A versus B space, we can investigate the locus of parameter solutions for each k, and examine the edges of these locus-filled areas for tight bounds.Reducing the space to 2D by constraining one asymptote or the other allows a graphical approach, but the results are not better than a one-parameter arctan interpolator.
For N ≥ 3, excellent approximations and bounds are possible, but with so many parameters are not very interesting.Figure 9: Relative errors of arctan-interpolated approximations (dash-dot curves) between the upper (dashed) and lower (solid) bounds.These are among the possibly interesting approximations suggested by the tight-bounds approach.The non-bounding approximations, optimized for different criteria, all have maximum relative errors below 1%.

Arctan interpolators
Getting a tighter bound or better approximation from the rational-function family requires at least a handful of parameters.An alternative approach is to find a one-parameter shape that fits better.We have found that the arctan shape with parameter b (like b 0 in the one-parameter rational-function interpolator, corresponding to the k value at the midpoint, g(b) = 0.5) does a good job: As Figure 7 shows, the shapes of the arctan interpolators are imperfect, but are much better than the first-order rational function, fitting better in some regions than in others.Note that both the N = 1 rational function and the arctan interpolators are symmetric about their centers (with log k as the independent variable); they are logistic sigmoid and Gudermannian shapes, respectively.The ideal that they are to bound or approximate, however, is not quite symmetric in log k.So a family of not-quite-symmetric interpolators can perhaps do better.
Several special b parameter values for the arctan interpolator can be derived to match each of the ideal properties mentioned above.We can constrain the π P −1 0 , which yields a lower bound to g.
To find an upper bound to g (lower bound to ν), we decrease b until the margin is nonnegative for all k, which is at about b = 0.205282.Finding an analytic formulation for that tight bound is a challenge left to others.
Several interpolated bounds and approximations are summarized in Table 2. Most are closed-form analytic expressions; see Figure 9 for the relative errors of various arctan interpolator versions.

Conclusions
Tight upper and lower bounds to the median of the gamma distribution are introduced.The simplest new lower bound is never below the 48th percentile, and the simplest new upper bound, of the same form 2 −1/k (A + Bk), is never above the 55th percentile, over the entire range of k > 0. Using arctan and rational-function interpolators between these bounds, two one-parameter fami-lies of closed-form bounds and approximations to the median of a gamma distribution are proposed.
The one-parameter rational-function family has simple closed-form formulae for tightest upper and lower bounds, staying below 50.85 and above 49.69 percentile, respectively; higher-order rational functions can provide tighter bounds or better approximations.
The one-parameter arctan family of interpolators is a better fit to the ideal interpolator, and includes a version that is most accurate in the low-k tail and provides a closed-from tight upper bound, staying below 50.18 percentile.With different b parameter, several approximations in the family, including the closedform version that is exact at k = 1, stay between 49.97 and 50.03 percentile.We have not found an analytic formula for the parameter that gives the tight lower bound, which stays above 49.96percentile, but have shown where to find it, graphically or numerically.
The approach of interpolating between tight bounds opens the way to finding tighter bounds and more accurate approximations, and to finding more such families of bounds and approximations via other interpolator forms.

3 Figure 2 :
Figure 2:The A-B parameter space is shaded with dash-dot lines where ν(k)2 1/k = A + Bk, for a set of very small to very large k values in geometric progression (using numerically computed ν(k) values).Key values of A and B are indicated.Points outside (or on the edge) of the shaded region represent bounds, while points inside the shaded region represent functions that cross the median function.There is an obvious uniquely tight upper bound ν U , and a curved locus of tight lower bounds from ν L0 , which is tightest near k = 0, to ν L∞ , which is tightest for high k.One point (pentagram) on the curved locus represents a lower bound ν L1 that is tight at k = 1, for which A + B = 2 log 2 (which is the equation of the dashed line).The point ν 1 represents a good asymptotic approximation close to ν L0 , but not a bound; see the next figure.The dotted line from ν U to ν L∞ at B = 1 intersects the lines for all k in monotonic order, with A decreasing while k increases.

Figure 4 :
Figure4: The lower bounds ν L0 , ν L1 , and ν L∞ (solid), and upper bound ν U (dashed) are shown in red over the ideal median (black heavy dots), with their absolute errors in blue, all premultiplied by 2 1/k to reduce the required plot range.The approximation ν 1 , which is not a bound, is also shown; note that its error curve changes from solid to dashed at the cusp, while the errors for ν L0 and ν L1 have log(0) cusps where the error grazes zero but does not change sign.The k parameters at these cusps correspond to the sloped lines indicated in the previous figures.

Figure 5 :
Figure5: The percentiles achieved by four new median bounds of the form 2 −1/k (A + Bk) (solid curves) are plotted, along with the linear bounds k and k − 1 3 (dotted), upper and lower bounds from Berg and Pedersen[3] (dash-dot), and a pair of closer bounds formed by interpolation between ν U and ν L∞ using a one-parameter rational function (dashed).The bounds 2 1/k k, ν U , ν L∞ , and the interpolated bounds converge on 50th percentile at both low and high k, while the other six do not.Both the upper and lower interpolated bounds are close to ν U at low k and close to ν L∞ at high k; tighter such interpolated bounds, developed in a later section, would crowd the center of the graph.

Figure 6 :
Figure 6: Functions A(k) and B(k), either of which can solve ν = 2 −1/k (A+Bk), with the other constant at the limiting values indicated by the circles, the parameters of ν U ; modeling these curves can lead to better bounds or approximations.

Figure 8 :
Figure 8: The absolute margins of the interpolated bounds are smaller than those of the bounds they started from.Compare Figure 4.The generally smallest margins are for the arctan interpolators, and the intermediate for the N = 1 rational-function interpolators.

Table 2 :
Comparison of several one-parameter interpolated bounds and approximations g(k)ν L∞ + (1 − g(k))ν U , some of which have closed-form parameters.Bounds are indicated by U or L. through the known value ν = log 2 at k = 1, with b = cot π 2 P 1 , but that does not give a bound.Or we match the true median at high k to within O(k −2 ) with b = π 2 P ∞ .That also does not give a bound.Or we can match at low k to within O(2 −1/k k 2 ), like ν1 does, with b = 2

Table 1 :
Comparison of several median bounds and asymptotes.