[Paper Review] Efficient algorithms for the inversion of the cumulative central beta distribution

Amparo Gil, Javier Segura+1 more|arXiv (Cornell University)|May 11, 2016

Fractional Differential Equations Solutions11 references3 citations

TL;DR

This paper presents highly accurate and efficient algorithms for inverting the cumulative central beta distribution using the Schwarzian-Newton method (SNM), a fourth-order fixed-point iteration with certified convergence. By combining SNM with improved asymptotic approximations and sharp tail bounds, the method achieves superior accuracy and robustness—especially for high-precision applications—outperforming standard Newton-based approaches used in software like R.

ABSTRACT

Accurate and efficient algorithms for the inversion of the cumulative central beta distribution are described. The algorithms are based on the combination of a fourth-order fixed point method with good non-local convergence properties (the Schwarzian-Newton method), asymptotic inversion methods and sharp bounds in the tails of the distribution function.

Motivation & Objective

• To develop more accurate and efficient algorithms for inverting the cumulative central beta distribution, a fundamental operation in statistical computing.
• To overcome convergence and accuracy limitations of existing Newton-based methods used in software such as R, particularly in extreme parameter regimes.
• To provide a reliable, certified-convergence algorithm for quantile computation in the beta distribution, which underpins other common distributions like F and Student’s t.
• To improve initial value estimation for the inversion process using advanced asymptotic expansions and sharp tail bounds, reducing iterations and enhancing speed.
• To enable high-precision computation (near double and single precision) across the full parameter space of the beta distribution, including extreme values of p and q.

Proposed method

• Employ the Schwarzian-Newton method (SNM), a fourth-order fixed-point iteration with superior non-local convergence properties, as the core iterative solver for $ I_x(p,q) = \alpha $.
• Use uniform asymptotic expansions in terms of the error function and the incomplete gamma function to generate high-quality initial approximations for large $ p $ or $ q $.
• Apply sharp upper and lower bounds for the tails of the beta distribution to directly estimate solutions when $ \alpha \to 0 $ or $ \alpha \to 1 $, avoiding iterative refinement.
• Leverage the symmetry relation $ I_x(p,q) = 1 - I_{1-x}(q,p) $ to reduce computation in the upper tail by transforming to the lower tail.
• Design two algorithmic schemes: one for double-precision accuracy using SNM with optimized starting values, and another for single-precision using asymptotic approximations directly when sufficient.
• Ensure certified convergence by combining SNM’s robustness with accurate initial estimates derived from asymptotic theory and tail bounds.

Experimental results

Research questions

• RQ1 Can the Schwarzian-Newton method replace Newton’s method in beta distribution inversion to achieve faster, more reliable convergence?
• RQ2 How can asymptotic approximations for large $ p $ or $ q $ be used to improve initial value estimation and reduce iteration count?
• RQ3 What are the tightest practical bounds for the tails of the beta distribution that can be used to directly estimate quantiles without iteration?
• RQ4 How does the performance of the proposed algorithm compare to existing implementations (e.g., in R) in terms of accuracy and convergence robustness?
• RQ5 Can a single algorithm be designed to handle the full range of $ p, q > 0 $, including extreme values, with guaranteed convergence and high precision?

Key findings

• The Schwarzian-Newton method achieves fourth-order convergence, significantly outperforming Newton’s method (quadratic convergence), especially in regions with poor initial estimates.
• For double-precision accuracy, the proposed algorithm using SNM with tailored initial values (asymptotic or bound-based) ensures convergence across all parameter regimes, including extreme values of $ p $ and $ q $.
• When $ \alpha \leq 0.01 $, the upper tail bound from §2.3 provides a valid solution for $ p < 0.3 $, eliminating the need for iteration in this region.
• For $ p > 30 $ and $ q < 0.5 $, the SNM with bound-based starting values achieves reliable inversion, even when standard methods fail.
• In the central region ($ 0.01 < \alpha \leq 0.5 $), asymptotic approximations based on the error function and incomplete gamma function provide sufficiently accurate starting points to reduce iterations.
• For single-precision applications, direct use of asymptotic approximations (e.g., in terms of the error function or gamma function) is sufficient in many cases, reducing computational cost without sacrificing accuracy.