What numerical methods are used to invert CDF functions that don't have a closed form inverse?

Im currently using numpy.percentile in a Python program and I assume that percentiles for probability distributions are calculated by inverting the Cumulative Distribution Function (CDF) of the random variable. Since many of the standard probability distributions CDFs don't have a closed form inverse (e.g. the normal distribution) I wonder how this inverse is approximated.

I realize there is a lot of literature on the subject and I do want to put in effort. My hope is that there is someone among us who is well versed in numerical analysis and can give a short overview on this matter that makes it easier to find the detailed algorithms (and derivations) and understand their corresponding importance.

The inverse of a CDF is called a quantile function $Q$. Wikipedia goes into a brief overview on how it can be computed. In short, there are 3 major points to be made:

• The derivative of the CDF, the PDF, is usually known and easily computed (compared to the CDF).
• The CDF is strictly increasing with known domain and range.
• The quantile function may be represented in terms of differential equations.

The derivative enables root-finding methods such as Newton's method, which converges rapidly provided a sufficiently close starting point and sufficiently well-behaved CDF.

Since the CDF is strictly increasing with known domain and range, bracketing methods can be used to ensure worst-case convergence if the solution is far from the initial estimate or in a region where Newton's method cannot reasonably approximate the CDF.

Note also that in the cases when the CDF is being numerically computed using the PDF, it is generally cheaper to compute $\operatorname{Pr}(x_{i-1}\le X<x_i)$ than it is to compute $\operatorname{Pr}(X<x_i)$. This makes root-finding algorithms cheaper than they usually are.

Although it is possible to only use the extreme points when considering bracketing methods, when many calculations are desired, results can be pre-tabulated e.g. to find $Q(p_i)$ for $p_i\in\{0,0.1,0.3,0.5,0.7,0.9,1\}$ beforehand (or to find $p_i=\operatorname{Pr}(X<x_i)$ for several $x_i$ to get points $Q(p_i)\equiv x_i$) and then work from there. For a rough approximation of the quantile function, an interpolation may be done instead of trying to solve for $Q(p)$ exactly.

In the analytic cases, one may also interpret the quantile function as a differential equation, either directly from the PDF or indirectly from the CDF. This approach opens up the possibility for a variety of methods, such as series expansions.