Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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]

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.

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
1  
@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 –  aukie 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
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 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]
   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.

share|improve this answer

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

Your Answer

 
discard

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.