Author: celestin
Date: Sun Nov 25 16:44:39 2012
New Revision: 1413369
URL: http://svn.apache.org/viewvc?rev=1413369&view=rev
Log:
MATH-738: implementation of Beta.bcorr(double, double), an auxiliary function
for the computation of Beta.beta(double, double).
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/special/Beta.java
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/special/BetaTest.java
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/special/Beta.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/special/Beta.java?rev=1413369&r1=1413368&r2=1413369&view=diff
==============================================================================
---
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/special/Beta.java
(original)
+++
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/special/Beta.java
Sun Nov 25 16:44:39 2012
@@ -16,12 +16,38 @@
*/
package org.apache.commons.math3.special;
+import org.apache.commons.math3.exception.NumberIsTooSmallException;
import org.apache.commons.math3.util.ContinuedFraction;
import org.apache.commons.math3.util.FastMath;
/**
+ * <p>
* This is a utility class that provides computation methods related to the
* Beta family of functions.
+ * </p>
+ * <p>
+ * Implementation of {@link #bcorr(double, double)} is based on the algorithms
+ * described in
+ * </p>
+ * <ul>
+ * <li><a href="http://dx.doi.org/10.1145/22721.23109";>Didonato and Morris
+ * (1986)</a>, <em>Computation of the Incomplete Gamma Function Ratios and
+ * their Inverse</em>, TOMS 12(4), 377-393,</li>
+ * <li><a href="http://dx.doi.org/10.1145/131766.131776";>Didonato and Morris
+ * (1992)</a>, <em>Algorithm 708: Significant Digit Computation of the
+ * Incomplete Beta Function Ratios</em>, TOMS 18(3), 360-373,</li>
+ * </ul>
+ * <p>
+ * and implemented in the
+ * <a href="http://www.dtic.mil/docs/citations/ADA476840";>NSWC Library of
Mathematical Functions</a>,
+ * available
+ * <a
href="http://www.ualberta.ca/CNS/RESEARCH/Software/NumericalNSWC/site.html";>here</a>.
+ * This library is "approved for public release", and the
+ * <a
href="http://www.dtic.mil/dtic/pdf/announcements/CopyrightGuidance.pdf";>Copyright
guidance</a>
+ * indicates that unless otherwise stated in the code, all FORTRAN functions in
+ * this library are license free. Since no such notice appears in the code
these
+ * functions can safely be ported to Commons-Math.
+ * </p>
*
* @version $Id$
*/
@@ -30,6 +56,47 @@ public class Beta {
private static final double DEFAULT_EPSILON = 10e-15;
/**
+ * <p>
+ * The coefficients of the series expansion of the Î function. This
+ * function is defined as follows
+ * </p>
+ * <center>Î(x) = log Î(x) - (x - 0.5) log a + a - 0.5 log 2Ï,</center>
+ * <p>
+ * see equation (23) in Didonato and Morris (1992). The series expansion
+ * reads
+ * </p>
+ * <pre>
+ * n
+ * ====
+ * 1 \ i
+ * Î(x) = --- > DELTA t
+ * x / i
+ * ====
+ * i = 0
+ * </pre>
+ * <p>
+ * where {@code t = (10 / x)^2}. This series applies for {@code x >= 10.0}.
+ * </p>
+ */
+ private static final double[] DELTA = {
+ .833333333333333333333333333333E-01,
+ -.277777777777777777777777752282E-04,
+ .793650793650793650791732130419E-07,
+ -.595238095238095232389839236182E-09,
+ .841750841750832853294451671990E-11,
+ -.191752691751854612334149171243E-12,
+ .641025640510325475730918472625E-14,
+ -.295506514125338232839867823991E-15,
+ .179643716359402238723287696452E-16,
+ -.139228964661627791231203060395E-17,
+ .133802855014020915603275339093E-18,
+ -.154246009867966094273710216533E-19,
+ .197701992980957427278370133333E-20,
+ -.234065664793997056856992426667E-21,
+ .171348014966398575409015466667E-22
+ };
+
+ /**
* Default constructor. Prohibit instantiation.
*/
private Beta() {}
@@ -204,4 +271,59 @@ public class Beta {
return ret;
}
+
+ /**
+ * Returns the value of Î(p) + Î(q) - Î(p + q), with p, q ⥠10. Based
on
+ * the <em>NSWC Library of Mathematics Subroutines</em> implementation,
+ * {@code BCORR}.
+ *
+ * @param p First argument.
+ * @param q Second argument.
+ * @return the value of {@code Delta(p) + Delta(q) - Delta(p + q)}.
+ * @throws NumberIsTooSmallException if {@code p < 10.0} or {@code q <
10.0}.
+ */
+ public static final double bcorr(final double p, final double q) {
+
+ if (p < 10.0) {
+ throw new NumberIsTooSmallException(p, 10.0, true);
+ }
+ if (q < 10.0) {
+ throw new NumberIsTooSmallException(q, 10.0, true);
+ }
+
+ final double a = FastMath.min(p, q);
+ final double b = FastMath.max(p, q);
+ final double h = a / b;
+ final double c = h / (1.0 + h);
+ final double x = 1.0 / (1.0 + h);
+ final double x2 = x * x;
+ /*
+ * Compute s[i] = (1 - x**(2 * i + 1)) / (1 - x)
+ */
+ final double[] s = new double[DELTA.length];
+ s[0] = 1.0;
+ for (int i = 1; i < s.length; i++) {
+ s[i] = 1.0 + (x + x2 * s[i - 1]);
+ }
+ /*
+ * Set w = Delta(b) - Delta(a + b)
+ */
+ double tmp = 10.0 / b;
+ final double tb = tmp * tmp;
+ double w = DELTA[DELTA.length - 1] * s[s.length - 1];
+ for (int i = DELTA.length - 2; i >= 0; i--) {
+ w = tb * w + DELTA[i] * s[i];
+ }
+ w *= c / b;
+ /*
+ * Compute Delta(a) + w
+ */
+ tmp = 10.0 / a;
+ final double ta = tmp * tmp;
+ double z = DELTA[DELTA.length - 1];
+ for (int i = DELTA.length - 2; i >= 0; i--) {
+ z = ta * z + DELTA[i];
+ }
+ return z / a + w;
+ }
}
Modified:
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/special/BetaTest.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/special/BetaTest.java?rev=1413369&r1=1413368&r2=1413369&view=diff
==============================================================================
---
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/special/BetaTest.java
(original)
+++
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/special/BetaTest.java
Sun Nov 25 16:44:39 2012
@@ -17,7 +17,10 @@
package org.apache.commons.math3.special;
import org.apache.commons.math3.TestUtils;
+import org.apache.commons.math3.exception.NumberIsTooSmallException;
+import org.apache.commons.math3.util.FastMath;
+import org.junit.Assert;
import org.junit.Test;
/**
@@ -119,4 +122,157 @@ public class BetaTest {
public void testLogBetaPositivePositive() {
testLogBeta(-0.693147180559945, 1.0, 2.0);
}
+
+ private static final double[][] BCORR_REF = {
+ { 10.0 , 10.0 , .01249480717472882 },
+ { 10.0 , 11.0 , .01193628470267385 },
+ { 10.0 , 12.0 , .01148578547212797 },
+ { 10.0 , 13.0 , .01111659739668398 },
+ { 10.0 , 14.0 , .01080991216314295 },
+ { 10.0 , 15.0 , .01055214134859758 },
+ { 10.0 , 16.0 , .01033324912491747 },
+ { 10.0 , 17.0 , .01014568069918883 },
+ { 10.0 , 18.0 , .009983653199146491 },
+ { 10.0 , 19.0 , .009842674320242729 },
+ { 10.0 , 20.0 , 0.0097192081956071 },
+ { 11.0 , 10.0 , .01193628470267385 },
+ { 11.0 , 11.0 , .01135973290745925 },
+ { 11.0 , 12.0 , .01089355537047828 },
+ { 11.0 , 13.0 , .01051064829297728 },
+ { 11.0 , 14.0 , 0.0101918899639826 },
+ { 11.0 , 15.0 , .009923438811859604 },
+ { 11.0 , 16.0 , .009695052724952705 },
+ { 11.0 , 17.0 , 0.00949900745283617 },
+ { 11.0 , 18.0 , .009329379874933402 },
+ { 11.0 , 19.0 , 0.00918156080743147 },
+ { 11.0 , 20.0 , 0.00905191635141762 },
+ { 12.0 , 10.0 , .01148578547212797 },
+ { 12.0 , 11.0 , .01089355537047828 },
+ { 12.0 , 12.0 , .01041365883144029 },
+ { 12.0 , 13.0 , .01001867865848564 },
+ { 12.0 , 14.0 , 0.00968923999191334 },
+ { 12.0 , 15.0 , .009411294976563555 },
+ { 12.0 , 16.0 , .009174432043268762 },
+ { 12.0 , 17.0 , .008970786693291802 },
+ { 12.0 , 18.0 , .008794318926790865 },
+ { 12.0 , 19.0 , .008640321527910711 },
+ { 12.0 , 20.0 , .008505077879954796 },
+ { 13.0 , 10.0 , .01111659739668398 },
+ { 13.0 , 11.0 , .01051064829297728 },
+ { 13.0 , 12.0 , .01001867865848564 },
+ { 13.0 , 13.0 , .009613018147953376 },
+ { 13.0 , 14.0 , .009274085618154277 },
+ { 13.0 , 15.0 , 0.0089876637564166 },
+ { 13.0 , 16.0 , .008743200745261382 },
+ { 13.0 , 17.0 , .008532715206686251 },
+ { 13.0 , 18.0 , .008350069108807093 },
+ { 13.0 , 19.0 , .008190472517984874 },
+ { 13.0 , 20.0 , .008050138630244345 },
+ { 14.0 , 10.0 , .01080991216314295 },
+ { 14.0 , 11.0 , 0.0101918899639826 },
+ { 14.0 , 12.0 , 0.00968923999191334 },
+ { 14.0 , 13.0 , .009274085618154277 },
+ { 14.0 , 14.0 , .008926676241967286 },
+ { 14.0 , 15.0 , .008632654302369184 },
+ { 14.0 , 16.0 , .008381351102615795 },
+ { 14.0 , 17.0 , .008164687232662443 },
+ { 14.0 , 18.0 , .007976441942841219 },
+ { 14.0 , 19.0 , .007811755112234388 },
+ { 14.0 , 20.0 , .007666780069317652 },
+ { 15.0 , 10.0 , .01055214134859758 },
+ { 15.0 , 11.0 , .009923438811859604 },
+ { 15.0 , 12.0 , .009411294976563555 },
+ { 15.0 , 13.0 , 0.0089876637564166 },
+ { 15.0 , 14.0 , .008632654302369184 },
+ { 15.0 , 15.0 , 0.00833179217417291 },
+ { 15.0 , 16.0 , .008074310643041299 },
+ { 15.0 , 17.0 , .007852047581145882 },
+ { 15.0 , 18.0 , .007658712051540045 },
+ { 15.0 , 19.0 , .007489384065757007 },
+ { 15.0 , 20.0 , .007340165635725612 },
+ { 16.0 , 10.0 , .01033324912491747 },
+ { 16.0 , 11.0 , .009695052724952705 },
+ { 16.0 , 12.0 , .009174432043268762 },
+ { 16.0 , 13.0 , .008743200745261382 },
+ { 16.0 , 14.0 , .008381351102615795 },
+ { 16.0 , 15.0 , .008074310643041299 },
+ { 16.0 , 16.0 , .007811229919967624 },
+ { 16.0 , 17.0 , .007583876618287594 },
+ { 16.0 , 18.0 , .007385899933505551 },
+ { 16.0 , 19.0 , .007212328560607852 },
+ { 16.0 , 20.0 , .007059220321091879 },
+ { 17.0 , 10.0 , .01014568069918883 },
+ { 17.0 , 11.0 , 0.00949900745283617 },
+ { 17.0 , 12.0 , .008970786693291802 },
+ { 17.0 , 13.0 , .008532715206686251 },
+ { 17.0 , 14.0 , .008164687232662443 },
+ { 17.0 , 15.0 , .007852047581145882 },
+ { 17.0 , 16.0 , .007583876618287594 },
+ { 17.0 , 17.0 , .007351882161431358 },
+ { 17.0 , 18.0 , .007149662089534654 },
+ { 17.0 , 19.0 , .006972200907152378 },
+ { 17.0 , 20.0 , .006815518216094137 },
+ { 18.0 , 10.0 , .009983653199146491 },
+ { 18.0 , 11.0 , .009329379874933402 },
+ { 18.0 , 12.0 , .008794318926790865 },
+ { 18.0 , 13.0 , .008350069108807093 },
+ { 18.0 , 14.0 , .007976441942841219 },
+ { 18.0 , 15.0 , .007658712051540045 },
+ { 18.0 , 16.0 , .007385899933505551 },
+ { 18.0 , 17.0 , .007149662089534654 },
+ { 18.0 , 18.0 , .006943552208153373 },
+ { 18.0 , 19.0 , .006762516574228829 },
+ { 18.0 , 20.0 , .006602541598043117 },
+ { 19.0 , 10.0 , .009842674320242729 },
+ { 19.0 , 11.0 , 0.00918156080743147 },
+ { 19.0 , 12.0 , .008640321527910711 },
+ { 19.0 , 13.0 , .008190472517984874 },
+ { 19.0 , 14.0 , .007811755112234388 },
+ { 19.0 , 15.0 , .007489384065757007 },
+ { 19.0 , 16.0 , .007212328560607852 },
+ { 19.0 , 17.0 , .006972200907152378 },
+ { 19.0 , 18.0 , .006762516574228829 },
+ { 19.0 , 19.0 , .006578188655176814 },
+ { 19.0 , 20.0 , .006415174623476747 },
+ { 20.0 , 10.0 , 0.0097192081956071 },
+ { 20.0 , 11.0 , 0.00905191635141762 },
+ { 20.0 , 12.0 , .008505077879954796 },
+ { 20.0 , 13.0 , .008050138630244345 },
+ { 20.0 , 14.0 , .007666780069317652 },
+ { 20.0 , 15.0 , .007340165635725612 },
+ { 20.0 , 16.0 , .007059220321091879 },
+ { 20.0 , 17.0 , .006815518216094137 },
+ { 20.0 , 18.0 , .006602541598043117 },
+ { 20.0 , 19.0 , .006415174623476747 },
+ { 20.0 , 20.0 , .006249349445691423 },
+ };
+
+ @Test
+ public void testBcorr() {
+
+ final int ulps = 3;
+ for (int i = 0; i < BCORR_REF.length; i++) {
+ final double[] ref = BCORR_REF[i];
+ final double a = ref[0];
+ final double b = ref[1];
+ final double expected = ref[2];
+ final double actual = Beta.bcorr(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);
+ }
+ }
+
+ @Test(expected = NumberIsTooSmallException.class)
+ public void testBcorrPrecondition1() {
+
+ Beta.bcorr(9.0, 10.0);
+ }
+
+ @Test(expected = NumberIsTooSmallException.class)
+ public void testBcorrPrecondition2() {
+
+ Beta.bcorr(10.0, 9.0);
+ }
}
