commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From t.@apache.org
Subject svn commit: r1341171 - in /commons/proper/math/trunk/src: changes/ main/java/org/apache/commons/math3/util/ test/java/org/apache/commons/math3/distribution/
Date Mon, 21 May 2012 19:55:31 GMT
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");
+        }
+    }
 }



Mime
View raw message