Efficient and stable computation of inverse CDF

Question

What is the most efficient and numerically stable algorithm for computing the inverse CDF $F^{-1}(y)$ of a probability function, assuming that both the PDF $f(x)$ and the CDF $F(x)$ are known analytically but the inverse CDF is not?

Clearly, this is the same as finding the root of the nonlinear function $G(x) \equiv F(x) - y$, with $x \in \mathbb{R}$ for a given $y \in (0, 1)$. We know/assume that $G(x)$ is:

• monotonically non-decreasing (being $F(x)$ a CDF);
• differentiable at least $k \ge 1$ times (its first derivative is $G^\prime(x) = f(x)$);
• continuous (we assume that $f(x)$ is not a delta function).
• I am particularly interested in the case in which the shape of the PDF is a generic mixture of Gaussian distributions multiplied by a polynomial of arbitrary degree (in this case, the CDF can be computed analytically).

There are several standard methods for root-finding of generic nonlinear functions (see for example Chapter 9 of Numerical Recipes). The method I would use is Brent's method, with perhaps a couple of steps of Newton-Raphson in the end to refine the root (if at all); but I wonder if there are better ways given the assumptions above.

Also, note that I need to compute the inverse CDF for a large number ($\sim 10^5$) of distinct CDFs within this class; a lookup table or pre-computation is not feasible.

Answer 1

The assumptions you have are no more specific than what you need to assume to make something like Newton's method work. In fact, you don't even assume enough to make the problem unique: you only assume that $F(x)$ is non-decreasing, when of course you need it to be monotonically increasing to make finding an answer unique.

As a consequence, there is really nothing very much that could be exploited here, other than the fact that if you assume strict monotonicity, then the root is unique. Thus, what you suggest in terms of methods makes perfect sense.