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 */
017package org.apache.commons.math3.special;
018
019import org.apache.commons.math3.exception.NumberIsTooSmallException;
020import org.apache.commons.math3.exception.OutOfRangeException;
021import org.apache.commons.math3.util.ContinuedFraction;
022import 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 1546350 2013-11-28 11:41:12Z erans $
053 */
054public 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 ||
193            b <= 0) {
194            ret = Double.NaN;
195        } else if (x > (a + 1) / (2 + b + a) &&
196                   1 - x <= (b + 1) / (2 + b + a)) {
197            ret = 1 - regularizedBeta(1 - x, b, a, epsilon, maxIterations);
198        } else {
199            ContinuedFraction fraction = new ContinuedFraction() {
200
201                @Override
202                protected double getB(int n, double x) {
203                    double ret;
204                    double m;
205                    if (n % 2 == 0) { // even
206                        m = n / 2.0;
207                        ret = (m * (b - m) * x) /
208                            ((a + (2 * m) - 1) * (a + (2 * m)));
209                    } else {
210                        m = (n - 1.0) / 2.0;
211                        ret = -((a + m) * (a + b + m) * x) /
212                                ((a + (2 * m)) * (a + (2 * m) + 1.0));
213                    }
214                    return ret;
215                }
216
217                @Override
218                protected double getA(int n, double x) {
219                    return 1.0;
220                }
221            };
222            ret = FastMath.exp((a * FastMath.log(x)) + (b * FastMath.log1p(-x)) -
223                FastMath.log(a) - logBeta(a, b)) *
224                1.0 / fraction.evaluate(x, epsilon, maxIterations);
225        }
226
227        return ret;
228    }
229
230    /**
231     * Returns the natural logarithm of the beta function B(a, b).
232     *
233     * The implementation of this method is based on:
234     * <ul>
235     * <li><a href="http://mathworld.wolfram.com/BetaFunction.html">
236     * Beta Function</a>, equation (1).</li>
237     * </ul>
238     *
239     * @param a Parameter {@code a}.
240     * @param b Parameter {@code b}.
241     * @param epsilon This parameter is ignored.
242     * @param maxIterations This parameter is ignored.
243     * @return log(B(a, b)).
244     * @deprecated as of version 3.1, this method is deprecated as the
245     * computation of the beta function is no longer iterative; it will be
246     * removed in version 4.0. Current implementation of this method
247     * internally calls {@link #logBeta(double, double)}.
248     */
249    @Deprecated
250    public static double logBeta(double a, double b,
251                                 double epsilon,
252                                 int maxIterations) {
253
254        return logBeta(a, b);
255    }
256
257
258    /**
259     * Returns the value of log ?(a + b) for 1 ? a, b ? 2. Based on the
260     * <em>NSWC Library of Mathematics Subroutines</em> double precision
261     * implementation, {@code DGSMLN}. In {@code BetaTest.testLogGammaSum()},
262     * this private method is accessed through reflection.
263     *
264     * @param a First argument.
265     * @param b Second argument.
266     * @return the value of {@code log(Gamma(a + b))}.
267     * @throws OutOfRangeException if {@code a} or {@code b} is lower than
268     * {@code 1.0} or greater than {@code 2.0}.
269     */
270    private static double logGammaSum(final double a, final double b)
271        throws OutOfRangeException {
272
273        if ((a < 1.0) || (a > 2.0)) {
274            throw new OutOfRangeException(a, 1.0, 2.0);
275        }
276        if ((b < 1.0) || (b > 2.0)) {
277            throw new OutOfRangeException(b, 1.0, 2.0);
278        }
279
280        final double x = (a - 1.0) + (b - 1.0);
281        if (x <= 0.5) {
282            return Gamma.logGamma1p(1.0 + x);
283        } else if (x <= 1.5) {
284            return Gamma.logGamma1p(x) + FastMath.log1p(x);
285        } else {
286            return Gamma.logGamma1p(x - 1.0) + FastMath.log(x * (1.0 + x));
287        }
288    }
289
290    /**
291     * Returns the value of log[?(b) / ?(a + b)] for a ? 0 and b ? 10. Based on
292     * the <em>NSWC Library of Mathematics Subroutines</em> double precision
293     * implementation, {@code DLGDIV}. In
294     * {@code BetaTest.testLogGammaMinusLogGammaSum()}, this private method is
295     * accessed through reflection.
296     *
297     * @param a First argument.
298     * @param b Second argument.
299     * @return the value of {@code log(Gamma(b) / Gamma(a + b))}.
300     * @throws NumberIsTooSmallException if {@code a < 0.0} or {@code b < 10.0}.
301     */
302    private static double logGammaMinusLogGammaSum(final double a,
303                                                   final double b)
304        throws NumberIsTooSmallException {
305
306        if (a < 0.0) {
307            throw new NumberIsTooSmallException(a, 0.0, true);
308        }
309        if (b < 10.0) {
310            throw new NumberIsTooSmallException(b, 10.0, true);
311        }
312
313        /*
314         * d = a + b - 0.5
315         */
316        final double d;
317        final double w;
318        if (a <= b) {
319            d = b + (a - 0.5);
320            w = deltaMinusDeltaSum(a, b);
321        } else {
322            d = a + (b - 0.5);
323            w = deltaMinusDeltaSum(b, a);
324        }
325
326        final double u = d * FastMath.log1p(a / b);
327        final double v = a * (FastMath.log(b) - 1.0);
328
329        return u <= v ? (w - u) - v : (w - v) - u;
330    }
331
332    /**
333     * Returns the value of ?(b) - ?(a + b), with 0 ? a ? b and b ? 10. Based
334     * on equations (26), (27) and (28) in Didonato and Morris (1992).
335     *
336     * @param a First argument.
337     * @param b Second argument.
338     * @return the value of {@code Delta(b) - Delta(a + b)}
339     * @throws OutOfRangeException if {@code a < 0} or {@code a > b}
340     * @throws NumberIsTooSmallException if {@code b < 10}
341     */
342    private static double deltaMinusDeltaSum(final double a,
343                                             final double b)
344        throws OutOfRangeException, NumberIsTooSmallException {
345
346        if ((a < 0) || (a > b)) {
347            throw new OutOfRangeException(a, 0, b);
348        }
349        if (b < 10) {
350            throw new NumberIsTooSmallException(b, 10, true);
351        }
352
353        final double h = a / b;
354        final double p = h / (1.0 + h);
355        final double q = 1.0 / (1.0 + h);
356        final double q2 = q * q;
357        /*
358         * s[i] = 1 + q + ... - q**(2 * i)
359         */
360        final double[] s = new double[DELTA.length];
361        s[0] = 1.0;
362        for (int i = 1; i < s.length; i++) {
363            s[i] = 1.0 + (q + q2 * s[i - 1]);
364        }
365        /*
366         * w = Delta(b) - Delta(a + b)
367         */
368        final double sqrtT = 10.0 / b;
369        final double t = sqrtT * sqrtT;
370        double w = DELTA[DELTA.length - 1] * s[s.length - 1];
371        for (int i = DELTA.length - 2; i >= 0; i--) {
372            w = t * w + DELTA[i] * s[i];
373        }
374        return w * p / b;
375    }
376
377    /**
378     * Returns the value of ?(p) + ?(q) - ?(p + q), with p, q ? 10. Based on
379     * the <em>NSWC Library of Mathematics Subroutines</em> double precision
380     * implementation, {@code DBCORR}. In
381     * {@code BetaTest.testSumDeltaMinusDeltaSum()}, this private method is
382     * accessed through reflection.
383     *
384     * @param p First argument.
385     * @param q Second argument.
386     * @return the value of {@code Delta(p) + Delta(q) - Delta(p + q)}.
387     * @throws NumberIsTooSmallException if {@code p < 10.0} or {@code q < 10.0}.
388     */
389    private static double sumDeltaMinusDeltaSum(final double p,
390                                                final double q) {
391
392        if (p < 10.0) {
393            throw new NumberIsTooSmallException(p, 10.0, true);
394        }
395        if (q < 10.0) {
396            throw new NumberIsTooSmallException(q, 10.0, true);
397        }
398
399        final double a = FastMath.min(p, q);
400        final double b = FastMath.max(p, q);
401        final double sqrtT = 10.0 / a;
402        final double t = sqrtT * sqrtT;
403        double z = DELTA[DELTA.length - 1];
404        for (int i = DELTA.length - 2; i >= 0; i--) {
405            z = t * z + DELTA[i];
406        }
407        return z / a + deltaMinusDeltaSum(a, b);
408    }
409
410    /**
411     * Returns the value of log B(p, q) for 0 ? x ? 1 and p, q > 0. Based on the
412     * <em>NSWC Library of Mathematics Subroutines</em> implementation,
413     * {@code DBETLN}.
414     *
415     * @param p First argument.
416     * @param q Second argument.
417     * @return the value of {@code log(Beta(p, q))}, {@code NaN} if
418     * {@code p <= 0} or {@code q <= 0}.
419     */
420    public static double logBeta(final double p, final double q) {
421        if (Double.isNaN(p) || Double.isNaN(q) || (p <= 0.0) || (q <= 0.0)) {
422            return Double.NaN;
423        }
424
425        final double a = FastMath.min(p, q);
426        final double b = FastMath.max(p, q);
427        if (a >= 10.0) {
428            final double w = sumDeltaMinusDeltaSum(a, b);
429            final double h = a / b;
430            final double c = h / (1.0 + h);
431            final double u = -(a - 0.5) * FastMath.log(c);
432            final double v = b * FastMath.log1p(h);
433            if (u <= v) {
434                return (((-0.5 * FastMath.log(b) + HALF_LOG_TWO_PI) + w) - u) - v;
435            } else {
436                return (((-0.5 * FastMath.log(b) + HALF_LOG_TWO_PI) + w) - v) - u;
437            }
438        } else if (a > 2.0) {
439            if (b > 1000.0) {
440                final int n = (int) FastMath.floor(a - 1.0);
441                double prod = 1.0;
442                double ared = a;
443                for (int i = 0; i < n; i++) {
444                    ared -= 1.0;
445                    prod *= ared / (1.0 + ared / b);
446                }
447                return (FastMath.log(prod) - n * FastMath.log(b)) +
448                        (Gamma.logGamma(ared) +
449                         logGammaMinusLogGammaSum(ared, b));
450            } else {
451                double prod1 = 1.0;
452                double ared = a;
453                while (ared > 2.0) {
454                    ared -= 1.0;
455                    final double h = ared / b;
456                    prod1 *= h / (1.0 + h);
457                }
458                if (b < 10.0) {
459                    double prod2 = 1.0;
460                    double bred = b;
461                    while (bred > 2.0) {
462                        bred -= 1.0;
463                        prod2 *= bred / (ared + bred);
464                    }
465                    return FastMath.log(prod1) +
466                           FastMath.log(prod2) +
467                           (Gamma.logGamma(ared) +
468                           (Gamma.logGamma(bred) -
469                            logGammaSum(ared, bred)));
470                } else {
471                    return FastMath.log(prod1) +
472                           Gamma.logGamma(ared) +
473                           logGammaMinusLogGammaSum(ared, b);
474                }
475            }
476        } else if (a >= 1.0) {
477            if (b > 2.0) {
478                if (b < 10.0) {
479                    double prod = 1.0;
480                    double bred = b;
481                    while (bred > 2.0) {
482                        bred -= 1.0;
483                        prod *= bred / (a + bred);
484                    }
485                    return FastMath.log(prod) +
486                           (Gamma.logGamma(a) +
487                            (Gamma.logGamma(bred) -
488                             logGammaSum(a, bred)));
489                } else {
490                    return Gamma.logGamma(a) +
491                           logGammaMinusLogGammaSum(a, b);
492                }
493            } else {
494                return Gamma.logGamma(a) +
495                       Gamma.logGamma(b) -
496                       logGammaSum(a, b);
497            }
498        } else {
499            if (b >= 10.0) {
500                return Gamma.logGamma(a) +
501                       logGammaMinusLogGammaSum(a, b);
502            } else {
503                // The following command is the original NSWC implementation.
504                // return Gamma.logGamma(a) +
505                // (Gamma.logGamma(b) - Gamma.logGamma(a + b));
506                // The following command turns out to be more accurate.
507                return FastMath.log(Gamma.gamma(a) * Gamma.gamma(b) /
508                                    Gamma.gamma(a + b));
509            }
510        }
511    }
512}