Author: celestin
Date: Mon Nov 19 20:07:25 2012
New Revision: 1411387
URL: http://svn.apache.org/viewvc?rev=1411387&view=rev
Log:
In class Gamma, implementation of (a, b) -> log(Gamma(a + b)).
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/special/Gamma.java
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/special/GammaTest.java
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/special/Gamma.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/special/Gamma.java?rev=1411387&r1=1411386&r2=1411387&view=diff
==============================================================================
---
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/special/Gamma.java
(original)
+++
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/special/Gamma.java
Mon Nov 19 20:07:25 2012
@@ -19,6 +19,7 @@ package org.apache.commons.math3.special
import org.apache.commons.math3.exception.MaxCountExceededException;
import org.apache.commons.math3.exception.NumberIsTooLargeException;
import org.apache.commons.math3.exception.NumberIsTooSmallException;
+import org.apache.commons.math3.exception.OutOfRangeException;
import org.apache.commons.math3.util.ContinuedFraction;
import org.apache.commons.math3.util.FastMath;
@@ -701,4 +702,33 @@ public class Gamma {
}
return ret;
}
+
+ /**
+ * Returns the value of log Î(a + b) for 1 ⤠a, b ⤠2.
+ *
+ * @param a First argument.
+ * @param b Second argument.
+ * @return the value of {@code log(Gamma(a + b))}.
+ * @throws OutOfRangeException if {@code a} or {@code b} is lower than
+ * {@code 1.0} or greater than {@code 2.0}.
+ */
+ public static double logGammaSum(final double a, final double b)
+ throws OutOfRangeException {
+
+ if ((a < 1.0) || (a > 2.0)) {
+ throw new OutOfRangeException(a, 1.0, 2.0);
+ }
+ if ((b < 1.0) || (b > 2.0)) {
+ throw new OutOfRangeException(b, 1.0, 2.0);
+ }
+
+ final double x = a + b - 2.0;
+ if (x <= 0.25) {
+ return Gamma.logGamma1p(1.0 + x);
+ } else if (x <= 1.25) {
+ return Gamma.logGamma1p(x) + FastMath.log1p(x);
+ } else {
+ return Gamma.logGamma1p(x - 1.0) + FastMath.log(x * (1.0 + x));
+ }
+ }
}
Modified:
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/special/GammaTest.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/special/GammaTest.java?rev=1411387&r1=1411386&r2=1411387&view=diff
==============================================================================
---
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/special/GammaTest.java
(original)
+++
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/special/GammaTest.java
Mon Nov 19 20:07:25 2012
@@ -465,7 +465,7 @@ public class GammaTest {
* print("{", float(x[i]), ",", float(gamma(x[i])), "},");
* </pre>
*/
- private final static double[][] GAMMA_REF = {
+ private static final double[][] GAMMA_REF = {
{ - 19.875 , 4.920331854832504e-18 },
{ - 19.75 , 3.879938752480031e-18 },
{ - 19.625 , 4.323498423815027e-18 },
@@ -954,6 +954,127 @@ public class GammaTest {
}
}
+ /**
+ * Reference data for the {@link Gamma#logGammaSum(double, double)}
+ * function. This data was generated with the following
+ * <a href="http://maxima.sourceforge.net/";>Maxima</a> script.
+ *
+ * <pre>
+ * kill(all);
+ *
+ * fpprec : 64;
+ * gsumln(a, b) := log(gamma(a + b));
+ *
+ * x : [1.0b0, 1.125b0, 1.25b0, 1.375b0, 1.5b0, 1.625b0, 1.75b0, 1.875b0,
2.0b0];
+ *
+ * for i : 1 while i <= length(x) do
+ * for j : 1 while j <= length(x) do block(
+ * a : x[i],
+ * b : x[j],
+ * print("{", float(a), ",", float(b), ",", float(gsumln(a, b)), "},")
+ * );
+ * </pre>
+ */
+ private static final double[][] LOG_GAMMA_SUM_REF = {
+ { 1.0 , 1.0 , 0.0 },
+ { 1.0 , 1.125 , .05775985153034387 },
+ { 1.0 , 1.25 , .1248717148923966 },
+ { 1.0 , 1.375 , .2006984603774558 },
+ { 1.0 , 1.5 , .2846828704729192 },
+ { 1.0 , 1.625 , .3763336820249054 },
+ { 1.0 , 1.75 , .4752146669149371 },
+ { 1.0 , 1.875 , .5809359740231859 },
+ { 1.0 , 2.0 , .6931471805599453 },
+ { 1.125 , 1.0 , .05775985153034387 },
+ { 1.125 , 1.125 , .1248717148923966 },
+ { 1.125 , 1.25 , .2006984603774558 },
+ { 1.125 , 1.375 , .2846828704729192 },
+ { 1.125 , 1.5 , .3763336820249054 },
+ { 1.125 , 1.625 , .4752146669149371 },
+ { 1.125 , 1.75 , .5809359740231859 },
+ { 1.125 , 1.875 , .6931471805599453 },
+ { 1.125 , 2.0 , 0.811531653906724 },
+ { 1.25 , 1.0 , .1248717148923966 },
+ { 1.25 , 1.125 , .2006984603774558 },
+ { 1.25 , 1.25 , .2846828704729192 },
+ { 1.25 , 1.375 , .3763336820249054 },
+ { 1.25 , 1.5 , .4752146669149371 },
+ { 1.25 , 1.625 , .5809359740231859 },
+ { 1.25 , 1.75 , .6931471805599453 },
+ { 1.25 , 1.875 , 0.811531653906724 },
+ { 1.25 , 2.0 , .9358019311087253 },
+ { 1.375 , 1.0 , .2006984603774558 },
+ { 1.375 , 1.125 , .2846828704729192 },
+ { 1.375 , 1.25 , .3763336820249054 },
+ { 1.375 , 1.375 , .4752146669149371 },
+ { 1.375 , 1.5 , .5809359740231859 },
+ { 1.375 , 1.625 , .6931471805599453 },
+ { 1.375 , 1.75 , 0.811531653906724 },
+ { 1.375 , 1.875 , .9358019311087253 },
+ { 1.375 , 2.0 , 1.06569589786406 },
+ { 1.5 , 1.0 , .2846828704729192 },
+ { 1.5 , 1.125 , .3763336820249054 },
+ { 1.5 , 1.25 , .4752146669149371 },
+ { 1.5 , 1.375 , .5809359740231859 },
+ { 1.5 , 1.5 , .6931471805599453 },
+ { 1.5 , 1.625 , 0.811531653906724 },
+ { 1.5 , 1.75 , .9358019311087253 },
+ { 1.5 , 1.875 , 1.06569589786406 },
+ { 1.5 , 2.0 , 1.200973602347074 },
+ { 1.625 , 1.0 , .3763336820249054 },
+ { 1.625 , 1.125 , .4752146669149371 },
+ { 1.625 , 1.25 , .5809359740231859 },
+ { 1.625 , 1.375 , .6931471805599453 },
+ { 1.625 , 1.5 , 0.811531653906724 },
+ { 1.625 , 1.625 , .9358019311087253 },
+ { 1.625 , 1.75 , 1.06569589786406 },
+ { 1.625 , 1.875 , 1.200973602347074 },
+ { 1.625 , 2.0 , 1.341414578068493 },
+ { 1.75 , 1.0 , .4752146669149371 },
+ { 1.75 , 1.125 , .5809359740231859 },
+ { 1.75 , 1.25 , .6931471805599453 },
+ { 1.75 , 1.375 , 0.811531653906724 },
+ { 1.75 , 1.5 , .9358019311087253 },
+ { 1.75 , 1.625 , 1.06569589786406 },
+ { 1.75 , 1.75 , 1.200973602347074 },
+ { 1.75 , 1.875 , 1.341414578068493 },
+ { 1.75 , 2.0 , 1.486815578593417 },
+ { 1.875 , 1.0 , .5809359740231859 },
+ { 1.875 , 1.125 , .6931471805599453 },
+ { 1.875 , 1.25 , 0.811531653906724 },
+ { 1.875 , 1.375 , .9358019311087253 },
+ { 1.875 , 1.5 , 1.06569589786406 },
+ { 1.875 , 1.625 , 1.200973602347074 },
+ { 1.875 , 1.75 , 1.341414578068493 },
+ { 1.875 , 1.875 , 1.486815578593417 },
+ { 1.875 , 2.0 , 1.6369886482725 },
+ { 2.0 , 1.0 , .6931471805599453 },
+ { 2.0 , 1.125 , 0.811531653906724 },
+ { 2.0 , 1.25 , .9358019311087253 },
+ { 2.0 , 1.375 , 1.06569589786406 },
+ { 2.0 , 1.5 , 1.200973602347074 },
+ { 2.0 , 1.625 , 1.341414578068493 },
+ { 2.0 , 1.75 , 1.486815578593417 },
+ { 2.0 , 1.875 , 1.6369886482725 },
+ { 2.0 , 2.0 , 1.791759469228055 },
+ };
+
+ @Test
+ public void testLogGammaSum() {
+ final int ulps = 5;
+ for (int i = 0; i < LOG_GAMMA_SUM_REF.length; i++) {
+ final double[] ref = LOG_GAMMA_SUM_REF[i];
+ final double a = ref[0];
+ final double b = ref[1];
+ final double expected = ref[2];
+ final double actual = Gamma.logGammaSum(a, b);
+ final double tol = ulps * FastMath.ulp(expected);
+ final StringBuilder builder = new StringBuilder();
+ builder.append(a).append(", ").append(b);
+ Assert.assertEquals(builder.toString(), expected, actual, tol);
+ }
+ }
+
private void checkRelativeError(String msg, double expected, double
actual, double tolerance) {
Assert.assertEquals(msg, expected, actual, FastMath.abs(tolerance *
actual));
}
