    001    /*
    002     * Licensed to the Apache Software Foundation (ASF) under one or more
    003     * contributor license agreements.  See the NOTICE file distributed with
    004     * this work for additional information regarding copyright ownership.
    005     * The ASF licenses this file to You under the Apache License, Version 2.0
    006     * (the "License"); you may not use this file except in compliance with
    007     * the License.  You may obtain a copy of the License at
    008     *
    009     *      http://www.apache.org/licenses/LICENSE-2.0
    010     *
    011     * Unless required by applicable law or agreed to in writing, software
    012     * distributed under the License is distributed on an "AS IS" BASIS,
    013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
    014     * See the License for the specific language governing permissions and
    015     * limitations under the License.
    016     */
    017    package org.apache.commons.math3.special;
    018
    019    import org.apache.commons.math3.exception.NumberIsTooSmallException;
    020    import org.apache.commons.math3.exception.OutOfRangeException;
    021    import org.apache.commons.math3.util.ContinuedFraction;
    022    import org.apache.commons.math3.util.FastMath;
    023
    024    /**
    025     * <p>
    026     * This is a utility class that provides computation methods related to the
    027     * Beta family of functions.
    028     * </p>
    029     * <p>
    030     * Implementation of {@link #logBeta(double, double)} is based on the
    031     * algorithms described in
    032     * <ul>
    033     * <li><a href="http://dx.doi.org/10.1145/22721.23109">Didonato and Morris
    034     *     (1986)</a>, <em>Computation of the Incomplete Gamma Function Ratios
    035     *     and their Inverse</em>, TOMS 12(4), 377-393,</li>
    036     * <li><a href="http://dx.doi.org/10.1145/131766.131776">Didonato and Morris
    037     *     (1992)</a>, <em>Algorithm 708: Significant Digit Computation of the
    038     *     Incomplete Beta Function Ratios</em>, TOMS 18(3), 360-373,</li>
    039     * </ul>
    040     * and implemented in the
    041     * <a href="http://www.dtic.mil/docs/citations/ADA476840">NSWC Library of Mathematical Functions</a>,
    042     * available
    043     * <a href="http://www.ualberta.ca/CNS/RESEARCH/Software/NumericalNSWC/site.html">here</a>.
    044     * This library is "approved for public release", and the
    045     * <a href="http://www.dtic.mil/dtic/pdf/announcements/CopyrightGuidance.pdf">Copyright guidance</a>
    046     * indicates that unless otherwise stated in the code, all FORTRAN functions in
    047     * this library are license free. Since no such notice appears in the code these
    048     * functions can safely be ported to Commons-Math.
    049     * </p>
    050     *
    051     *
    052     * @version $Id: Beta.java 1420669 2012-12-12 13:40:35Z erans $
    053     */
    054    public class Beta {
    055        /** Maximum allowed numerical error. */
    056        private static final double DEFAULT_EPSILON = 1E-14;
    057
    058        /** The constant value of ½log 2π. */
    059        private static final double HALF_LOG_TWO_PI = .9189385332046727;
    060
    061        /**
    062         * <p>
    063         * The coefficients of the series expansion of the Δ function. This function
    064         * is defined as follows
    065         * </p>
    066         * <center>Δ(x) = log Γ(x) - (x - 0.5) log a + a - 0.5 log 2π,</center>
    067         * <p>
    068         * see equation (23) in Didonato and Morris (1992). The series expansion,
    069         * which applies for x ≥ 10, reads
    070         * </p>
    071         * <pre>
    072         *                 14
    073         *                ====
    074         *             1  \                2 n
    075         *     Δ(x) = ---  >    d  (10 / x)
    076         *             x  /      n
    077         *                ====
    078         *                n = 0
    079         * <pre>
    080         */
    081        private static final double[] DELTA = {
    082            .833333333333333333333333333333E-01,
    083            -.277777777777777777777777752282E-04,
    084            .793650793650793650791732130419E-07,
    085            -.595238095238095232389839236182E-09,
    086            .841750841750832853294451671990E-11,
    087            -.191752691751854612334149171243E-12,
    088            .641025640510325475730918472625E-14,
    089            -.295506514125338232839867823991E-15,
    090            .179643716359402238723287696452E-16,
    091            -.139228964661627791231203060395E-17,
    092            .133802855014020915603275339093E-18,
    093            -.154246009867966094273710216533E-19,
    094            .197701992980957427278370133333E-20,
    095            -.234065664793997056856992426667E-21,
    096            .171348014966398575409015466667E-22
    097        };
    098
    099        /**
    100         * Default constructor.  Prohibit instantiation.
    101         */
    102        private Beta() {}
    103
    104        /**
    105         * Returns the
    106         * <a href="http://mathworld.wolfram.com/RegularizedBetaFunction.html">
    107         * regularized beta function</a> I(x, a, b).
    108         *
    109         * @param x Value.
    110         * @param a Parameter {@code a}.
    111         * @param b Parameter {@code b}.
    112         * @return the regularized beta function I(x, a, b).
    113         * @throws org.apache.commons.math3.exception.MaxCountExceededException
    114         * if the algorithm fails to converge.
    115         */
    116        public static double regularizedBeta(double x, double a, double b) {
    117            return regularizedBeta(x, a, b, DEFAULT_EPSILON, Integer.MAX_VALUE);
    118        }
    119
    120        /**
    121         * Returns the
    122         * <a href="http://mathworld.wolfram.com/RegularizedBetaFunction.html">
    123         * regularized beta function</a> I(x, a, b).
    124         *
    125         * @param x Value.
    126         * @param a Parameter {@code a}.
    127         * @param b Parameter {@code b}.
    128         * @param epsilon When the absolute value of the nth item in the
    129         * series is less than epsilon the approximation ceases to calculate
    130         * further elements in the series.
    131         * @return the regularized beta function I(x, a, b)
    132         * @throws org.apache.commons.math3.exception.MaxCountExceededException
    133         * if the algorithm fails to converge.
    134         */
    135        public static double regularizedBeta(double x,
    136                                             double a, double b,
    137                                             double epsilon) {
    138            return regularizedBeta(x, a, b, epsilon, Integer.MAX_VALUE);
    139        }
    140
    141        /**
    142         * Returns the regularized beta function I(x, a, b).
    143         *
    144         * @param x the value.
    145         * @param a Parameter {@code a}.
    146         * @param b Parameter {@code b}.
    147         * @param maxIterations Maximum number of "iterations" to complete.
    148         * @return the regularized beta function I(x, a, b)
    149         * @throws org.apache.commons.math3.exception.MaxCountExceededException
    150         * if the algorithm fails to converge.
    151         */
    152        public static double regularizedBeta(double x,
    153                                             double a, double b,
    154                                             int maxIterations) {
    155            return regularizedBeta(x, a, b, DEFAULT_EPSILON, maxIterations);
    156        }
    157
    158        /**
    159         * Returns the regularized beta function I(x, a, b).
    160         *
    161         * The implementation of this method is based on:
    162         * <ul>
    163         * <li>
    164         * <a href="http://mathworld.wolfram.com/RegularizedBetaFunction.html">
    165         * Regularized Beta Function</a>.</li>
    166         * <li>
    167         * <a href="http://functions.wolfram.com/06.21.10.0001.01">
    168         * Regularized Beta Function</a>.</li>
    169         * </ul>
    170         *
    171         * @param x the value.
    172         * @param a Parameter {@code a}.
    173         * @param b Parameter {@code b}.
    174         * @param epsilon When the absolute value of the nth item in the
    175         * series is less than epsilon the approximation ceases to calculate
    176         * further elements in the series.
    177         * @param maxIterations Maximum number of "iterations" to complete.
    178         * @return the regularized beta function I(x, a, b)
    179         * @throws org.apache.commons.math3.exception.MaxCountExceededException
    180         * if the algorithm fails to converge.
    181         */
    182        public static double regularizedBeta(double x,
    183                                             final double a, final double b,
    184                                             double epsilon, int maxIterations) {
    185            double ret;
    186
    187            if (Double.isNaN(x) ||
    188                Double.isNaN(a) ||
    189                Double.isNaN(b) ||
    190                x < 0 ||
    191                x > 1 ||
    192                a <= 0.0 ||
    193                b <= 0.0) {
    194                ret = Double.NaN;
    195            } else if (x > (a + 1.0) / (a + b + 2.0)) {
    196                ret = 1.0 - regularizedBeta(1.0 - x, b, a, epsilon, maxIterations);
    197            } else {
    198                ContinuedFraction fraction = new ContinuedFraction() {
    199
    200                    @Override
    201                    protected double getB(int n, double x) {
    202                        double ret;
    203                        double m;
    204                        if (n % 2 == 0) { // even
    205                            m = n / 2.0;
    206                            ret = (m * (b - m) * x) /
    207                                ((a + (2 * m) - 1) * (a + (2 * m)));
    208                        } else {
    209                            m = (n - 1.0) / 2.0;
    210                            ret = -((a + m) * (a + b + m) * x) /
    211                                    ((a + (2 * m)) * (a + (2 * m) + 1.0));
    212                        }
    213                        return ret;
    214                    }
    215
    216                    @Override
    217                    protected double getA(int n, double x) {
    218                        return 1.0;
    219                    }
    220                };
    221                ret = FastMath.exp((a * FastMath.log(x)) + (b * FastMath.log(1.0 - x)) -
    222                    FastMath.log(a) - logBeta(a, b)) *
    223                    1.0 / fraction.evaluate(x, epsilon, maxIterations);
    224            }
    225
    226            return ret;
    227        }
    228
    229        /**
    230         * Returns the natural logarithm of the beta function B(a, b).
    231         *
    232         * The implementation of this method is based on:
    233         * <ul>
    234         * <li><a href="http://mathworld.wolfram.com/BetaFunction.html">
    235         * Beta Function</a>, equation (1).</li>
    236         * </ul>
    237         *
    238         * @param a Parameter {@code a}.
    239         * @param b Parameter {@code b}.
    240         * @param epsilon This parameter is ignored.
    241         * @param maxIterations This parameter is ignored.
    242         * @return log(B(a, b)).
    243         * @deprecated as of version 3.1, this method is deprecated as the
    244         * computation of the beta function is no longer iterative; it will be
    245         * removed in version 4.0. Current implementation of this method
    246         * internally calls {@link #logBeta(double, double)}.
    247         */
    248        @Deprecated
    249        public static double logBeta(double a, double b,
    250                                     double epsilon,
    251                                     int maxIterations) {
    252
    253            return logBeta(a, b);
    254        }
    255
    256
    257        /**
    258         * Returns the value of log Γ(a + b) for 1 ≤ a, b ≤ 2. Based on the
    259         * <em>NSWC Library of Mathematics Subroutines</em> double precision
    260         * implementation, {@code DGSMLN}. In {@link BetaTest#testLogGammaSum()},
    261         * this private method is accessed through reflection.
    262         *
    263         * @param a First argument.
    264         * @param b Second argument.
    265         * @return the value of {@code log(Gamma(a + b))}.
    266         * @throws OutOfRangeException if {@code a} or {@code b} is lower than
    267         * {@code 1.0} or greater than {@code 2.0}.
    268         */
    269        private static double logGammaSum(final double a, final double b)
    270            throws OutOfRangeException {
    271
    272            if ((a < 1.0) || (a > 2.0)) {
    273                throw new OutOfRangeException(a, 1.0, 2.0);
    274            }
    275            if ((b < 1.0) || (b > 2.0)) {
    276                throw new OutOfRangeException(b, 1.0, 2.0);
    277            }
    278
    279            final double x = (a - 1.0) + (b - 1.0);
    280            if (x <= 0.5) {
    281                return Gamma.logGamma1p(1.0 + x);
    282            } else if (x <= 1.5) {
    283                return Gamma.logGamma1p(x) + FastMath.log1p(x);
    284            } else {
    285                return Gamma.logGamma1p(x - 1.0) + FastMath.log(x * (1.0 + x));
    286            }
    287        }
    288
    289        /**
    290         * Returns the value of log[Γ(b) / Γ(a + b)] for a ≥ 0 and b ≥ 10. Based on
    291         * the <em>NSWC Library of Mathematics Subroutines</em> double precision
    292         * implementation, {@code DLGDIV}. In
    293         * {@link BetaTest#testLogGammaMinusLogGammaSum()}, this private method is
    294         * accessed through reflection.
    295         *
    296         * @param a First argument.
    297         * @param b Second argument.
    298         * @return the value of {@code log(Gamma(b) / Gamma(a + b))}.
    299         * @throws NumberIsTooSmallException if {@code a < 0.0} or {@code b < 10.0}.
    300         */
    301        private static double logGammaMinusLogGammaSum(final double a,
    302                                                       final double b)
    303            throws NumberIsTooSmallException {
    304
    305            if (a < 0.0) {
    306                throw new NumberIsTooSmallException(a, 0.0, true);
    307            }
    308            if (b < 10.0) {
    309                throw new NumberIsTooSmallException(b, 10.0, true);
    310            }
    311
    312            /*
    313             * d = a + b - 0.5
    314             */
    315            final double d;
    316            final double w;
    317            if (a <= b) {
    318                d = b + (a - 0.5);
    319                w = deltaMinusDeltaSum(a, b);
    320            } else {
    321                d = a + (b - 0.5);
    322                w = deltaMinusDeltaSum(b, a);
    323            }
    324
    325            final double u = d * FastMath.log1p(a / b);
    326            final double v = a * (FastMath.log(b) - 1.0);
    327
    328            return u <= v ? (w - u) - v : (w - v) - u;
    329        }
    330
    331        /**
    332         * Returns the value of Δ(b) - Δ(a + b), with 0 ≤ a ≤ b and b ≥ 10. Based
    333         * on equations (26), (27) and (28) in Didonato and Morris (1992).
    334         *
    335         * @param a First argument.
    336         * @param b Second argument.
    337         * @return the value of {@code Delta(b) - Delta(a + b)}
    338         * @throws OutOfRangeException if {@code a < 0} or {@code a > b}
    339         * @throws NumberIsTooSmallException if {@code b < 10}
    340         */
    341        private static double deltaMinusDeltaSum(final double a,
    342                                                 final double b)
    343            throws OutOfRangeException, NumberIsTooSmallException {
    344
    345            if ((a < 0) || (a > b)) {
    346                throw new OutOfRangeException(a, 0, b);
    347            }
    348            if (b < 10) {
    349                throw new NumberIsTooSmallException(b, 10, true);
    350            }
    351
    352            final double h = a / b;
    353            final double p = h / (1.0 + h);
    354            final double q = 1.0 / (1.0 + h);
    355            final double q2 = q * q;
    356            /*
    357             * s[i] = 1 + q + ... - q**(2 * i)
    358             */
    359            final double[] s = new double[DELTA.length];
    360            s[0] = 1.0;
    361            for (int i = 1; i < s.length; i++) {
    362                s[i] = 1.0 + (q + q2 * s[i - 1]);
    363            }
    364            /*
    365             * w = Delta(b) - Delta(a + b)
    366             */
    367            final double sqrtT = 10.0 / b;
    368            final double t = sqrtT * sqrtT;
    369            double w = DELTA[DELTA.length - 1] * s[s.length - 1];
    370            for (int i = DELTA.length - 2; i >= 0; i--) {
    371                w = t * w + DELTA[i] * s[i];
    372            }
    373            return w * p / b;
    374        }
    375
    376        /**
    377         * Returns the value of Δ(p) + Δ(q) - Δ(p + q), with p, q ≥ 10. Based on
    378         * the <em>NSWC Library of Mathematics Subroutines</em> double precision
    379         * implementation, {@code DBCORR}. In
    380         * {@link BetaTest#testSumDeltaMinusDeltaSum()}, this private method is
    381         * accessed through reflection.
    382         *
    383         * @param p First argument.
    384         * @param q Second argument.
    385         * @return the value of {@code Delta(p) + Delta(q) - Delta(p + q)}.
    386         * @throws NumberIsTooSmallException if {@code p < 10.0} or {@code q < 10.0}.
    387         */
    388        private static double sumDeltaMinusDeltaSum(final double p,
    389                                                    final double q) {
    390
    391            if (p < 10.0) {
    392                throw new NumberIsTooSmallException(p, 10.0, true);
    393            }
    394            if (q < 10.0) {
    395                throw new NumberIsTooSmallException(q, 10.0, true);
    396            }
    397
    398            final double a = FastMath.min(p, q);
    399            final double b = FastMath.max(p, q);
    400            final double sqrtT = 10.0 / a;
    401            final double t = sqrtT * sqrtT;
    402            double z = DELTA[DELTA.length - 1];
    403            for (int i = DELTA.length - 2; i >= 0; i--) {
    404                z = t * z + DELTA[i];
    405            }
    406            return z / a + deltaMinusDeltaSum(a, b);
    407        }
    408
    409        /**
    410         * Returns the value of log B(p, q) for 0 ≤ x ≤ 1 and p, q > 0. Based on the
    411         * <em>NSWC Library of Mathematics Subroutines</em> implementation,
    412         * {@code DBETLN}.
    413         *
    414         * @param p First argument.
    415         * @param q Second argument.
    416         * @return the value of {@code log(Beta(p, q))}, {@code NaN} if
    417         * {@code p <= 0} or {@code q <= 0}.
    418         */
    419        public static double logBeta(final double p, final double q) {
    420            if (Double.isNaN(p) || Double.isNaN(q) || (p <= 0.0) || (q <= 0.0)) {
    421                return Double.NaN;
    422            }
    423
    424            final double a = FastMath.min(p, q);
    425            final double b = FastMath.max(p, q);
    426            if (a >= 10.0) {
    427                final double w = sumDeltaMinusDeltaSum(a, b);
    428                final double h = a / b;
    429                final double c = h / (1.0 + h);
    430                final double u = -(a - 0.5) * FastMath.log(c);
    431                final double v = b * FastMath.log1p(h);
    432                if (u <= v) {
    433                    return (((-0.5 * FastMath.log(b) + HALF_LOG_TWO_PI) + w) - u) - v;
    434                } else {
    435                    return (((-0.5 * FastMath.log(b) + HALF_LOG_TWO_PI) + w) - v) - u;
    436                }
    437            } else if (a > 2.0) {
    438                if (b > 1000.0) {
    439                    final int n = (int) FastMath.floor(a - 1.0);
    440                    double prod = 1.0;
    441                    double ared = a;
    442                    for (int i = 0; i < n; i++) {
    443                        ared -= 1.0;
    444                        prod *= ared / (1.0 + ared / b);
    445                    }
    446                    return (FastMath.log(prod) - n * FastMath.log(b)) +
    447                            (Gamma.logGamma(ared) +
    448                             logGammaMinusLogGammaSum(ared, b));
    449                } else {
    450                    double prod1 = 1.0;
    451                    double ared = a;
    452                    while (ared > 2.0) {
    453                        ared -= 1.0;
    454                        final double h = ared / b;
    455                        prod1 *= h / (1.0 + h);
    456                    }
    457                    if (b < 10.0) {
    458                        double prod2 = 1.0;
    459                        double bred = b;
    460                        while (bred > 2.0) {
    461                            bred -= 1.0;
    462                            prod2 *= bred / (ared + bred);
    463                        }
    464                        return FastMath.log(prod1) +
    465                               FastMath.log(prod2) +
    466                               (Gamma.logGamma(ared) +
    467                               (Gamma.logGamma(bred) -
    468                                logGammaSum(ared, bred)));
    469                    } else {
    470                        return FastMath.log(prod1) +
    471                               Gamma.logGamma(ared) +
    472                               logGammaMinusLogGammaSum(ared, b);
    473                    }
    474                }
    475            } else if (a >= 1.0) {
    476                if (b > 2.0) {
    477                    if (b < 10.0) {
    478                        double prod = 1.0;
    479                        double bred = b;
    480                        while (bred > 2.0) {
    481                            bred -= 1.0;
    482                            prod *= bred / (a + bred);
    483                        }
    484                        return FastMath.log(prod) +
    485                               (Gamma.logGamma(a) +
    486                                (Gamma.logGamma(bred) -
    487                                 logGammaSum(a, bred)));
    488                    } else {
    489                        return Gamma.logGamma(a) +
    490                               logGammaMinusLogGammaSum(a, b);
    491                    }
    492                } else {
    493                    return Gamma.logGamma(a) +
    494                           Gamma.logGamma(b) -
    495                           logGammaSum(a, b);
    496                }
    497            } else {
    498                if (b >= 10.0) {
    499                    return Gamma.logGamma(a) +
    500                           logGammaMinusLogGammaSum(a, b);
    501                } else {
    502                    // The following command is the original NSWC implementation.
    503                    // return Gamma.logGamma(a) +
    504                    // (Gamma.logGamma(b) - Gamma.logGamma(a + b));
    505                    // The following command turns out to be more accurate.
    506                    return FastMath.log(Gamma.gamma(a) * Gamma.gamma(b) /
    507                                        Gamma.gamma(a + b));
    508                }
    509            }
    510        }
    511    }