I'm interested in implementing an algorithm to find the quantile of the beta distribution, and I'm looking at this paper: Journal of the Royal Statistical Society Series C (Applied Statistics). 1973, Vol. 22 Issue 3, p411. 4p, K. L. Majumder and G. P. Bhattacharjee. The paper presents an algorithm, together with a FORTRAN implementation, to find $x$ satisfying
$$ I_x(p,q) = \frac1{B(p,q)}\int_0^xt^{p-1}(1-t)^{q-1}\,\mathrm{d}t \;=\; \alpha$$
where the inputs are
- the quantile $\alpha \in [0,1]$;
- and the parameters of the Beta distribution $p,q>0$.
The algorithm comprises
- Initial approximation, followed by
- Newton-Raphson iteration.
(note that the algorithm has been refined by subsequent papers)
I'm interested in understanding the initial approximation a little better. The paper starts by quoting some results
where $Y_\alpha$ is Hasting's approximation for the upper $100\alpha\%$ point of the standard normal distribution. If $\chi^2_\alpha < 0$, then
$$ x_0 = 1- [(1-\alpha)q\,B(p,q)]^{1/q}\;\;\;(\dagger)$$
Again if $(4p+2q-2)/\chi^2_\alpha$ does not exceed $1$, $x_0$ is obtained from
$$ x_0 = [\alpha p \, B(p,q)]^{1/p}\;\;\;(\ddagger)$$
What I'm interested in is where do $(\dagger)$ and $(\ddagger)$ come from? Thanks for any help.