# How to compute the inverse CDF properly?

Note: this has been fixed in version 9.

I want to compute the CDF and inverse CDF of the hyperbolic distribution:

α = 2; β = 3/2;
x = -3;
u = N[CDF[HyperbolicDistribution[α, β, 1, 0], x], 30]


The output of the above code is

5.36058384200167863956651004148*10^-6


A small, but not too small number. Now, since the CDF and inverse CDF are compositional inverses, I expect the output from the following line to be approximately $-3$:

xA = N[InverseCDF[HyperbolicDistribution[α, β, 1, 0], u], 30]


-23.0838707895835608942828900958


which is just way off. So, am I doing something wrong here? Or are Mathematica's CDF[] and InverseCDF[] functions not to be trusted for this distribution?

When the above experiment is repeated with the normal distribution, there are no problems:

x = -12;
u = N[CDF[NormalDistribution[0, 1], x], 30]
xA = N[InverseCDF[NormalDistribution[0, 1], u], 30]


The output is $-12$ as expected.

-
Looks like something wonky in the numerical evaluation... I'll note that Quantile[] also returns the same erroneous value. –  Ｊ. Ｍ. Oct 29 '12 at 3:40
@J. M. I have a suspicion this is a bug in Mathemica, I have observed this inconsistency with the cdf and quantile for other distributions also. StudentT works as expected though –  mark Oct 29 '12 at 5:17
@aukie looks like a bug. Works for x=-1 or 1 though. –  chris Oct 29 '12 at 8:42
@aukie This bug has been fixed. It would be helpful if you could provide additional instances of incorrect behavior you mention you encountered and sent them to support@wolfram.com –  Sasha Nov 9 '12 at 14:01

One hack-ish method for evaluating an inverse CDF is to use the event location functionality of NDSolve[]. As an example:

dist = HyperbolicDistribution[2, 3/2, 1, 0];
c0 = N[CDF[dist, 0], 20]
0.058032099055437685722


Suppose that we want to evaluate the inverse CDF for this particular hyperbolic distribution at the $p$-value $7\times10^{-6}$. Since this value is smaller than c0, we integrate the ODE for the CDF,

$$c^\prime(x)=\mathtt{PDF[dist},\;x\mathtt{]},\quad c^\prime(0)=\mathtt{c0}$$

backwards, like so:

With[{p = 7*^-6},
x0 = Extract[c["Domain"] /.
First @ NDSolve[{c'[u] == PDF[dist, u], c[0] == c0}, c, {u, -Infinity, 0},
AccuracyGoal -> Infinity,
Method -> {"EventLocator", "Event" -> c[u] == p,
Method -> "StiffnessSwitching"}], {1, 1}]]
-2.9217514743315895


However, the value thus returned is not sufficiently accurate:

CDF[dist, x0]
7.00002*10^-6


The cure is easy: treat x0 as a starting value for further polishing with Newton-Raphson:

With[{p = 7*^-6},
xx = FixedPoint[(# - (CDF[dist, #] - p)/PDF[dist, #]) &, x0]]
-2.9217521320606825


which gives a much better result:

CDF[dist, xx]
7.*10^-6


If one wanted to evaluate the inverse CDF for values greater than c0, one should of course integrate forward (that is, change the integration interval to {u, 0, Infinity}) and then extract the right end of the "Domain" of the InterpolatingFunction[] (using, say, Extract[(* stuff *), {1, 2}]).

Another ODE-based method that can be used is based on the usual formula for the derivative of an inverse function:

$$\frac{\mathrm d}{\mathrm dx}f^{(-1)}(x)=\frac1{f^\prime(f^{(-1)}(x))}$$

which can be treated as an ODE for the inverse CDF. For the case of the hyperbolic distribution, Leobacher and Pillichshammer, in their paper, suggest the use of the initial condition

$$f^{(-1)}(\mathtt{CDF[HyperbolicDistribution[}\alpha,\beta,\delta,\mu\mathtt{]},x_0\mathtt{]})=x_0,\qquad x_0=\mu+\frac{\beta\delta}{\sqrt{\alpha^2-\beta^2}}$$

Here's a Mathematica demonstration of this strategy, using the $p$-value $2\times10^{-8}$:

dist = HyperbolicDistribution[2, 3/2, 1, 0];
With[{α = 2, β = 3/2, δ = 1, μ = 0}, xs = N[μ + β δ/Sqrt[α^2 - β^2], 25]];
cs = CDF[dist, xs];

With[{p = 2*^-8},
x0 = c[p] /.
First @ NDSolve[{c'[u] == 1/PDF[dist, c[u]], c[cs] == xs}, c, {u, cs, p},
AccuracyGoal -> Infinity, Method -> "StiffnessSwitching"]]
-4.630068132899461


One can again use Newton-Raphson to further polish the value obtained:

With[{p = 2*^-8}, xx = FixedPoint[(# - (CDF[dist, #] - p)/PDF[dist, #]) &, x0]]
-4.62518360657725


This strategy will run into numerical difficulties for $p$-values very near to $0$ or very near to $1$; see the Leobacher/Pillichshammer paper for details on how to proceed in these cases.

Encapsulating any of these two strategies into a general procedure is left as an exercise to the interested reader.

-

Works with version 9.0.1.0 even without higher precision

\$Version

"9.0 for Mac OS X x86 (64-bit) (January 24, 2013)"

a = 2; b = 3/2;
x = -3;
u = CDF[HyperbolicDistribution[a, b, 1, 0], x] // N

(* 5.360583842001691*^-6 *)

x == InverseCDF[HyperbolicDistribution[a, b, 1, 0], u]

(* True *)

-