diff options
Diffstat (limited to 'src/main/java/org/apache/commons/math/special/Gamma.java')
| -rw-r--r-- | src/main/java/org/apache/commons/math/special/Gamma.java | 339 |
1 files changed, 339 insertions, 0 deletions
diff --git a/src/main/java/org/apache/commons/math/special/Gamma.java b/src/main/java/org/apache/commons/math/special/Gamma.java new file mode 100644 index 0000000..327aa3a --- /dev/null +++ b/src/main/java/org/apache/commons/math/special/Gamma.java @@ -0,0 +1,339 @@ +/* + * 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); + } +} |