mirror of
https://github.com/apache/commons-math.git
synced 2025-02-08 02:59:36 +00:00
[MATH-718] Use modified Lentz-Thompson algorithm for continued fraction evaluation.
git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1341171 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Thomas Neidhart
parent
282f7175a4
commit
3a08bfa6d8
@ -52,6 +52,10 @@ If the output is not quite correct, check for invisible trailing spaces!
|
||||
<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.
|
||||
|
||||
@ -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);
|
||||
|
||||
// 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;
|
||||
}
|
||||
|
||||
if (infinite) {
|
||||
// Scaling failed
|
||||
throw new ConvergenceException(LocalizedFormats.CONTINUED_FRACTION_INFINITY_DIVERGENCE,
|
||||
x);
|
||||
dN = 1 / dN;
|
||||
final double deltaN = cN * dN;
|
||||
hN = hPrev * deltaN;
|
||||
|
||||
if (Double.isInfinite(hN)) {
|
||||
throw new ConvergenceException(LocalizedFormats.CONTINUED_FRACTION_INFINITY_DIVERGENCE,
|
||||
x);
|
||||
}
|
||||
|
||||
double r = p2 / q2;
|
||||
|
||||
if (Double.isNaN(r)) {
|
||||
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;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
@ -129,4 +129,17 @@ public class BinomialDistributionTest extends IntegerDistributionAbstractTest {
|
||||
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);
|
||||
}
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
@ -142,4 +142,18 @@ public class FDistributionTest extends RealDistributionAbstractTest {
|
||||
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");
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user