commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From celes...@apache.org
Subject svn commit: r1413369 - in /commons/proper/math/trunk/src: main/java/org/apache/commons/math3/special/Beta.java test/java/org/apache/commons/math3/special/BetaTest.java
Date Sun, 25 Nov 2012 16:44:40 GMT
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);
+    }
 }



Mime
View raw message