| /* |
| * Licensed to the Apache Software Foundation (ASF) under one or more |
| * contributor license agreements. See the NOTICE file distributed with |
| * this work for additional information regarding copyright ownership. |
| * The ASF licenses this file to You under the Apache License, Version 2.0 |
| * (the "License"); you may not use this file except in compliance with |
| * the License. You may obtain a copy of the License at |
| * |
| * http://www.apache.org/licenses/LICENSE-2.0 |
| * |
| * Unless required by applicable law or agreed to in writing, software |
| * distributed under the License is distributed on an "AS IS" BASIS, |
| * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
| * See the License for the specific language governing permissions and |
| * limitations under the License. |
| */ |
| package org.apache.commons.math.special; |
| |
| import org.apache.commons.math.MathException; |
| import org.apache.commons.math.MaxIterationsExceededException; |
| import org.apache.commons.math.util.ContinuedFraction; |
| import org.apache.commons.math.util.FastMath; |
| |
| /** |
| * This is a utility class that provides computation methods related to the |
| * Gamma family of functions. |
| * |
| * @version $Revision: 1042510 $ $Date: 2010-12-06 02:54:18 +0100 (lun. 06 déc. 2010) $ |
| */ |
| public class Gamma { |
| |
| /** |
| * <a href="http://en.wikipedia.org/wiki/Euler-Mascheroni_constant">Euler-Mascheroni constant</a> |
| * @since 2.0 |
| */ |
| public static final double GAMMA = 0.577215664901532860606512090082; |
| |
| /** Maximum allowed numerical error. */ |
| private static final double DEFAULT_EPSILON = 10e-15; |
| |
| /** Lanczos coefficients */ |
| private static final double[] LANCZOS = |
| { |
| 0.99999999999999709182, |
| 57.156235665862923517, |
| -59.597960355475491248, |
| 14.136097974741747174, |
| -0.49191381609762019978, |
| .33994649984811888699e-4, |
| .46523628927048575665e-4, |
| -.98374475304879564677e-4, |
| .15808870322491248884e-3, |
| -.21026444172410488319e-3, |
| .21743961811521264320e-3, |
| -.16431810653676389022e-3, |
| .84418223983852743293e-4, |
| -.26190838401581408670e-4, |
| .36899182659531622704e-5, |
| }; |
| |
| /** Avoid repeated computation of log of 2 PI in logGamma */ |
| private static final double HALF_LOG_2_PI = 0.5 * FastMath.log(2.0 * FastMath.PI); |
| |
| // limits for switching algorithm in digamma |
| /** C limit. */ |
| private static final double C_LIMIT = 49; |
| |
| /** S limit. */ |
| private static final double S_LIMIT = 1e-5; |
| |
| /** |
| * Default constructor. Prohibit instantiation. |
| */ |
| private Gamma() { |
| super(); |
| } |
| |
| /** |
| * Returns the natural logarithm of the gamma function Γ(x). |
| * |
| * The implementation of this method is based on: |
| * <ul> |
| * <li><a href="http://mathworld.wolfram.com/GammaFunction.html"> |
| * Gamma Function</a>, equation (28).</li> |
| * <li><a href="http://mathworld.wolfram.com/LanczosApproximation.html"> |
| * Lanczos Approximation</a>, equations (1) through (5).</li> |
| * <li><a href="http://my.fit.edu/~gabdo/gamma.txt">Paul Godfrey, A note on |
| * the computation of the convergent Lanczos complex Gamma approximation |
| * </a></li> |
| * </ul> |
| * |
| * @param x the value. |
| * @return log(Γ(x)) |
| */ |
| public static double logGamma(double x) { |
| double ret; |
| |
| if (Double.isNaN(x) || (x <= 0.0)) { |
| ret = Double.NaN; |
| } else { |
| double g = 607.0 / 128.0; |
| |
| double sum = 0.0; |
| for (int i = LANCZOS.length - 1; i > 0; --i) { |
| sum = sum + (LANCZOS[i] / (x + i)); |
| } |
| sum = sum + LANCZOS[0]; |
| |
| double tmp = x + g + .5; |
| ret = ((x + .5) * FastMath.log(tmp)) - tmp + |
| HALF_LOG_2_PI + FastMath.log(sum / x); |
| } |
| |
| return ret; |
| } |
| |
| /** |
| * Returns the regularized gamma function P(a, x). |
| * |
| * @param a the a parameter. |
| * @param x the value. |
| * @return the regularized gamma function P(a, x) |
| * @throws MathException if the algorithm fails to converge. |
| */ |
| public static double regularizedGammaP(double a, double x) |
| throws MathException |
| { |
| return regularizedGammaP(a, x, DEFAULT_EPSILON, Integer.MAX_VALUE); |
| } |
| |
| |
| /** |
| * Returns the regularized gamma function P(a, x). |
| * |
| * The implementation of this method is based on: |
| * <ul> |
| * <li> |
| * <a href="http://mathworld.wolfram.com/RegularizedGammaFunction.html"> |
| * Regularized Gamma Function</a>, equation (1).</li> |
| * <li> |
| * <a href="http://mathworld.wolfram.com/IncompleteGammaFunction.html"> |
| * Incomplete Gamma Function</a>, equation (4).</li> |
| * <li> |
| * <a href="http://mathworld.wolfram.com/ConfluentHypergeometricFunctionoftheFirstKind.html"> |
| * Confluent Hypergeometric Function of the First Kind</a>, equation (1). |
| * </li> |
| * </ul> |
| * |
| * @param a the a parameter. |
| * @param x the value. |
| * @param epsilon When the absolute value of the nth item in the |
| * series is less than epsilon the approximation ceases |
| * to calculate further elements in the series. |
| * @param maxIterations Maximum number of "iterations" to complete. |
| * @return the regularized gamma function P(a, x) |
| * @throws MathException if the algorithm fails to converge. |
| */ |
| public static double regularizedGammaP(double a, |
| double x, |
| double epsilon, |
| int maxIterations) |
| throws MathException |
| { |
| double ret; |
| |
| if (Double.isNaN(a) || Double.isNaN(x) || (a <= 0.0) || (x < 0.0)) { |
| ret = Double.NaN; |
| } else if (x == 0.0) { |
| ret = 0.0; |
| } else if (x >= a + 1) { |
| // use regularizedGammaQ because it should converge faster in this |
| // case. |
| ret = 1.0 - regularizedGammaQ(a, x, epsilon, maxIterations); |
| } else { |
| // calculate series |
| double n = 0.0; // current element index |
| double an = 1.0 / a; // n-th element in the series |
| double sum = an; // partial sum |
| while (FastMath.abs(an/sum) > epsilon && n < maxIterations && sum < Double.POSITIVE_INFINITY) { |
| // compute next element in the series |
| n = n + 1.0; |
| an = an * (x / (a + n)); |
| |
| // update partial sum |
| sum = sum + an; |
| } |
| if (n >= maxIterations) { |
| throw new MaxIterationsExceededException(maxIterations); |
| } else if (Double.isInfinite(sum)) { |
| ret = 1.0; |
| } else { |
| ret = FastMath.exp(-x + (a * FastMath.log(x)) - logGamma(a)) * sum; |
| } |
| } |
| |
| return ret; |
| } |
| |
| /** |
| * Returns the regularized gamma function Q(a, x) = 1 - P(a, x). |
| * |
| * @param a the a parameter. |
| * @param x the value. |
| * @return the regularized gamma function Q(a, x) |
| * @throws MathException if the algorithm fails to converge. |
| */ |
| public static double regularizedGammaQ(double a, double x) |
| throws MathException |
| { |
| return regularizedGammaQ(a, x, DEFAULT_EPSILON, Integer.MAX_VALUE); |
| } |
| |
| /** |
| * Returns the regularized gamma function Q(a, x) = 1 - P(a, x). |
| * |
| * The implementation of this method is based on: |
| * <ul> |
| * <li> |
| * <a href="http://mathworld.wolfram.com/RegularizedGammaFunction.html"> |
| * Regularized Gamma Function</a>, equation (1).</li> |
| * <li> |
| * <a href="http://functions.wolfram.com/GammaBetaErf/GammaRegularized/10/0003/"> |
| * Regularized incomplete gamma function: Continued fraction representations (formula 06.08.10.0003)</a></li> |
| * </ul> |
| * |
| * @param a the a parameter. |
| * @param x the value. |
| * @param epsilon When the absolute value of the nth item in the |
| * series is less than epsilon the approximation ceases |
| * to calculate further elements in the series. |
| * @param maxIterations Maximum number of "iterations" to complete. |
| * @return the regularized gamma function P(a, x) |
| * @throws MathException if the algorithm fails to converge. |
| */ |
| public static double regularizedGammaQ(final double a, |
| double x, |
| double epsilon, |
| int maxIterations) |
| throws MathException |
| { |
| double ret; |
| |
| if (Double.isNaN(a) || Double.isNaN(x) || (a <= 0.0) || (x < 0.0)) { |
| ret = Double.NaN; |
| } else if (x == 0.0) { |
| ret = 1.0; |
| } else if (x < a + 1.0) { |
| // use regularizedGammaP because it should converge faster in this |
| // case. |
| ret = 1.0 - regularizedGammaP(a, x, epsilon, maxIterations); |
| } else { |
| // create continued fraction |
| ContinuedFraction cf = new ContinuedFraction() { |
| |
| @Override |
| protected double getA(int n, double x) { |
| return ((2.0 * n) + 1.0) - a + x; |
| } |
| |
| @Override |
| protected double getB(int n, double x) { |
| return n * (a - n); |
| } |
| }; |
| |
| ret = 1.0 / cf.evaluate(x, epsilon, maxIterations); |
| ret = FastMath.exp(-x + (a * FastMath.log(x)) - logGamma(a)) * ret; |
| } |
| |
| return ret; |
| } |
| |
| |
| /** |
| * <p>Computes the digamma function of x.</p> |
| * |
| * <p>This is an independently written implementation of the algorithm described in |
| * Jose Bernardo, Algorithm AS 103: Psi (Digamma) Function, Applied Statistics, 1976.</p> |
| * |
| * <p>Some of the constants have been changed to increase accuracy at the moderate expense |
| * of run-time. The result should be accurate to within 10^-8 absolute tolerance for |
| * x >= 10^-5 and within 10^-8 relative tolerance for x > 0.</p> |
| * |
| * <p>Performance for large negative values of x will be quite expensive (proportional to |
| * |x|). Accuracy for negative values of x should be about 10^-8 absolute for results |
| * less than 10^5 and 10^-8 relative for results larger than that.</p> |
| * |
| * @param x the argument |
| * @return digamma(x) to within 10-8 relative or absolute error whichever is smaller |
| * @see <a href="http://en.wikipedia.org/wiki/Digamma_function"> Digamma at wikipedia </a> |
| * @see <a href="http://www.uv.es/~bernardo/1976AppStatist.pdf"> Bernardo's original article </a> |
| * @since 2.0 |
| */ |
| public static double digamma(double x) { |
| if (x > 0 && x <= S_LIMIT) { |
| // use method 5 from Bernardo AS103 |
| // accurate to O(x) |
| return -GAMMA - 1 / x; |
| } |
| |
| if (x >= C_LIMIT) { |
| // use method 4 (accurate to O(1/x^8) |
| double inv = 1 / (x * x); |
| // 1 1 1 1 |
| // log(x) - --- - ------ + ------- - ------- |
| // 2 x 12 x^2 120 x^4 252 x^6 |
| return FastMath.log(x) - 0.5 / x - inv * ((1.0 / 12) + inv * (1.0 / 120 - inv / 252)); |
| } |
| |
| return digamma(x + 1) - 1 / x; |
| } |
| |
| /** |
| * <p>Computes the trigamma function of x. This function is derived by taking the derivative of |
| * the implementation of digamma.</p> |
| * |
| * @param x the argument |
| * @return trigamma(x) to within 10-8 relative or absolute error whichever is smaller |
| * @see <a href="http://en.wikipedia.org/wiki/Trigamma_function"> Trigamma at wikipedia </a> |
| * @see Gamma#digamma(double) |
| * @since 2.0 |
| */ |
| public static double trigamma(double x) { |
| if (x > 0 && x <= S_LIMIT) { |
| return 1 / (x * x); |
| } |
| |
| if (x >= C_LIMIT) { |
| double inv = 1 / (x * x); |
| // 1 1 1 1 1 |
| // - + ---- + ---- - ----- + ----- |
| // x 2 3 5 7 |
| // 2 x 6 x 30 x 42 x |
| return 1 / x + inv / 2 + inv / x * (1.0 / 6 - inv * (1.0 / 30 + inv / 42)); |
| } |
| |
| return trigamma(x + 1) + 1 / (x * x); |
| } |
| } |
