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.ConvergenceException;
    21  import org.apache.commons.math.MathException;
    22  import org.apache.commons.math.util.ContinuedFraction;
    23
    24  /***
    25   * This is a utility class that provides computation methods related to the
    26   * Gamma family of functions.
    27   *
    28   * @version $Revision: 1.22 $ $Date: 2004/10/08 05:53:18 $
    29   */
    30  public class Gamma implements Serializable {
    31
    32      /*** Maximum allowed numerical error. */
    33      private static final double DEFAULT_EPSILON = 10e-9;
    34
    35      /*** Lanczos coefficients */
    36      private static double[] lanczos =
    37      {
    38          0.99999999999999709182,
    39          57.156235665862923517,
    40          -59.597960355475491248,
    41          14.136097974741747174,
    42          -0.49191381609762019978,
    43          .33994649984811888699e-4,
    44          .46523628927048575665e-4,
    45          -.98374475304879564677e-4,
    46          .15808870322491248884e-3,
    47          -.21026444172410488319e-3,
    48          .21743961811521264320e-3,
    49          -.16431810653676389022e-3,
    50          .84418223983852743293e-4,
    51          -.26190838401581408670e-4,
    52          .36899182659531622704e-5,
    53      };
    54
    55      /*** Avoid repeated computation of log of 2 PI in logGamma */
    56      private static final double HALF_LOG_2_PI = 0.5 * Math.log(2.0 * Math.PI);
    57
    58
    59      /***
    60       * Default constructor.  Prohibit instantiation.
    61       */
    62      private Gamma() {
    63          super();
    64      }
    65
    66      /***
    67       * Returns the natural logarithm of the gamma function &#915;(x).
    68       *
    69       * The implementation of this method is based on:
    70       * <ul>
    71       * <li><a href="http://mathworld.wolfram.com/GammaFunction.html">
    72       * Gamma Function</a>, equation (28).</li>
    73       * <li><a href="http://mathworld.wolfram.com/LanczosApproximation.html">
    74       * Lanczos Approximation</a>, equations (1) through (5).</li>
    75       * <li><a href="http://my.fit.edu/~gabdo/gamma.txt">Paul Godfrey, A note on
    76       * the computation of the convergent Lanczos complex Gamma approximation
    77       * </a></li>
    78       * </ul>
    79       *
    80       * @param x the value.
    81       * @return log(&#915;(x))
    82       */
    83      public static double logGamma(double x) {
    84          double ret;
    85
    86          if (Double.isNaN(x) || (x <= 0.0)) {
    87              ret = Double.NaN;
    88          } else {
    89              double g = 607.0 / 128.0;
    90
    91              double sum = 0.0;
    92              for (int i = lanczos.length - 1; i > 0; --i) {
    93                  sum = sum + (lanczos[i] / (x + i));
    94              }
    95              sum = sum + lanczos[0];
    96
    97              double tmp = x + g + .5;
    98              ret = ((x + .5) * Math.log(tmp)) - tmp +
    99                  HALF_LOG_2_PI + Math.log(sum / x);
    100         }
    101
    102         return ret;
    103     }
    104
    105     /***
    106      * Returns the regularized gamma function P(a, x).
    107      *
    108      * @param a the a parameter.
    109      * @param x the value.
    110      * @return the regularized gamma function P(a, x)
    111      * @throws MathException if the algorithm fails to converge.
    112      */
    113     public static double regularizedGammaP(double a, double x)
    114         throws MathException
    115     {
    116         return regularizedGammaP(a, x, DEFAULT_EPSILON, Integer.MAX_VALUE);
    117     }
    118
    119
    120     /***
    121      * Returns the regularized gamma function P(a, x).
    122      *
    123      * The implementation of this method is based on:
    124      * <ul>
    125      * <li>
    126      * <a href="http://mathworld.wolfram.com/RegularizedGammaFunction.html">
    127      * Regularized Gamma Function</a>, equation (1).</li>
    128      * <li>
    129      * <a href="http://mathworld.wolfram.com/IncompleteGammaFunction.html">
    130      * Incomplete Gamma Function</a>, equation (4).</li>
    131      * <li>
    132      * <a href="http://mathworld.wolfram.com/ConfluentHypergeometricFunctionoftheFirstKind.html">
    133      * Confluent Hypergeometric Function of the First Kind</a>, equation (1).
    134      * </li>
    135      * </ul>
    136      *
    137      * @param a the a parameter.
    138      * @param x the value.
    139      * @param epsilon When the absolute value of the nth item in the
    140      *                series is less than epsilon the approximation ceases
    141      *                to calculate further elements in the series.
    142      * @param maxIterations Maximum number of "iterations" to complete.
    143      * @return the regularized gamma function P(a, x)
    144      * @throws MathException if the algorithm fails to converge.
    145      */
    146     public static double regularizedGammaP(double a,
    147                                            double x,
    148                                            double epsilon,
    149                                            int maxIterations)
    150         throws MathException
    151     {
    152         double ret;
    153
    154         if (Double.isNaN(a) || Double.isNaN(x) || (a <= 0.0) || (x < 0.0)) {
    155             ret = Double.NaN;
    156         } else if (x == 0.0) {
    157             ret = 0.0;
    158         } else if (a > 1.0 && x > a) {
    159             // use regularizedGammaQ because it should converge faster in this
    160             // case.
    161             ret = 1.0 - regularizedGammaQ(a, x, epsilon, maxIterations);
    162         } else {
    163             // calculate series
    164             double n = 0.0; // current element index
    165             double an = 1.0 / a; // n-th element in the series
    166             double sum = an; // partial sum
    167             while (Math.abs(an) > epsilon && n < maxIterations) {
    168                 // compute next element in the series
    169                 n = n + 1.0;
    170                 an = an * (x / (a + n));
    171
    172                 // update partial sum
    173                 sum = sum + an;
    174             }
    175             if (n >= maxIterations) {
    176                 throw new ConvergenceException(
    177                     "maximum number of iterations reached");
    178             } else {
    179                 ret = Math.exp(-x + (a * Math.log(x)) - logGamma(a)) * sum;
    180             }
    181         }
    182
    183         return ret;
    184     }
    185
    186     /***
    187      * Returns the regularized gamma function Q(a, x) = 1 - P(a, x).
    188      *
    189      * @param a the a parameter.
    190      * @param x the value.
    191      * @return the regularized gamma function Q(a, x)
    192      * @throws MathException if the algorithm fails to converge.
    193      */
    194     public static double regularizedGammaQ(double a, double x)
    195         throws MathException
    196     {
    197         return regularizedGammaQ(a, x, DEFAULT_EPSILON, Integer.MAX_VALUE);
    198     }
    199
    200     /***
    201      * Returns the regularized gamma function Q(a, x) = 1 - P(a, x).
    202      *
    203      * The implementation of this method is based on:
    204      * <ul>
    205      * <li>
    206      * <a href="http://mathworld.wolfram.com/RegularizedGammaFunction.html">
    207      * Regularized Gamma Function</a>, equation (1).</li>
    208      * <li>
    209      * <a href="    http://functions.wolfram.com/GammaBetaErf/GammaRegularized/10/0003/">
    210      * Regularized incomplete gamma function: Continued fraction representations  (formula 06.08.10.0003)</a></li>
    211      * </ul>
    212      *
    213      * @param a the a parameter.
    214      * @param x the value.
    215      * @param epsilon When the absolute value of the nth item in the
    216      *                series is less than epsilon the approximation ceases
    217      *                to calculate further elements in the series.
    218      * @param maxIterations Maximum number of "iterations" to complete.
    219      * @return the regularized gamma function P(a, x)
    220      * @throws MathException if the algorithm fails to converge.
    221      */
    222     public static double regularizedGammaQ(final double a,
    223                                            double x,
    224                                            double epsilon,
    225                                            int maxIterations)
    226         throws MathException
    227     {
    228         double ret;
    229
    230         if (Double.isNaN(a) || Double.isNaN(x) || (a <= 0.0) || (x < 0.0)) {
    231             ret = Double.NaN;
    232         } else if (x == 0.0) {
    233             ret = 1.0;
    234         } else if (x < a || a <= 1.0) {
    235             // use regularizedGammaP because it should converge faster in this
    236             // case.
    237             ret = 1.0 - regularizedGammaP(a, x, epsilon, maxIterations);
    238         } else {
    239             // create continued fraction
    240             ContinuedFraction cf = new ContinuedFraction() {
    241                 protected double getA(int n, double x) {
    242                     return ((2.0 * n) + 1.0) - a + x;
    243                 }
    244
    245                 protected double getB(int n, double x) {
    246                     return n * (a - n);
    247                 }
    248             };
    249
    250             ret = 1.0 / cf.evaluate(x, epsilon, maxIterations);
    251             ret = Math.exp(-x + (a * Math.log(x)) - logGamma(a)) * ret;
    252         }
    253
    254         return ret;
    255     }
    256 }