brentworden 2004/06/10 11:27:47
Modified: math/src/test/org/apache/commons/math/distribution
FDistributionTest.java
math/src/java/org/apache/commons/math/special Beta.java
Log:
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.
Revision Changes Path
1.15 +9 -1
jakarta-commons/math/src/test/org/apache/commons/math/distribution/FDistributionTest.java
Index: FDistributionTest.java
===================================================================
RCS file:
/home/cvs/jakarta-commons/math/src/test/org/apache/commons/math/distribution/FDistributionTest.java,v
retrieving revision 1.14
retrieving revision 1.15
diff -u -r1.14 -r1.15
--- FDistributionTest.java 30 May 2004 01:39:33 -0000 1.14
+++ FDistributionTest.java 10 Jun 2004 18:27:47 -0000 1.15
@@ -97,4 +97,12 @@
}
}
+ 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);
+ }
}
1.20 +12 -28
jakarta-commons/math/src/java/org/apache/commons/math/special/Beta.java
Index: Beta.java
===================================================================
RCS file:
/home/cvs/jakarta-commons/math/src/java/org/apache/commons/math/special/Beta.java,v
retrieving revision 1.19
retrieving revision 1.20
diff -u -r1.19 -r1.20
--- Beta.java 27 Apr 2004 04:37:59 -0000 1.19
+++ Beta.java 10 Jun 2004 18:27:47 -0000 1.20
@@ -122,48 +122,32 @@
(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;
---------------------------------------------------------------------
To unsubscribe, e-mail: [EMAIL PROTECTED]
For additional commands, e-mail: [EMAIL PROTECTED]
