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    }