# Efficient and stable computation of inverse CDF

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.

• If computing $F^{-1}(y)$ is difficult, have you considered precomputing an interpolating function that approximates it well on $(0,1)$, so that repeated evaluations of $F^{-1}(y)$ become much cheaper? That would bypass the issue of root-finding. Also, I believe a typical black-box method is Alefeld-Potra-Shi (e.g. github.com/JuliaLang/Roots.jl), which is perfectly suitable for this problem. I also think this isn't so different from a typical "generic" root-finding problem. Apr 5, 2016 at 22:39
• Thanks for the answer. Precomputation is not feasible -- I added a point to explain why. I didn't know the Alefeld-Potra-Shi method, I will have a look! Apr 5, 2016 at 22:43

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.
• (+1) Thanks. Yeah, I was a bit stingy on the general assumptions... monotonicity is okay (in fact one of the cases I care about is when $f(x)$ is a mixture of Gaussians times a polynomial). Well, good to know that I am not missing anything then. Apr 8, 2016 at 2:42