diff options
Diffstat (limited to 'src/main/java/org/apache/commons/math3/special/Beta.java')
| -rw-r--r-- | src/main/java/org/apache/commons/math3/special/Beta.java | 478 |
1 files changed, 478 insertions, 0 deletions
diff --git a/src/main/java/org/apache/commons/math3/special/Beta.java b/src/main/java/org/apache/commons/math3/special/Beta.java new file mode 100644 index 0000000..2f6b6da --- /dev/null +++ b/src/main/java/org/apache/commons/math3/special/Beta.java @@ -0,0 +1,478 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.commons.math3.special; + +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; + +/** + * This is a utility class that provides computation methods related to the Beta family of + * functions. + * + * <p>Implementation of {@link #logBeta(double, double)} is based on the algorithms described in + * + * <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><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, + * </ul> + * + * 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. + */ +public class Beta { + /** Maximum allowed numerical error. */ + private static final double DEFAULT_EPSILON = 1E-14; + + /** The constant value of ½log 2π. */ + private static final double HALF_LOG_TWO_PI = .9189385332046727; + + /** + * <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, + * which applies for x ≥ 10, reads + * </p> + * <pre> + * 14 + * ==== + * 1 \ 2 n + * Δ(x) = --- > d (10 / x) + * x / n + * ==== + * n = 0 + * <pre> + */ + 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() {} + + /** + * Returns the <a href="http://mathworld.wolfram.com/RegularizedBetaFunction.html">regularized + * beta function</a> I(x, a, b). + * + * @param x Value. + * @param a Parameter {@code a}. + * @param b Parameter {@code b}. + * @return the regularized beta function I(x, a, b). + * @throws org.apache.commons.math3.exception.MaxCountExceededException if the algorithm fails + * to converge. + */ + public static double regularizedBeta(double x, double a, double b) { + return regularizedBeta(x, a, b, DEFAULT_EPSILON, Integer.MAX_VALUE); + } + + /** + * Returns the <a href="http://mathworld.wolfram.com/RegularizedBetaFunction.html">regularized + * beta function</a> I(x, a, b). + * + * @param x Value. + * @param a Parameter {@code a}. + * @param b Parameter {@code b}. + * @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. + * @return the regularized beta function I(x, a, b) + * @throws org.apache.commons.math3.exception.MaxCountExceededException if the algorithm fails + * to converge. + */ + public static double regularizedBeta(double x, double a, double b, double epsilon) { + return regularizedBeta(x, a, b, epsilon, Integer.MAX_VALUE); + } + + /** + * Returns the regularized beta function I(x, a, b). + * + * @param x the value. + * @param a Parameter {@code a}. + * @param b Parameter {@code b}. + * @param maxIterations Maximum number of "iterations" to complete. + * @return the regularized beta function I(x, a, b) + * @throws org.apache.commons.math3.exception.MaxCountExceededException if the algorithm fails + * to converge. + */ + public static double regularizedBeta(double x, double a, double b, int maxIterations) { + return regularizedBeta(x, a, b, DEFAULT_EPSILON, maxIterations); + } + + /** + * Returns the regularized beta function I(x, a, b). + * + * <p>The implementation of this method is based on: + * + * <ul> + * <li><a href="http://mathworld.wolfram.com/RegularizedBetaFunction.html">Regularized Beta + * Function</a>. + * <li><a href="http://functions.wolfram.com/06.21.10.0001.01">Regularized Beta Function</a>. + * </ul> + * + * @param x the value. + * @param a Parameter {@code a}. + * @param b Parameter {@code b}. + * @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 beta function I(x, a, b) + * @throws org.apache.commons.math3.exception.MaxCountExceededException if the algorithm fails + * to converge. + */ + public static double regularizedBeta( + double x, final double a, final double b, double epsilon, int maxIterations) { + double ret; + + if (Double.isNaN(x) + || Double.isNaN(a) + || Double.isNaN(b) + || x < 0 + || x > 1 + || a <= 0 + || b <= 0) { + ret = Double.NaN; + } else if (x > (a + 1) / (2 + b + a) && 1 - x <= (b + 1) / (2 + b + a)) { + ret = 1 - regularizedBeta(1 - x, b, a, epsilon, maxIterations); + } else { + ContinuedFraction fraction = + new ContinuedFraction() { + + /** {@inheritDoc} */ + @Override + protected double getB(int n, double x) { + double ret; + double m; + if (n % 2 == 0) { // even + m = n / 2.0; + ret = (m * (b - m) * x) / ((a + (2 * m) - 1) * (a + (2 * m))); + } else { + m = (n - 1.0) / 2.0; + ret = + -((a + m) * (a + b + m) * x) + / ((a + (2 * m)) * (a + (2 * m) + 1.0)); + } + return ret; + } + + /** {@inheritDoc} */ + @Override + protected double getA(int n, double x) { + return 1.0; + } + }; + ret = + FastMath.exp( + (a * FastMath.log(x)) + + (b * FastMath.log1p(-x)) + - FastMath.log(a) + - logBeta(a, b)) + * 1.0 + / fraction.evaluate(x, epsilon, maxIterations); + } + + return ret; + } + + /** + * Returns the natural logarithm of the beta function B(a, b). + * + * <p>The implementation of this method is based on: + * + * <ul> + * <li><a href="http://mathworld.wolfram.com/BetaFunction.html">Beta Function</a>, equation + * (1). + * </ul> + * + * @param a Parameter {@code a}. + * @param b Parameter {@code b}. + * @param epsilon This parameter is ignored. + * @param maxIterations This parameter is ignored. + * @return log(B(a, b)). + * @deprecated as of version 3.1, this method is deprecated as the computation of the beta + * function is no longer iterative; it will be removed in version 4.0. Current + * implementation of this method internally calls {@link #logBeta(double, double)}. + */ + @Deprecated + public static double logBeta(double a, double b, double epsilon, int maxIterations) { + + return logBeta(a, b); + } + + /** + * Returns the value of log Γ(a + b) for 1 ≤ a, b ≤ 2. Based on the <em>NSWC Library of + * Mathematics Subroutines</em> double precision implementation, {@code DGSMLN}. In {@code + * BetaTest.testLogGammaSum()}, this private method is accessed through reflection. + * + * @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}. + */ + private 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 - 1.0) + (b - 1.0); + if (x <= 0.5) { + return Gamma.logGamma1p(1.0 + x); + } else if (x <= 1.5) { + return Gamma.logGamma1p(x) + FastMath.log1p(x); + } else { + return Gamma.logGamma1p(x - 1.0) + FastMath.log(x * (1.0 + x)); + } + } + + /** + * Returns the value of log[Γ(b) / Γ(a + b)] for a ≥ 0 and b ≥ 10. Based on the <em>NSWC Library + * of Mathematics Subroutines</em> double precision implementation, {@code DLGDIV}. In {@code + * BetaTest.testLogGammaMinusLogGammaSum()}, this private method is accessed through reflection. + * + * @param a First argument. + * @param b Second argument. + * @return the value of {@code log(Gamma(b) / Gamma(a + b))}. + * @throws NumberIsTooSmallException if {@code a < 0.0} or {@code b < 10.0}. + */ + private static double logGammaMinusLogGammaSum(final double a, final double b) + throws NumberIsTooSmallException { + + if (a < 0.0) { + throw new NumberIsTooSmallException(a, 0.0, true); + } + if (b < 10.0) { + throw new NumberIsTooSmallException(b, 10.0, true); + } + + /* + * d = a + b - 0.5 + */ + final double d; + final double w; + if (a <= b) { + d = b + (a - 0.5); + w = deltaMinusDeltaSum(a, b); + } else { + d = a + (b - 0.5); + w = deltaMinusDeltaSum(b, a); + } + + final double u = d * FastMath.log1p(a / b); + final double v = a * (FastMath.log(b) - 1.0); + + return u <= v ? (w - u) - v : (w - v) - u; + } + + /** + * Returns the value of Δ(b) - Δ(a + b), with 0 ≤ a ≤ b and b ≥ 10. Based on equations (26), + * (27) and (28) in Didonato and Morris (1992). + * + * @param a First argument. + * @param b Second argument. + * @return the value of {@code Delta(b) - Delta(a + b)} + * @throws OutOfRangeException if {@code a < 0} or {@code a > b} + * @throws NumberIsTooSmallException if {@code b < 10} + */ + private static double deltaMinusDeltaSum(final double a, final double b) + throws OutOfRangeException, NumberIsTooSmallException { + + if ((a < 0) || (a > b)) { + throw new OutOfRangeException(a, 0, b); + } + if (b < 10) { + throw new NumberIsTooSmallException(b, 10, true); + } + + final double h = a / b; + final double p = h / (1.0 + h); + final double q = 1.0 / (1.0 + h); + final double q2 = q * q; + /* + * s[i] = 1 + q + ... - q**(2 * i) + */ + final double[] s = new double[DELTA.length]; + s[0] = 1.0; + for (int i = 1; i < s.length; i++) { + s[i] = 1.0 + (q + q2 * s[i - 1]); + } + /* + * w = Delta(b) - Delta(a + b) + */ + final double sqrtT = 10.0 / b; + final double t = sqrtT * sqrtT; + double w = DELTA[DELTA.length - 1] * s[s.length - 1]; + for (int i = DELTA.length - 2; i >= 0; i--) { + w = t * w + DELTA[i] * s[i]; + } + return w * p / b; + } + + /** + * Returns the value of Δ(p) + Δ(q) - Δ(p + q), with p, q ≥ 10. Based on the <em>NSWC Library of + * Mathematics Subroutines</em> double precision implementation, {@code DBCORR}. In {@code + * BetaTest.testSumDeltaMinusDeltaSum()}, this private method is accessed through reflection. + * + * @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}. + */ + private static double sumDeltaMinusDeltaSum(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 sqrtT = 10.0 / a; + final double t = sqrtT * sqrtT; + double z = DELTA[DELTA.length - 1]; + for (int i = DELTA.length - 2; i >= 0; i--) { + z = t * z + DELTA[i]; + } + return z / a + deltaMinusDeltaSum(a, b); + } + + /** + * Returns the value of log B(p, q) for 0 ≤ x ≤ 1 and p, q > 0. Based on the <em>NSWC Library of + * Mathematics Subroutines</em> implementation, {@code DBETLN}. + * + * @param p First argument. + * @param q Second argument. + * @return the value of {@code log(Beta(p, q))}, {@code NaN} if {@code p <= 0} or {@code q <= + * 0}. + */ + public static double logBeta(final double p, final double q) { + if (Double.isNaN(p) || Double.isNaN(q) || (p <= 0.0) || (q <= 0.0)) { + return Double.NaN; + } + + final double a = FastMath.min(p, q); + final double b = FastMath.max(p, q); + if (a >= 10.0) { + final double w = sumDeltaMinusDeltaSum(a, b); + final double h = a / b; + final double c = h / (1.0 + h); + final double u = -(a - 0.5) * FastMath.log(c); + final double v = b * FastMath.log1p(h); + if (u <= v) { + return (((-0.5 * FastMath.log(b) + HALF_LOG_TWO_PI) + w) - u) - v; + } else { + return (((-0.5 * FastMath.log(b) + HALF_LOG_TWO_PI) + w) - v) - u; + } + } else if (a > 2.0) { + if (b > 1000.0) { + final int n = (int) FastMath.floor(a - 1.0); + double prod = 1.0; + double ared = a; + for (int i = 0; i < n; i++) { + ared -= 1.0; + prod *= ared / (1.0 + ared / b); + } + return (FastMath.log(prod) - n * FastMath.log(b)) + + (Gamma.logGamma(ared) + logGammaMinusLogGammaSum(ared, b)); + } else { + double prod1 = 1.0; + double ared = a; + while (ared > 2.0) { + ared -= 1.0; + final double h = ared / b; + prod1 *= h / (1.0 + h); + } + if (b < 10.0) { + double prod2 = 1.0; + double bred = b; + while (bred > 2.0) { + bred -= 1.0; + prod2 *= bred / (ared + bred); + } + return FastMath.log(prod1) + + FastMath.log(prod2) + + (Gamma.logGamma(ared) + + (Gamma.logGamma(bred) - logGammaSum(ared, bred))); + } else { + return FastMath.log(prod1) + + Gamma.logGamma(ared) + + logGammaMinusLogGammaSum(ared, b); + } + } + } else if (a >= 1.0) { + if (b > 2.0) { + if (b < 10.0) { + double prod = 1.0; + double bred = b; + while (bred > 2.0) { + bred -= 1.0; + prod *= bred / (a + bred); + } + return FastMath.log(prod) + + (Gamma.logGamma(a) + (Gamma.logGamma(bred) - logGammaSum(a, bred))); + } else { + return Gamma.logGamma(a) + logGammaMinusLogGammaSum(a, b); + } + } else { + return Gamma.logGamma(a) + Gamma.logGamma(b) - logGammaSum(a, b); + } + } else { + if (b >= 10.0) { + return Gamma.logGamma(a) + logGammaMinusLogGammaSum(a, b); + } else { + // The following command is the original NSWC implementation. + // return Gamma.logGamma(a) + + // (Gamma.logGamma(b) - Gamma.logGamma(a + b)); + // The following command turns out to be more accurate. + return FastMath.log(Gamma.gamma(a) * Gamma.gamma(b) / Gamma.gamma(a + b)); + } + } + } +} |