View Javadoc

    1   /*
    2    * Copyright 2003-2004 The Apache Software Foundation.
    3    *
    4    * Licensed under the Apache License, Version 2.0 (the "License");
    5    * you may not use this file except in compliance with the License.
    6    * You may obtain a copy of the License at
    7    *
    8    *      http://www.apache.org/licenses/LICENSE-2.0
    9    *
    10   * Unless required by applicable law or agreed to in writing, software
    11   * distributed under the License is distributed on an "AS IS" BASIS,
    12   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
    13   * See the License for the specific language governing permissions and
    14   * limitations under the License.
    15   */
    16  package org.apache.commons.math.special;
    17
    18  import java.io.Serializable;
    19
    20  import org.apache.commons.math.MathException;
    21  import org.apache.commons.math.util.ContinuedFraction;
    22
    23  /***
    24   * This is a utility class that provides computation methods related to the
    25   * Beta family of functions.
    26   *
    27   * @version $Revision: 1.21 $ $Date: 2004/06/23 16:26:17 $
    28   */
    29  public class Beta implements Serializable {
    30      /*** Maximum allowed numerical error. */
    31      private static final double DEFAULT_EPSILON = 10e-9;
    32
    33      /***
    34       * Default constructor.  Prohibit instantiation.
    35       */
    36      private Beta() {
    37          super();
    38      }
    39
    40      /***
    41       * Returns the
    42       * <a href="http://mathworld.wolfram.com/RegularizedBetaFunction.html">
    43       * regularized beta function</a> I(x, a, b).
    44       *
    45       * @param x the value.
    46       * @param a the a parameter.
    47       * @param b the b parameter.
    48       * @return the regularized beta function I(x, a, b)
    49       * @throws MathException if the algorithm fails to converge.
    50       */
    51      public static double regularizedBeta(double x, double a, double b)
    52          throws MathException
    53      {
    54          return regularizedBeta(x, a, b, DEFAULT_EPSILON, Integer.MAX_VALUE);
    55      }
    56
    57      /***
    58       * Returns the
    59       * <a href="http://mathworld.wolfram.com/RegularizedBetaFunction.html">
    60       * regularized beta function</a> I(x, a, b).
    61       *
    62       * @param x the value.
    63       * @param a the a parameter.
    64       * @param b the b parameter.
    65       * @param epsilon When the absolute value of the nth item in the
    66       *                series is less than epsilon the approximation ceases
    67       *                to calculate further elements in the series.
    68       * @return the regularized beta function I(x, a, b)
    69       * @throws MathException if the algorithm fails to converge.
    70       */
    71      public static double regularizedBeta(double x, double a, double b,
    72          double epsilon) throws MathException
    73      {
    74          return regularizedBeta(x, a, b, epsilon, Integer.MAX_VALUE);
    75      }
    76
    77      /***
    78       * Returns the regularized beta function I(x, a, b).
    79       *
    80       * @param x the value.
    81       * @param a the a parameter.
    82       * @param b the b parameter.
    83       * @param maxIterations Maximum number of "iterations" to complete.
    84       * @return the regularized beta function I(x, a, b)
    85       * @throws MathException if the algorithm fails to converge.
    86       */
    87      public static double regularizedBeta(double x, double a, double b,
    88          int maxIterations) throws MathException
    89      {
    90          return regularizedBeta(x, a, b, DEFAULT_EPSILON, maxIterations);
    91      }
    92
    93      /***
    94       * Returns the regularized beta function I(x, a, b).
    95       *
    96       * The implementation of this method is based on:
    97       * <ul>
    98       * <li>
    99       * <a href="http://mathworld.wolfram.com/RegularizedBetaFunction.html">
    100      * Regularized Beta Function</a>.</li>
    101      * <li>
    102      * <a href="http://functions.wolfram.com/06.21.10.0001.01">
    103      * Regularized Beta Function</a>.</li>
    104      * </ul>
    105      *
    106      * @param x the value.
    107      * @param a the a parameter.
    108      * @param b the b parameter.
    109      * @param epsilon When the absolute value of the nth item in the
    110      *                series is less than epsilon the approximation ceases
    111      *                to calculate further elements in the series.
    112      * @param maxIterations Maximum number of "iterations" to complete.
    113      * @return the regularized beta function I(x, a, b)
    114      * @throws MathException if the algorithm fails to converge.
    115      */
    116     public static double regularizedBeta(double x, final double a,
    117         final double b, double epsilon, int maxIterations) throws MathException
    118     {
    119         double ret;
    120
    121         if (Double.isNaN(x) || Double.isNaN(a) || Double.isNaN(b) || (x < 0) ||
    122             (x > 1) || (a <= 0.0) || (b <= 0.0))
    123         {
    124             ret = Double.NaN;
    125         } else if (x > (a + 1.0) / (a + b + 2.0)) {
    126             ret = 1.0 - regularizedBeta(1.0 - x, b, a, epsilon, maxIterations);
    127         } else {
    128             ContinuedFraction fraction = new ContinuedFraction() {
    129                 protected double getB(int n, double x) {
    130                     double ret;
    131                     double m;
    132                     if (n % 2 == 0) { // even
    133                         m = n / 2.0;
    134                         ret = (m * (b - m) * x) /
    135                             ((a + (2 * m) - 1) * (a + (2 * m)));
    136                     } else {
    137                         m = (n - 1.0) / 2.0;
    138                         ret = -((a + m) * (a + b + m) * x) /
    139                                 ((a + (2 * m)) * (a + (2 * m) + 1.0));
    140                     }
    141                     return ret;
    142                 }
    143
    144                 protected double getA(int n, double x) {
    145                     return 1.0;
    146                 }
    147             };
    148             ret = Math.exp((a * Math.log(x)) + (b * Math.log(1.0 - x)) -
    149                 Math.log(a) - logBeta(a, b, epsilon, maxIterations)) *
    150                 1.0 / fraction.evaluate(x, epsilon, maxIterations);
    151         }
    152
    153         return ret;
    154     }
    155
    156     /***
    157      * Returns the natural logarithm of the beta function B(a, b).
    158      *
    159      * @param a the a parameter.
    160      * @param b the b parameter.
    161      * @return log(B(a, b))
    162      */
    163     public static double logBeta(double a, double b) {
    164         return logBeta(a, b, DEFAULT_EPSILON, Integer.MAX_VALUE);
    165     }
    166
    167     /***
    168      * Returns the natural logarithm of the beta function B(a, b).
    169      *
    170      * The implementation of this method is based on:
    171      * <ul>
    172      * <li><a href="http://mathworld.wolfram.com/BetaFunction.html">
    173      * Beta Function</a>, equation (1).</li>
    174      * </ul>
    175      *
    176      * @param a the a parameter.
    177      * @param b the b parameter.
    178      * @param epsilon When the absolute value of the nth item in the
    179      *                series is less than epsilon the approximation ceases
    180      *                to calculate further elements in the series.
    181      * @param maxIterations Maximum number of "iterations" to complete.
    182      * @return log(B(a, b))
    183      */
    184     public static double logBeta(double a, double b, double epsilon,
    185         int maxIterations) {
    186
    187         double ret;
    188
    189         if (Double.isNaN(a) || Double.isNaN(b) || (a <= 0.0) || (b <= 0.0)) {
    190             ret = Double.NaN;
    191         } else {
    192             ret = Gamma.logGamma(a) + Gamma.logGamma(b) -
    193                 Gamma.logGamma(a + b);
    194         }
    195
    196         return ret;
    197     }
    198 }