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 */
053public class Beta {
054    /** Maximum allowed numerical error. */
055    private static final double DEFAULT_EPSILON = 1E-14;
056
057    /** The constant value of ½log 2π. */
058    private static final double HALF_LOG_TWO_PI = .9189385332046727;
059
060    /**
061     * <p>
062     * The coefficients of the series expansion of the Δ function. This function
063     * is defined as follows
064     * </p>
065     * <center>Δ(x) = log Γ(x) - (x - 0.5) log a + a - 0.5 log 2π,</center>
066     * <p>
067     * see equation (23) in Didonato and Morris (1992). The series expansion,
068     * which applies for x ≥ 10, reads
069     * </p>
070     * <pre>
071     *                 14
072     *                ====
073     *             1  \                2 n
074     *     Δ(x) = ---  >    d  (10 / x)
075     *             x  /      n
076     *                ====
077     *                n = 0
078     * <pre>
079     */
080    private static final double[] DELTA = {
081        .833333333333333333333333333333E-01,
082        -.277777777777777777777777752282E-04,
083        .793650793650793650791732130419E-07,
084        -.595238095238095232389839236182E-09,
085        .841750841750832853294451671990E-11,
086        -.191752691751854612334149171243E-12,
087        .641025640510325475730918472625E-14,
088        -.295506514125338232839867823991E-15,
089        .179643716359402238723287696452E-16,
090        -.139228964661627791231203060395E-17,
091        .133802855014020915603275339093E-18,
092        -.154246009867966094273710216533E-19,
093        .197701992980957427278370133333E-20,
094        -.234065664793997056856992426667E-21,
095        .171348014966398575409015466667E-22
096    };
097
098    /**
099     * Default constructor.  Prohibit instantiation.
100     */
101    private Beta() {}
102
103    /**
104     * Returns the
105     * <a href="http://mathworld.wolfram.com/RegularizedBetaFunction.html">
106     * regularized beta function</a> I(x, a, b).
107     *
108     * @param x Value.
109     * @param a Parameter {@code a}.
110     * @param b Parameter {@code b}.
111     * @return the regularized beta function I(x, a, b).
112     * @throws org.apache.commons.math3.exception.MaxCountExceededException
113     * if the algorithm fails to converge.
114     */
115    public static double regularizedBeta(double x, double a, double b) {
116        return regularizedBeta(x, a, b, DEFAULT_EPSILON, Integer.MAX_VALUE);
117    }
118
119    /**
120     * Returns the
121     * <a href="http://mathworld.wolfram.com/RegularizedBetaFunction.html">
122     * regularized beta function</a> I(x, a, b).
123     *
124     * @param x Value.
125     * @param a Parameter {@code a}.
126     * @param b Parameter {@code b}.
127     * @param epsilon When the absolute value of the nth item in the
128     * series is less than epsilon the approximation ceases to calculate
129     * further elements in the series.
130     * @return the regularized beta function I(x, a, b)
131     * @throws org.apache.commons.math3.exception.MaxCountExceededException
132     * if the algorithm fails to converge.
133     */
134    public static double regularizedBeta(double x,
135                                         double a, double b,
136                                         double epsilon) {
137        return regularizedBeta(x, a, b, epsilon, Integer.MAX_VALUE);
138    }
139
140    /**
141     * Returns the regularized beta function I(x, a, b).
142     *
143     * @param x the value.
144     * @param a Parameter {@code a}.
145     * @param b Parameter {@code b}.
146     * @param maxIterations Maximum number of "iterations" to complete.
147     * @return the regularized beta function I(x, a, b)
148     * @throws org.apache.commons.math3.exception.MaxCountExceededException
149     * if the algorithm fails to converge.
150     */
151    public static double regularizedBeta(double x,
152                                         double a, double b,
153                                         int maxIterations) {
154        return regularizedBeta(x, a, b, DEFAULT_EPSILON, maxIterations);
155    }
156
157    /**
158     * Returns the regularized beta function I(x, a, b).
159     *
160     * The implementation of this method is based on:
161     * <ul>
162     * <li>
163     * <a href="http://mathworld.wolfram.com/RegularizedBetaFunction.html">
164     * Regularized Beta Function</a>.</li>
165     * <li>
166     * <a href="http://functions.wolfram.com/06.21.10.0001.01">
167     * Regularized Beta Function</a>.</li>
168     * </ul>
169     *
170     * @param x the value.
171     * @param a Parameter {@code a}.
172     * @param b Parameter {@code b}.
173     * @param epsilon When the absolute value of the nth item in the
174     * series is less than epsilon the approximation ceases to calculate
175     * further elements in the series.
176     * @param maxIterations Maximum number of "iterations" to complete.
177     * @return the regularized beta function I(x, a, b)
178     * @throws org.apache.commons.math3.exception.MaxCountExceededException
179     * if the algorithm fails to converge.
180     */
181    public static double regularizedBeta(double x,
182                                         final double a, final double b,
183                                         double epsilon, int maxIterations) {
184        double ret;
185
186        if (Double.isNaN(x) ||
187            Double.isNaN(a) ||
188            Double.isNaN(b) ||
189            x < 0 ||
190            x > 1 ||
191            a <= 0 ||
192            b <= 0) {
193            ret = Double.NaN;
194        } else if (x > (a + 1) / (2 + b + a) &&
195                   1 - x <= (b + 1) / (2 + b + a)) {
196            ret = 1 - regularizedBeta(1 - x, b, a, epsilon, maxIterations);
197        } else {
198            ContinuedFraction fraction = new ContinuedFraction() {
199
200                /** {@inheritDoc} */
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                /** {@inheritDoc} */
218                @Override
219                protected double getA(int n, double x) {
220                    return 1.0;
221                }
222            };
223            ret = FastMath.exp((a * FastMath.log(x)) + (b * FastMath.log1p(-x)) -
224                FastMath.log(a) - logBeta(a, b)) *
225                1.0 / fraction.evaluate(x, epsilon, maxIterations);
226        }
227
228        return ret;
229    }
230
231    /**
232     * Returns the natural logarithm of the beta function B(a, b).
233     *
234     * The implementation of this method is based on:
235     * <ul>
236     * <li><a href="http://mathworld.wolfram.com/BetaFunction.html">
237     * Beta Function</a>, equation (1).</li>
238     * </ul>
239     *
240     * @param a Parameter {@code a}.
241     * @param b Parameter {@code b}.
242     * @param epsilon This parameter is ignored.
243     * @param maxIterations This parameter is ignored.
244     * @return log(B(a, b)).
245     * @deprecated as of version 3.1, this method is deprecated as the
246     * computation of the beta function is no longer iterative; it will be
247     * removed in version 4.0. Current implementation of this method
248     * internally calls {@link #logBeta(double, double)}.
249     */
250    @Deprecated
251    public static double logBeta(double a, double b,
252                                 double epsilon,
253                                 int maxIterations) {
254
255        return logBeta(a, b);
256    }
257
258
259    /**
260     * Returns the value of log Γ(a + b) for 1 ≤ a, b ≤ 2. Based on the
261     * <em>NSWC Library of Mathematics Subroutines</em> double precision
262     * implementation, {@code DGSMLN}. In {@code BetaTest.testLogGammaSum()},
263     * this private method is accessed through reflection.
264     *
265     * @param a First argument.
266     * @param b Second argument.
267     * @return the value of {@code log(Gamma(a + b))}.
268     * @throws OutOfRangeException if {@code a} or {@code b} is lower than
269     * {@code 1.0} or greater than {@code 2.0}.
270     */
271    private static double logGammaSum(final double a, final double b)
272        throws OutOfRangeException {
273
274        if ((a < 1.0) || (a > 2.0)) {
275            throw new OutOfRangeException(a, 1.0, 2.0);
276        }
277        if ((b < 1.0) || (b > 2.0)) {
278            throw new OutOfRangeException(b, 1.0, 2.0);
279        }
280
281        final double x = (a - 1.0) + (b - 1.0);
282        if (x <= 0.5) {
283            return Gamma.logGamma1p(1.0 + x);
284        } else if (x <= 1.5) {
285            return Gamma.logGamma1p(x) + FastMath.log1p(x);
286        } else {
287            return Gamma.logGamma1p(x - 1.0) + FastMath.log(x * (1.0 + x));
288        }
289    }
290
291    /**
292     * Returns the value of log[Γ(b) / Γ(a + b)] for a ≥ 0 and b ≥ 10. Based on
293     * the <em>NSWC Library of Mathematics Subroutines</em> double precision
294     * implementation, {@code DLGDIV}. In
295     * {@code BetaTest.testLogGammaMinusLogGammaSum()}, this private method is
296     * accessed through reflection.
297     *
298     * @param a First argument.
299     * @param b Second argument.
300     * @return the value of {@code log(Gamma(b) / Gamma(a + b))}.
301     * @throws NumberIsTooSmallException if {@code a < 0.0} or {@code b < 10.0}.
302     */
303    private static double logGammaMinusLogGammaSum(final double a,
304                                                   final double b)
305        throws NumberIsTooSmallException {
306
307        if (a < 0.0) {
308            throw new NumberIsTooSmallException(a, 0.0, true);
309        }
310        if (b < 10.0) {
311            throw new NumberIsTooSmallException(b, 10.0, true);
312        }
313
314        /*
315         * d = a + b - 0.5
316         */
317        final double d;
318        final double w;
319        if (a <= b) {
320            d = b + (a - 0.5);
321            w = deltaMinusDeltaSum(a, b);
322        } else {
323            d = a + (b - 0.5);
324            w = deltaMinusDeltaSum(b, a);
325        }
326
327        final double u = d * FastMath.log1p(a / b);
328        final double v = a * (FastMath.log(b) - 1.0);
329
330        return u <= v ? (w - u) - v : (w - v) - u;
331    }
332
333    /**
334     * Returns the value of Δ(b) - Δ(a + b), with 0 ≤ a ≤ b and b ≥ 10. Based
335     * on equations (26), (27) and (28) in Didonato and Morris (1992).
336     *
337     * @param a First argument.
338     * @param b Second argument.
339     * @return the value of {@code Delta(b) - Delta(a + b)}
340     * @throws OutOfRangeException if {@code a < 0} or {@code a > b}
341     * @throws NumberIsTooSmallException if {@code b < 10}
342     */
343    private static double deltaMinusDeltaSum(final double a,
344                                             final double b)
345        throws OutOfRangeException, NumberIsTooSmallException {
346
347        if ((a < 0) || (a > b)) {
348            throw new OutOfRangeException(a, 0, b);
349        }
350        if (b < 10) {
351            throw new NumberIsTooSmallException(b, 10, true);
352        }
353
354        final double h = a / b;
355        final double p = h / (1.0 + h);
356        final double q = 1.0 / (1.0 + h);
357        final double q2 = q * q;
358        /*
359         * s[i] = 1 + q + ... - q**(2 * i)
360         */
361        final double[] s = new double[DELTA.length];
362        s[0] = 1.0;
363        for (int i = 1; i < s.length; i++) {
364            s[i] = 1.0 + (q + q2 * s[i - 1]);
365        }
366        /*
367         * w = Delta(b) - Delta(a + b)
368         */
369        final double sqrtT = 10.0 / b;
370        final double t = sqrtT * sqrtT;
371        double w = DELTA[DELTA.length - 1] * s[s.length - 1];
372        for (int i = DELTA.length - 2; i >= 0; i--) {
373            w = t * w + DELTA[i] * s[i];
374        }
375        return w * p / b;
376    }
377
378    /**
379     * Returns the value of Δ(p) + Δ(q) - Δ(p + q), with p, q ≥ 10. Based on
380     * the <em>NSWC Library of Mathematics Subroutines</em> double precision
381     * implementation, {@code DBCORR}. In
382     * {@code BetaTest.testSumDeltaMinusDeltaSum()}, this private method is
383     * accessed through reflection.
384     *
385     * @param p First argument.
386     * @param q Second argument.
387     * @return the value of {@code Delta(p) + Delta(q) - Delta(p + q)}.
388     * @throws NumberIsTooSmallException if {@code p < 10.0} or {@code q < 10.0}.
389     */
390    private static double sumDeltaMinusDeltaSum(final double p,
391                                                final double q) {
392
393        if (p < 10.0) {
394            throw new NumberIsTooSmallException(p, 10.0, true);
395        }
396        if (q < 10.0) {
397            throw new NumberIsTooSmallException(q, 10.0, true);
398        }
399
400        final double a = FastMath.min(p, q);
401        final double b = FastMath.max(p, q);
402        final double sqrtT = 10.0 / a;
403        final double t = sqrtT * sqrtT;
404        double z = DELTA[DELTA.length - 1];
405        for (int i = DELTA.length - 2; i >= 0; i--) {
406            z = t * z + DELTA[i];
407        }
408        return z / a + deltaMinusDeltaSum(a, b);
409    }
410
411    /**
412     * Returns the value of log B(p, q) for 0 ≤ x ≤ 1 and p, q > 0. Based on the
413     * <em>NSWC Library of Mathematics Subroutines</em> implementation,
414     * {@code DBETLN}.
415     *
416     * @param p First argument.
417     * @param q Second argument.
418     * @return the value of {@code log(Beta(p, q))}, {@code NaN} if
419     * {@code p <= 0} or {@code q <= 0}.
420     */
421    public static double logBeta(final double p, final double q) {
422        if (Double.isNaN(p) || Double.isNaN(q) || (p <= 0.0) || (q <= 0.0)) {
423            return Double.NaN;
424        }
425
426        final double a = FastMath.min(p, q);
427        final double b = FastMath.max(p, q);
428        if (a >= 10.0) {
429            final double w = sumDeltaMinusDeltaSum(a, b);
430            final double h = a / b;
431            final double c = h / (1.0 + h);
432            final double u = -(a - 0.5) * FastMath.log(c);
433            final double v = b * FastMath.log1p(h);
434            if (u <= v) {
435                return (((-0.5 * FastMath.log(b) + HALF_LOG_TWO_PI) + w) - u) - v;
436            } else {
437                return (((-0.5 * FastMath.log(b) + HALF_LOG_TWO_PI) + w) - v) - u;
438            }
439        } else if (a > 2.0) {
440            if (b > 1000.0) {
441                final int n = (int) FastMath.floor(a - 1.0);
442                double prod = 1.0;
443                double ared = a;
444                for (int i = 0; i < n; i++) {
445                    ared -= 1.0;
446                    prod *= ared / (1.0 + ared / b);
447                }
448                return (FastMath.log(prod) - n * FastMath.log(b)) +
449                        (Gamma.logGamma(ared) +
450                         logGammaMinusLogGammaSum(ared, b));
451            } else {
452                double prod1 = 1.0;
453                double ared = a;
454                while (ared > 2.0) {
455                    ared -= 1.0;
456                    final double h = ared / b;
457                    prod1 *= h / (1.0 + h);
458                }
459                if (b < 10.0) {
460                    double prod2 = 1.0;
461                    double bred = b;
462                    while (bred > 2.0) {
463                        bred -= 1.0;
464                        prod2 *= bred / (ared + bred);
465                    }
466                    return FastMath.log(prod1) +
467                           FastMath.log(prod2) +
468                           (Gamma.logGamma(ared) +
469                           (Gamma.logGamma(bred) -
470                            logGammaSum(ared, bred)));
471                } else {
472                    return FastMath.log(prod1) +
473                           Gamma.logGamma(ared) +
474                           logGammaMinusLogGammaSum(ared, b);
475                }
476            }
477        } else if (a >= 1.0) {
478            if (b > 2.0) {
479                if (b < 10.0) {
480                    double prod = 1.0;
481                    double bred = b;
482                    while (bred > 2.0) {
483                        bred -= 1.0;
484                        prod *= bred / (a + bred);
485                    }
486                    return FastMath.log(prod) +
487                           (Gamma.logGamma(a) +
488                            (Gamma.logGamma(bred) -
489                             logGammaSum(a, bred)));
490                } else {
491                    return Gamma.logGamma(a) +
492                           logGammaMinusLogGammaSum(a, b);
493                }
494            } else {
495                return Gamma.logGamma(a) +
496                       Gamma.logGamma(b) -
497                       logGammaSum(a, b);
498            }
499        } else {
500            if (b >= 10.0) {
501                return Gamma.logGamma(a) +
502                       logGammaMinusLogGammaSum(a, b);
503            } else {
504                // The following command is the original NSWC implementation.
505                // return Gamma.logGamma(a) +
506                // (Gamma.logGamma(b) - Gamma.logGamma(a + b));
507                // The following command turns out to be more accurate.
508                return FastMath.log(Gamma.gamma(a) * Gamma.gamma(b) /
509                                    Gamma.gamma(a + b));
510            }
511        }
512    }
513}