numerics - How to compute the inverse CDF of HyperbolicDistribution properly?

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]

But Mathematica instead returns

-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.

Answer

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.