Sign up ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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


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


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.

share|improve this question
Looks like something wonky in the numerical evaluation... I'll note that Quantile[] also returns the same erroneous value. – J. M. 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 – Sasha Nov 9 '12 at 14:01

2 Answers 2

up vote 12 down vote accepted

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]

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}]]

However, the value thus returned is not sufficiently accurate:

CDF[dist, x0]

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

which gives a much better result:

CDF[dist, xx]

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"]]

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

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

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.

share|improve this answer

Works with version even without higher precision


"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 *)
share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.