PR: 29414

I changed the continued fraction used in regularizedBeta resulting in faster 
convergence.  I added the test case provided by scott and ran all units tests 
with all of them passing.


git-svn-id: https://svn.apache.org/repos/asf/jakarta/commons/proper/math/trunk@141289 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Brent Worden 2004-06-10 18:27:47 +00:00
parent 15f4e86e93
commit d83250a890
2 changed files with 21 additions and 29 deletions

View File

@ -24,7 +24,7 @@ import org.apache.commons.math.util.ContinuedFraction;
* This is a utility class that provides computation methods related to the
* Beta family of functions.
*
* @version $Revision: 1.19 $ $Date: 2004/04/27 04:37:59 $
* @version $Revision: 1.20 $ $Date: 2004/06/10 18:27:47 $
*/
public class Beta implements Serializable {
/** Maximum allowed numerical error. */
@ -122,48 +122,32 @@ public class Beta implements Serializable {
(x > 1) || (a <= 0.0) || (b <= 0.0))
{
ret = Double.NaN;
} else if (x > (a + 1.0) / (a + b + 1.0)) {
} else if (x > (a + 1.0) / (a + b + 2.0)) {
ret = 1.0 - regularizedBeta(1.0 - x, b, a, epsilon, maxIterations);
} else {
ContinuedFraction fraction = new ContinuedFraction() {
protected double getB(int n, double x) {
double ret;
double m;
switch (n) {
case 1 :
ret = 1.0;
break;
default :
if (n % 2 == 0) { // even
m = (n - 2.0) / 2.0;
ret = -((a + m) * (a + b + m) * x) /
((a + (2 * m)) * (a + (2 * m) + 1.0));
} else {
m = (n - 1.0) / 2.0;
ret = (m * (b - m) * x) /
((a + (2 * m) - 1) * (a + (2 * m)));
}
break;
if (n % 2 == 0) { // even
m = n / 2.0;
ret = (m * (b - m) * x) /
((a + (2 * m) - 1) * (a + (2 * m)));
} else {
m = (n - 1.0) / 2.0;
ret = -((a + m) * (a + b + m) * x) /
((a + (2 * m)) * (a + (2 * m) + 1.0));
}
return ret;
}
protected double getA(int n, double x) {
double ret;
switch (n) {
case 0 :
ret = 0.0;
break;
default :
ret = 1.0;
break;
}
return ret;
return 1.0;
}
};
ret = Math.exp((a * Math.log(x)) + (b * Math.log(1.0 - x)) -
Math.log(a) - logBeta(a, b, epsilon, maxIterations)) *
fraction.evaluate(x, epsilon, maxIterations);
1.0 / fraction.evaluate(x, epsilon, maxIterations);
}
return ret;

View File

@ -21,7 +21,7 @@ package org.apache.commons.math.distribution;
* Extends ContinuousDistributionAbstractTest. See class javadoc for
* ContinuousDistributionAbstractTest for details.
*
* @version $Revision: 1.14 $ $Date: 2004/05/30 01:39:33 $
* @version $Revision: 1.15 $ $Date: 2004/06/10 18:27:47 $
*/
public class FDistributionTest extends ContinuousDistributionAbstractTest {
@ -97,4 +97,12 @@ public class FDistributionTest extends ContinuousDistributionAbstractTest {
}
}
public void testLargeDegreesOfFreedom() throws Exception {
org.apache.commons.math.distribution.FDistributionImpl fd =
new org.apache.commons.math.distribution.FDistributionImpl(
100000., 100000.);
double p = fd.cumulativeProbability(.999);
double x = fd.inverseCumulativeProbability(p);
assertEquals(.999, x, 1.0e-5);
}
}