    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.math.special;
    018
    019    import org.apache.commons.math.MathException;
    020    import org.apache.commons.math.util.ContinuedFraction;
    021    import org.apache.commons.math.util.FastMath;
    022
    023    /**
    024     * This is a utility class that provides computation methods related to the
    025     * Beta family of functions.
    026     *
    027     * @version $Revision: 990655 $ $Date: 2010-08-29 23:49:40 +0200 (dim. 29 ao??t 2010) $
    028     */
    029    public class Beta {
    030
    031        /** Maximum allowed numerical error. */
    032        private static final double DEFAULT_EPSILON = 10e-15;
    033
    034        /**
    035         * Default constructor.  Prohibit instantiation.
    036         */
    037        private Beta() {
    038            super();
    039        }
    040
    041        /**
    042         * Returns the
    043         * <a href="http://mathworld.wolfram.com/RegularizedBetaFunction.html">
    044         * regularized beta function</a> I(x, a, b).
    045         *
    046         * @param x the value.
    047         * @param a the a parameter.
    048         * @param b the b parameter.
    049         * @return the regularized beta function I(x, a, b)
    050         * @throws MathException if the algorithm fails to converge.
    051         */
    052        public static double regularizedBeta(double x, double a, double b)
    053            throws MathException
    054        {
    055            return regularizedBeta(x, a, b, DEFAULT_EPSILON, Integer.MAX_VALUE);
    056        }
    057
    058        /**
    059         * Returns the
    060         * <a href="http://mathworld.wolfram.com/RegularizedBetaFunction.html">
    061         * regularized beta function</a> I(x, a, b).
    062         *
    063         * @param x the value.
    064         * @param a the a parameter.
    065         * @param b the b parameter.
    066         * @param epsilon When the absolute value of the nth item in the
    067         *                series is less than epsilon the approximation ceases
    068         *                to calculate further elements in the series.
    069         * @return the regularized beta function I(x, a, b)
    070         * @throws MathException if the algorithm fails to converge.
    071         */
    072        public static double regularizedBeta(double x, double a, double b,
    073            double epsilon) throws MathException
    074        {
    075            return regularizedBeta(x, a, b, epsilon, Integer.MAX_VALUE);
    076        }
    077
    078        /**
    079         * Returns the regularized beta function I(x, a, b).
    080         *
    081         * @param x the value.
    082         * @param a the a parameter.
    083         * @param b the b parameter.
    084         * @param maxIterations Maximum number of "iterations" to complete.
    085         * @return the regularized beta function I(x, a, b)
    086         * @throws MathException if the algorithm fails to converge.
    087         */
    088        public static double regularizedBeta(double x, double a, double b,
    089            int maxIterations) throws MathException
    090        {
    091            return regularizedBeta(x, a, b, DEFAULT_EPSILON, maxIterations);
    092        }
    093
    094        /**
    095         * Returns the regularized beta function I(x, a, b).
    096         *
    097         * The implementation of this method is based on:
    098         * <ul>
    099         * <li>
    100         * <a href="http://mathworld.wolfram.com/RegularizedBetaFunction.html">
    101         * Regularized Beta Function</a>.</li>
    102         * <li>
    103         * <a href="http://functions.wolfram.com/06.21.10.0001.01">
    104         * Regularized Beta Function</a>.</li>
    105         * </ul>
    106         *
    107         * @param x the value.
    108         * @param a the a parameter.
    109         * @param b the b parameter.
    110         * @param epsilon When the absolute value of the nth item in the
    111         *                series is less than epsilon the approximation ceases
    112         *                to calculate further elements in the series.
    113         * @param maxIterations Maximum number of "iterations" to complete.
    114         * @return the regularized beta function I(x, a, b)
    115         * @throws MathException if the algorithm fails to converge.
    116         */
    117        public static double regularizedBeta(double x, final double a,
    118            final double b, double epsilon, int maxIterations) throws MathException
    119        {
    120            double ret;
    121
    122            if (Double.isNaN(x) || Double.isNaN(a) || Double.isNaN(b) || (x < 0) ||
    123                (x > 1) || (a <= 0.0) || (b <= 0.0))
    124            {
    125                ret = Double.NaN;
    126            } else if (x > (a + 1.0) / (a + b + 2.0)) {
    127                ret = 1.0 - regularizedBeta(1.0 - x, b, a, epsilon, maxIterations);
    128            } else {
    129                ContinuedFraction fraction = new ContinuedFraction() {
    130
    131                    @Override
    132                    protected double getB(int n, double x) {
    133                        double ret;
    134                        double m;
    135                        if (n % 2 == 0) { // even
    136                            m = n / 2.0;
    137                            ret = (m * (b - m) * x) /
    138                                ((a + (2 * m) - 1) * (a + (2 * m)));
    139                        } else {
    140                            m = (n - 1.0) / 2.0;
    141                            ret = -((a + m) * (a + b + m) * x) /
    142                                    ((a + (2 * m)) * (a + (2 * m) + 1.0));
    143                        }
    144                        return ret;
    145                    }
    146
    147                    @Override
    148                    protected double getA(int n, double x) {
    149                        return 1.0;
    150                    }
    151                };
    152                ret = FastMath.exp((a * FastMath.log(x)) + (b * FastMath.log(1.0 - x)) -
    153                    FastMath.log(a) - logBeta(a, b, epsilon, maxIterations)) *
    154                    1.0 / fraction.evaluate(x, epsilon, maxIterations);
    155            }
    156
    157            return ret;
    158        }
    159
    160        /**
    161         * Returns the natural logarithm of the beta function B(a, b).
    162         *
    163         * @param a the a parameter.
    164         * @param b the b parameter.
    165         * @return log(B(a, b))
    166         */
    167        public static double logBeta(double a, double b) {
    168            return logBeta(a, b, DEFAULT_EPSILON, Integer.MAX_VALUE);
    169        }
    170
    171        /**
    172         * Returns the natural logarithm of the beta function B(a, b).
    173         *
    174         * The implementation of this method is based on:
    175         * <ul>
    176         * <li><a href="http://mathworld.wolfram.com/BetaFunction.html">
    177         * Beta Function</a>, equation (1).</li>
    178         * </ul>
    179         *
    180         * @param a the a parameter.
    181         * @param b the b parameter.
    182         * @param epsilon When the absolute value of the nth item in the
    183         *                series is less than epsilon the approximation ceases
    184         *                to calculate further elements in the series.
    185         * @param maxIterations Maximum number of "iterations" to complete.
    186         * @return log(B(a, b))
    187         */
    188        public static double logBeta(double a, double b, double epsilon,
    189            int maxIterations) {
    190
    191            double ret;
    192
    193            if (Double.isNaN(a) || Double.isNaN(b) || (a <= 0.0) || (b <= 0.0)) {
    194                ret = Double.NaN;
    195            } else {
    196                ret = Gamma.logGamma(a) + Gamma.logGamma(b) -
    197                    Gamma.logGamma(a + b);
    198            }
    199
    200            return ret;
    201        }
    202    }