In class Beta, removed auxiliary function bcorr, as it is not ready to be
included in version 3.1 of Commons-Math. git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1414525 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Sebastien Brisard
parent
2786244d09
commit
2fde62aff6
|
|
@ -16,38 +16,12 @@
|
|||
*/
|
||||
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$
|
||||
*/
|
||||
|
|
@ -55,47 +29,6 @@ public class Beta {
|
|||
/** Maximum allowed numerical error. */
|
||||
private static final double DEFAULT_EPSILON = 1E-14;
|
||||
|
||||
/**
|
||||
* <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.
|
||||
*/
|
||||
|
|
@ -271,59 +204,4 @@ 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}.
|
||||
*/
|
||||
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;
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -17,10 +17,6 @@
|
|||
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;
|
||||
|
||||
/**
|
||||
|
|
@ -122,157 +118,4 @@ 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);
|
||||
}
|
||||
}
|
||||
|
|
|
|||
Loading…
Reference in New Issue