[MATH-785] Numerical Underflow in ContinuedFraction

• Type: Bug
• Status: Closed
• Priority: Major
• Resolution: Fixed
• Affects Version/s: 3.0
• Fix Version/s: 3.1

Environment

Issue seen in both 3.0 release binary version as well as a fresh checkout of the subversion trunk.

java -version output:

java version "1.6.0_26"
Java(TM) SE Runtime Environment (build 1.6.0_26-b03)
Java HotSpot(TM) 64-Bit Server VM (build 20.1-b02, mixed mode)

(On Ubuntu 12.04)

• Assignee: Unassigned
• Reporter: Colin J. Fuller
• Created: 01/May/12 02:54
• Updated: 04/Mar/13 18:53
• Resolved: 21/May/12 20:01

Description

The ContinuedFraction calculation can underflow in the evaluate method, similar to the overflow case already dealt with. I encountered this problem while trying to evaluate the inverse cumulative probability of an F distribution with a large number of degrees of freedom.

I would guess this has the same cause as MATH-718 and MATH-738, though I am not experiencing inaccurate results but rather an exception.

For instance, the following test case fails:

double prob = 0.01;
FDistribution f = new FDistribution(200000, 200000);
double fails = f.inverseCumulativeProbability(prob);

This produces a NoBracketingException with the following stack trace:

org.apache.commons.math3.exception.NoBracketingException: function values at endpoints do not have different signs, endpoints: [0, 1], values: [-0.01, -∞]
at org.apache.commons.math3.analysis.solvers.BrentSolver.doSolve(BrentSolver.java:118)
at org.apache.commons.math3.analysis.solvers.BaseAbstractUnivariateSolver.solve(BaseAbstractUnivariateSolver.java:190)
at org.apache.commons.math3.analysis.solvers.BaseAbstractUnivariateSolver.solve(BaseAbstractUnivariateSolver.java:195)
at org.apache.commons.math3.analysis.solvers.UnivariateSolverUtils.solve(UnivariateSolverUtils.java:77)
at org.apache.commons.math3.distribution.AbstractRealDistribution.inverseCumulativeProbability(AbstractRealDistribution.java:156)

I could avoid the issue as in the comment to MATH-718 by relaxing the default value of epsilon in ContinuedFraction, although in my test case I can't see any reason the current default precision shouldn't be attainable.

I fixed the issue by implementing underflow detection in ContinuedFraction and rescaling to larger values similarly to how the overflow detection that is already there works. I will attach a patch shortly.

One possible issue with this fix is that if there exists a case where there is a legitimate reason for p2 or q2 to be zero (I cannot think of one), it might break that case.

Issue Links

• is related to: MATH-718 inverseCumulativeProbability of BinomialDistribution returns wrong value for large trials.