# Approximating the median of a $\Gamma(\alpha,1)$ distribution with $0<\alpha<1$

Is there a good approximation (or useful bounds) for the median $$\nu_\alpha$$ of a $$\Gamma(\alpha,1)$$ distribution with $$0<\alpha<1$$?

I have only been able to find things like Berg & Pedersen (2006, Methods and Applications of Analysis), who generalize earlier work of Choi (1994, Proceedings of the American Mathematical Society) and give a good asymptotic expression for $$\nu_\alpha$$ as $$\alpha\to\infty$$, which works pretty well for $$\alpha\geq 1$$ - but nothing whatsoever for $$0<\alpha<1$$.

(I am interested in this because would like to update my earlier Q&A about minimizing MAPEs and other forecasting errors for gamma distributed future outcomes. So far, it only works if $$\alpha\geq 2$$, and the case $$1<\alpha<2$$, which leads to getting $$\nu_\alpha$$ for $$0<\alpha<1$$ as above, is missing.)

• Since most stats packages implement an inverse-cdf (quantile function) for the gamma, even when $\alpha<1$, presumably that could be used. Or were you after a closed-form expression? – Glen_b -Reinstate Monica Nov 27 '19 at 11:28
• @Glen_b-ReinstateMonica: yes, and that works well numerically, and yes, I was hoping for a closed form similar to the one in section 6 in Berg & Pedersen (2006). – Stephan Kolassa Nov 27 '19 at 11:34
• This is probably extremely naive, but couldn'' you calculate the medians numerically with qgamma for $\alpha$ between $0$ and $1$ and then fit some sort of polynomial model or something else to these? – COOLSerdash Nov 27 '19 at 19:20
• @COOLSerdash: that would be a possibility, yes. The convexity of the median as a function of $\alpha$ (Berg & Pedersen, 2008) would likely help. What I am looking for is error guarantees, so perhaps someone has already done something along these lines using interpolation. – Stephan Kolassa Nov 27 '19 at 20:06
• @whuber: That sounds interesting. Do you have pointers to any of those results? Do the simple approximations come with bounds? The Berg & Pedersen expansion seems to work well even non-asymptotically, it's the missing bounds that bothers me. – Stephan Kolassa Feb 25 at 7:57

Let $$\mu_\alpha$$ be the median of a $$\Gamma(\alpha)$$ distribution. This means the area under the density $$f_{\Gamma(\alpha)}(x) = \frac{x^{\alpha-1}}{\Gamma(\alpha)} \,e^{-x}$$ between $$x=0$$ and $$x=\mu_\alpha$$ equals $$1/2.$$ A graph of $$f_{\Gamma(\alpha)}$$ is sketched here in black (for $$\alpha=0.3$$), understanding the graph extends infinitely upwards as $$x$$ approaches $$0$$ and flattens down to $$0$$ as $$x$$ grows large: The median $$\mu_\alpha$$ separates the left half of the area (darkened) from the right half. The dotted red curve is an upper bound for the Gamma density, enabling a lower bound for the median to be found by using the area under the red curve instead of the gray area.

$$1 - x \le e^{-x} \le 1$$ for $$x \ge 0$$ implies

$$\frac{x^{\alpha-1}}{\Gamma(\alpha)} \,(1-x) \le f_{\Gamma(\alpha)}(x) \le \frac{x^{\alpha-1}}{\Gamma(\alpha)}$$

which in turn gives $$\int_0^{\mu_\alpha}\frac{x^{\alpha-1}}{\Gamma(\alpha)} \,(1-x)\,\mathrm{d}x \le \int_0^{\mu_\alpha}f_{\Gamma(\alpha)}(x)\,\mathrm{d}x =\frac{1}{2} \le \int_0^{\mu_\alpha}\frac{x^{\alpha-1}}{\Gamma(\alpha)}\,\mathrm{d}x.$$

Evaluating the integrals produces

$$\frac{\mu_\alpha^\alpha}{\alpha\Gamma(\alpha)}\left(1 - \frac{\alpha}{\alpha+1}\mu_\alpha\right) \le \frac{1}{2} \le \frac{\mu_\alpha^\alpha}{\alpha\Gamma(\alpha)}.$$

This can be simplified a little by solving for $$\mu_\alpha$$ and recalling $$z\Gamma(z) = \Gamma(z+1)$$ for any $$z:$$

$$\mu_\alpha \ge \left(\frac{\Gamma(\alpha+1)}{2}\right)^{1/\alpha} \ge \mu_\alpha\left(1 - \frac{\alpha}{\alpha+1}\,\mu_\alpha\right)^{1/\alpha}.$$

One way to exploit the second inequality is to find a bound for $$\mu_\alpha/\alpha$$ when $$0 \lt \alpha \le 1$$. This expression is increasing: we have

$$\frac{\mu_\alpha}{\alpha} = e^{-\varphi(\alpha)}$$

for $$\varphi$$ as defined in equation (3) in Berg & Pedersen (2006), and Proposition 3.6 in the same paper shows that $$\varphi$$ is decreasing.

Therefore the value of $$\frac{\mu_\alpha}{\alpha}$$ at $$\alpha=1$$, equal to $$\log(2)$$, provides an upper bound

$$\mu_\alpha \le \alpha \log(2).$$

Crude as this is, it enables us to eliminate $$\mu_\alpha$$ from the denominator of the right hand side by replacing it with its upper bound, yielding the bounds

$$\left(\frac{\Gamma(\alpha+1)}{2}\right)^{1/\alpha} \le \mu_\alpha \le \left(\frac{\alpha + 1}{\alpha + 1 - \log(2) \alpha^2}\right)^{1/\alpha}\,\left(\frac{\Gamma(\alpha+1)}{2}\right)^{1/\alpha}.$$

Dividing these bounds by $$\mu_\alpha$$ gives the relative error plotted here: The lower bound clearly is an accurate approximation for $$\alpha \lt 0.2,$$ yielding at least two significant decimal digits, while the upper bound gives at least one significant digit throughout.

Having obtained definite bounds (useful for analysis), we may adjust them to produce even more accurate estimates. For instance, $$0.9075$$ times the upper bound approximates $$\mu_\alpha$$ to within one percent relative accuracy when $$0.16\le\alpha\le 1$$ and $$1.0035$$ times the lower bound achieves $$0.4\%$$ relative accuracy for $$\alpha \lt 0.16.$$