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.

Looks like something wonky in the numerical evaluation... I'll note that

`Quantile[]`

also returns the same erroneous value. – J. M.'s ennui – 2012-10-29T03:40:29.4001@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 – 2012-10-29T05:17:05.077

@aukie looks like a bug. Works for x=-1 or 1 though. – chris – 2012-10-29T08:42:35.647

2@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 – 2012-11-09T14:01:12.513