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    }