You are viewing a plain text version of this content. The canonical link for it is here.
Posted to dev@commons.apache.org by br...@apache.org on 2004/06/07 22:30:16 UTC
cvs commit: jakarta-commons/math/src/java/org/apache/commons/math/special Gamma.java
brentworden 2004/06/07 13:30:16
Modified: math/src/test/org/apache/commons/math/special GammaTest.java
math/src/test/org/apache/commons/math/stat/inference
ChiSquareTestTest.java
math/src/java/org/apache/commons/math/special Gamma.java
Log:
PR: 29419
Added an implementation of regularized gamma function, Q(a, x) = 1 - P(a,x), based on a continued fraction. This converges much faster for the large x case. I added the example submitted by Scott as a test case and ran all the test cases with everything passing.
Revision Changes Path
1.10 +5 -3 jakarta-commons/math/src/test/org/apache/commons/math/special/GammaTest.java
Index: GammaTest.java
===================================================================
RCS file: /home/cvs/jakarta-commons/math/src/test/org/apache/commons/math/special/GammaTest.java,v
retrieving revision 1.9
retrieving revision 1.10
diff -u -r1.9 -r1.10
--- GammaTest.java 21 Feb 2004 21:35:17 -0000 1.9
+++ GammaTest.java 7 Jun 2004 20:30:16 -0000 1.10
@@ -34,8 +34,10 @@
private void testRegularizedGamma(double expected, double a, double x) {
try {
- double actual = Gamma.regularizedGammaP(a, x);
- TestUtils.assertEquals(expected, actual, 10e-5);
+ double actualP = Gamma.regularizedGammaP(a, x);
+ double actualQ = Gamma.regularizedGammaQ(a, x);
+ TestUtils.assertEquals(expected, actualP, 10e-5);
+ TestUtils.assertEquals(actualP, 1.0 - actualQ, 10e-5);
} catch(MathException ex){
fail(ex.getMessage());
}
1.2 +17 -2 jakarta-commons/math/src/test/org/apache/commons/math/stat/inference/ChiSquareTestTest.java
Index: ChiSquareTestTest.java
===================================================================
RCS file: /home/cvs/jakarta-commons/math/src/test/org/apache/commons/math/stat/inference/ChiSquareTestTest.java,v
retrieving revision 1.1
retrieving revision 1.2
diff -u -r1.1 -r1.2
--- ChiSquareTestTest.java 3 May 2004 03:04:54 -0000 1.1
+++ ChiSquareTestTest.java 7 Jun 2004 20:30:16 -0000 1.2
@@ -126,5 +126,20 @@
} catch (IllegalArgumentException ex) {
// expected
}
- }
+ }
+
+ public void testChiSquareLargeTestStatistic() throws Exception {
+ double[] exp = new double[] {
+ 3389119.5, 649136.6, 285745.4, 25357364.76, 11291189.78, 543628.0,
+ 232921.0, 437665.75
+ };
+
+ long[] obs = new long[] {
+ 2372383, 584222, 257170, 17750155, 7903832, 489265, 209628, 393899
+ };
+ org.apache.commons.math.stat.inference.ChiSquareTestImpl csti =
+ new org.apache.commons.math.stat.inference.ChiSquareTestImpl();
+ double cst = csti.chiSquareTest(exp, obs);
+ assertEquals("chi-square p-value", 0.0, cst, 1E-3);
+ }
}
1.19 +120 -27 jakarta-commons/math/src/java/org/apache/commons/math/special/Gamma.java
Index: Gamma.java
===================================================================
RCS file: /home/cvs/jakarta-commons/math/src/java/org/apache/commons/math/special/Gamma.java,v
retrieving revision 1.18
retrieving revision 1.19
diff -u -r1.18 -r1.19
--- Gamma.java 23 Apr 2004 19:30:47 -0000 1.18
+++ Gamma.java 7 Jun 2004 20:30:16 -0000 1.19
@@ -19,6 +19,7 @@
import org.apache.commons.math.ConvergenceException;
import org.apache.commons.math.MathException;
+import org.apache.commons.math.util.ContinuedFraction;
/**
* This is a utility class that provides computation methods related to the
@@ -26,7 +27,8 @@
*
* @version $Revision$ $Date$
*/
-public class Gamma implements Serializable{
+public class Gamma implements Serializable {
+
/** Maximum allowed numerical error. */
private static final double DEFAULT_EPSILON = 10e-9;
@@ -59,6 +61,45 @@
}
/**
+ * Returns the natural logarithm of the gamma function Γ(x).
+ *
+ * The implementation of this method is based on:
+ * <ul>
+ * <li><a href="http://mathworld.wolfram.com/GammaFunction.html">
+ * Gamma Function</a>, equation (28).</li>
+ * <li><a href="http://mathworld.wolfram.com/LanczosApproximation.html">
+ * Lanczos Approximation</a>, equations (1) through (5).</li>
+ * <li><a href="http://my.fit.edu/~gabdo/gamma.txt">Paul Godfrey, A note on
+ * the computation of the convergent Lanczos complex Gamma approximation
+ * </a></li>
+ * </ul>
+ *
+ * @param x the value.
+ * @return log(Γ(x))
+ */
+ public static double logGamma(double x) {
+ double ret;
+
+ if (Double.isNaN(x) || (x <= 0.0)) {
+ ret = Double.NaN;
+ } else {
+ double g = 607.0 / 128.0;
+
+ double sum = 0.0;
+ for (int i = 1; i < lanczos.length; ++i) {
+ sum = sum + (lanczos[i] / (x + i));
+ }
+ sum = sum + lanczos[0];
+
+ double tmp = x + g + .5;
+ ret = ((x + .5) * Math.log(tmp)) - tmp +
+ (.5 * Math.log(2.0 * Math.PI)) + Math.log(sum) - Math.log(x);
+ }
+
+ return ret;
+ }
+
+ /**
* Returns the regularized gamma function P(a, x).
*
* @param a the a parameter.
@@ -71,7 +112,8 @@
{
return regularizedGammaP(a, x, DEFAULT_EPSILON, Integer.MAX_VALUE);
}
-
+
+
/**
* Returns the regularized gamma function P(a, x).
*
@@ -110,6 +152,10 @@
ret = Double.NaN;
} else if (x == 0.0) {
ret = 0.0;
+ } else if (a > 1.0 && x > a) {
+ // use regularizedGammaQ because it should converge faster in this
+ // case.
+ ret = 1.0 - regularizedGammaQ(a, x, epsilon, maxIterations);
} else {
// calculate series
double n = 0.0; // current element index
@@ -133,41 +179,88 @@
return ret;
}
-
+
/**
- * Returns the natural logarithm of the gamma function Γ(x).
- *
+ * Returns the regularized gamma function Q(a, x) = 1 - P(a, x).
+ *
+ * @param a the a parameter.
+ * @param x the value.
+ * @return the regularized gamma function Q(a, x)
+ * @throws MathException if the algorithm fails to converge.
+ */
+ public static double regularizedGammaQ(double a, double x)
+ throws MathException
+ {
+ return regularizedGammaQ(a, x, DEFAULT_EPSILON, Integer.MAX_VALUE);
+ }
+
+ /**
+ * Returns the regularized gamma function Q(a, x) = 1 - P(a, x).
+ *
* The implementation of this method is based on:
* <ul>
- * <li><a href="http://mathworld.wolfram.com/GammaFunction.html">
- * Gamma Function</a>, equation (28).</li>
- * <li><a href="http://mathworld.wolfram.com/LanczosApproximation.html">
- * Lanczos Approximation</a>, equations (1) through (5).</li>
- * <li><a href="http://my.fit.edu/~gabdo/gamma.txt">Paul Godfrey, A note on
- * the computation of the convergent Lanczos complex Gamma approximation
- * </a></li>
+ * <li>
+ * <a href="http://mathworld.wolfram.com/RegularizedGammaFunction.html">
+ * Regularized Gamma Function</a>, equation (1).</li>
+ * <li>
+ * <a href=" http://functions.wolfram.com/GammaBetaErf/GammaRegularized/10/0003/">
+ * Regularized incomplete gamma function: Continued fraction representations (formula 06.08.10.0003)</a></li>
* </ul>
*
+ * @param a the a parameter.
* @param x the value.
- * @return log(Γ(x))
+ * @param epsilon When the absolute value of the nth item in the
+ * series is less than epsilon the approximation ceases
+ * to calculate further elements in the series.
+ * @param maxIterations Maximum number of "iterations" to complete.
+ * @return the regularized gamma function P(a, x)
+ * @throws MathException if the algorithm fails to converge.
*/
- public static double logGamma(double x) {
+ public static double regularizedGammaQ(final double a,
+ double x,
+ double epsilon,
+ int maxIterations)
+ throws MathException
+ {
double ret;
- if (Double.isNaN(x) || (x <= 0.0)) {
+ if (Double.isNaN(a) || Double.isNaN(x) || (a <= 0.0) || (x < 0.0)) {
ret = Double.NaN;
+ } else if (x == 0.0) {
+ ret = 1.0;
+ } else if (x < a || a <= 1.0) {
+ // use regularizedGammaP because it should converge faster in this
+ // case.
+ ret = 1.0 - regularizedGammaP(a, x, epsilon, maxIterations);
} else {
- double g = 607.0 / 128.0;
-
- double sum = 0.0;
- for (int i = 1; i < lanczos.length; ++i) {
- sum = sum + (lanczos[i] / (x + i));
- }
- sum = sum + lanczos[0];
-
- double tmp = x + g + .5;
- ret = ((x + .5) * Math.log(tmp)) - tmp +
- (.5 * Math.log(2.0 * Math.PI)) + Math.log(sum) - Math.log(x);
+ // create continued fraction
+ ContinuedFraction cf = new ContinuedFraction() {
+ protected double getA(int n, double x) {
+ double ret;
+ switch(n) {
+ case 0: ret = 0.0; break;
+ default:
+ ret = ((2.0 * n) - 1.0) - a + x; break;
+ }
+ return ret;
+ }
+
+ protected double getB(int n, double x) {
+ double ret;
+ double t;
+ switch(n) {
+ case 1: ret = 1.0; break;
+ default:
+ t = n - 1.0;
+ ret = t * (a - t);
+ break;
+ }
+ return ret;
+ }
+ };
+
+ ret = cf.evaluate(x, epsilon, maxIterations);
+ ret = Math.exp(-x + (a * Math.log(x)) - logGamma(a)) * ret;
}
return ret;
---------------------------------------------------------------------
To unsubscribe, e-mail: commons-dev-unsubscribe@jakarta.apache.org
For additional commands, e-mail: commons-dev-help@jakarta.apache.org