Author: tn
Date: Mon May 21 19:55:30 2012
New Revision: 1341171
URL: http://svn.apache.org/viewvc?rev=1341171&view=rev
Log:
[MATH-718] Use modified Lentz-Thompson algorithm for continued fraction
evaluation.
Modified:
commons/proper/math/trunk/src/changes/changes.xml
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/util/ContinuedFraction.java
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/distribution/BinomialDistributionTest.java
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/distribution/FDistributionTest.java
Modified: commons/proper/math/trunk/src/changes/changes.xml
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/changes/changes.xml?rev=1341171&r1=1341170&r2=1341171&view=diff
==============================================================================
--- commons/proper/math/trunk/src/changes/changes.xml (original)
+++ commons/proper/math/trunk/src/changes/changes.xml Mon May 21 19:55:30 2012
@@ -52,6 +52,10 @@ If the output is not quite correct, chec
<body>
<release version="3.1" date="TBD" description="
">
+ <action dev="tn" type="fix" issue="MATH-718" >
+ Use modified Lentz-Thompson algorithm for continued fraction
evaluation to avoid
+ underflows.
+ </action>
<action dev="luc" type="fix" issue="MATH-780" >
Fixed a wrong assumption on BSP tree attributes when boundary
collapses to a too
small polygon at a non-leaf node.
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/util/ContinuedFraction.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/util/ContinuedFraction.java?rev=1341171&r1=1341170&r2=1341171&view=diff
==============================================================================
---
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/util/ContinuedFraction.java
(original)
+++
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/util/ContinuedFraction.java
Mon May 21 19:55:30 2012
@@ -101,19 +101,18 @@ public abstract class ContinuedFraction
* </p>
*
* <p>
- * The implementation of this method is based on equations 14-17 of:
+ * The implementation of this method is based on the modified Lentz
algorithm as described
+ * on page 18 ff. in:
* <ul>
* <li>
- * Eric W. Weisstein. "Continued Fraction." From MathWorld--A Wolfram Web
- * Resource. <a target="_blank"
- * href="http://mathworld.wolfram.com/ContinuedFraction.html";>
- * http://mathworld.wolfram.com/ContinuedFraction.html</a>
+ * I. J. Thompson, A. R. Barnett. "Coulomb and Bessel Functions of
Complex Arguments and Order."
+ * <a target="_blank"
href="http://www.fresco.org.uk/papers/Thompson-JCP64p490.pdf";>
+ * http://www.fresco.org.uk/papers/Thompson-JCP64p490.pdf</a>
* </li>
* </ul>
- * The recurrence relationship defined in those equations can result in
- * very large intermediate results which can result in numerical overflow.
- * As a means to combat these overflow conditions, the intermediate results
- * are scaled whenever they threaten to become numerically unstable.</p>
+ * Note: the implementation uses the terms a<sub>i</sub> and b<sub>i</sub>
as defined in
+ * <a href="http://mathworld.wolfram.com/ContinuedFraction.html";>Continued
Fraction / MathWorld</a>.
+ * </p>
*
* @param x the evaluation point.
* @param epsilon maximum error allowed.
@@ -122,72 +121,53 @@ public abstract class ContinuedFraction
* @throws ConvergenceException if the algorithm fails to converge.
*/
public double evaluate(double x, double epsilon, int maxIterations) {
- double p0 = 1.0;
- double p1 = getA(0, x);
- double q0 = 0.0;
- double q1 = 1.0;
- double c = p1 / q1;
- int n = 0;
- double relativeError = Double.MAX_VALUE;
- while (n < maxIterations && relativeError > epsilon) {
- ++n;
- double a = getA(n, x);
- double b = getB(n, x);
- double p2 = a * p1 + b * p0;
- double q2 = a * q1 + b * q0;
- boolean infinite = false;
- if (Double.isInfinite(p2) || Double.isInfinite(q2)) {
- /*
- * Need to scale. Try successive powers of the larger of a or b
- * up to 5th power. Throw ConvergenceException if one or both
- * of p2, q2 still overflow.
- */
- double scaleFactor = 1d;
- double lastScaleFactor = 1d;
- final int maxPower = 5;
- final double scale = FastMath.max(a,b);
- if (scale <= 0) { // Can't scale
- throw new
ConvergenceException(LocalizedFormats.CONTINUED_FRACTION_INFINITY_DIVERGENCE,
- x);
- }
- infinite = true;
- for (int i = 0; i < maxPower; i++) {
- lastScaleFactor = scaleFactor;
- scaleFactor *= scale;
- if (a != 0.0 && a > b) {
- p2 = p1 / lastScaleFactor + (b / scaleFactor * p0);
- q2 = q1 / lastScaleFactor + (b / scaleFactor * q0);
- } else if (b != 0) {
- p2 = (a / scaleFactor * p1) + p0 / lastScaleFactor;
- q2 = (a / scaleFactor * q1) + q0 / lastScaleFactor;
- }
- infinite = Double.isInfinite(p2) || Double.isInfinite(q2);
- if (!infinite) {
- break;
- }
- }
- }
+ final double small = 1e-50;
+ double hPrev = getA(0, x);
- if (infinite) {
- // Scaling failed
- throw new
ConvergenceException(LocalizedFormats.CONTINUED_FRACTION_INFINITY_DIVERGENCE,
- x);
+ // use the value of small as epsilon criteria for zero checks
+ if (Precision.equals(hPrev, 0.0, small)) {
+ hPrev = small;
+ }
+
+ int n = 1;
+ double dPrev = 0.0;
+ double cPrev = hPrev;
+ double hN = hPrev;
+
+ while (n < maxIterations) {
+ final double a = getA(n, x);
+ final double b = getB(n, x);
+
+ double dN = a + b * dPrev;
+ if (Precision.equals(dN, 0.0, small)) {
+ dN = small;
+ }
+ double cN = a + b / cPrev;
+ if (Precision.equals(cN, 0.0, small)) {
+ cN = small;
}
- double r = p2 / q2;
+ dN = 1 / dN;
+ final double deltaN = cN * dN;
+ hN = hPrev * deltaN;
- if (Double.isNaN(r)) {
+ if (Double.isInfinite(hN)) {
+ throw new
ConvergenceException(LocalizedFormats.CONTINUED_FRACTION_INFINITY_DIVERGENCE,
+ x);
+ }
+ if (Double.isNaN(hN)) {
throw new
ConvergenceException(LocalizedFormats.CONTINUED_FRACTION_NAN_DIVERGENCE,
x);
}
- relativeError = FastMath.abs(r / c - 1.0);
- // prepare for next iteration
- c = p2 / q2;
- p0 = p1;
- p1 = p2;
- q0 = q1;
- q1 = q2;
+ if (FastMath.abs(deltaN - 1.0) < epsilon) {
+ break;
+ }
+
+ dPrev = dN;
+ cPrev = cN;
+ hPrev = hN;
+ n++;
}
if (n >= maxIterations) {
@@ -195,6 +175,7 @@ public abstract class ContinuedFraction
maxIterations, x);
}
- return c;
+ return hN;
}
+
}
Modified:
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/distribution/BinomialDistributionTest.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/distribution/BinomialDistributionTest.java?rev=1341171&r1=1341170&r2=1341171&view=diff
==============================================================================
---
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/distribution/BinomialDistributionTest.java
(original)
+++
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/distribution/BinomialDistributionTest.java
Mon May 21 19:55:30 2012
@@ -129,4 +129,17 @@ public class BinomialDistributionTest ex
Assert.assertEquals(dist.getNumericalVariance(), 30d * 0.3d * (1d -
0.3d), tol);
}
+ @Test
+ public void testMath718() {
+ // for large trials the evaluation of ContinuedFraction was inaccurate
+ // do a sweep over several large trials to test if the current
implementation is
+ // numerically stable.
+
+ for (int trials = 500000; trials < 20000000; trials += 100000) {
+ BinomialDistribution dist = new BinomialDistribution(trials, 0.5);
+ int p = dist.inverseCumulativeProbability(0.5);
+ Assert.assertEquals(trials / 2, p);
+ }
+
+ }
}
Modified:
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/distribution/FDistributionTest.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/distribution/FDistributionTest.java?rev=1341171&r1=1341170&r2=1341171&view=diff
==============================================================================
---
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/distribution/FDistributionTest.java
(original)
+++
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/distribution/FDistributionTest.java
Mon May 21 19:55:30 2012
@@ -142,4 +142,18 @@ public class FDistributionTest extends R
Assert.assertEquals(dist.getNumericalMean(), 5d / (5d - 2d), tol);
Assert.assertEquals(dist.getNumericalVariance(), (2d * 5d * 5d * 4d) /
9d, tol);
}
+
+ @Test
+ public void testMath785() {
+ // this test was failing due to inaccurate results from
ContinuedFraction.
+
+ try {
+ double prob = 0.01;
+ FDistribution f = new FDistribution(200000, 200000);
+ double result = f.inverseCumulativeProbability(prob);
+ Assert.assertTrue(result < 1.0);
+ } catch (Exception e) {
+ Assert.fail("Failing to calculate inverse cumulative probability");
+ }
+ }
}
