1 /* 2 * Licensed to the Apache Software Foundation (ASF) under one or more 3 * contributor license agreements. See the NOTICE file distributed with 4 * this work for additional information regarding copyright ownership. 5 * The ASF licenses this file to You under the Apache License, Version 2.0 6 * (the "License"); you may not use this file except in compliance with 7 * the License. You may obtain a copy of the License at 8 * 9 * http://www.apache.org/licenses/LICENSE-2.0 10 * 11 * Unless required by applicable law or agreed to in writing, software 12 * distributed under the License is distributed on an "AS IS" BASIS, 13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 14 * See the License for the specific language governing permissions and 15 * limitations under the License. 16 */ 17 package org.apache.commons.math.special; 18 19 import org.apache.commons.math.MathException; 20 import org.apache.commons.math.MaxIterationsExceededException; 21 import org.apache.commons.math.util.ContinuedFraction; 22 import org.apache.commons.math.util.FastMath; 23 24 /** 25 * This is a utility class that provides computation methods related to the 26 * Gamma family of functions. 27 * 28 * @version $Revision: 1042510 $ $Date: 2010-12-06 02:54:18 +0100 (lun. 06 dc. 2010) $ 29 */ 30 public class Gamma { 31 32 /** 33 * <a href="http://en.wikipedia.org/wiki/Euler-Mascheroni_constant">Euler-Mascheroni constant</a> 34 * @since 2.0 35 */ 36 public static final double GAMMA = 0.577215664901532860606512090082; 37 38 /** Maximum allowed numerical error. */ 39 private static final double DEFAULT_EPSILON = 10e-15; 40 41 /** Lanczos coefficients */ 42 private static final double[] LANCZOS = 43 { 44 0.99999999999999709182, 45 57.156235665862923517, 46 -59.597960355475491248, 47 14.136097974741747174, 48 -0.49191381609762019978, 49 .33994649984811888699e-4, 50 .46523628927048575665e-4, 51 -.98374475304879564677e-4, 52 .15808870322491248884e-3, 53 -.21026444172410488319e-3, 54 .21743961811521264320e-3, 55 -.16431810653676389022e-3, 56 .84418223983852743293e-4, 57 -.26190838401581408670e-4, 58 .36899182659531622704e-5, 59 }; 60 61 /** Avoid repeated computation of log of 2 PI in logGamma */ 62 private static final double HALF_LOG_2_PI = 0.5 * FastMath.log(2.0 * FastMath.PI); 63 64 // limits for switching algorithm in digamma 65 /** C limit. */ 66 private static final double C_LIMIT = 49; 67 68 /** S limit. */ 69 private static final double S_LIMIT = 1e-5; 70 71 /** 72 * Default constructor. Prohibit instantiation. 73 */ 74 private Gamma() { 75 super(); 76 } 77 78 /** 79 * Returns the natural logarithm of the gamma function Γ(x). 80 * 81 * The implementation of this method is based on: 82 * <ul> 83 * <li><a href="http://mathworld.wolfram.com/GammaFunction.html"> 84 * Gamma Function</a>, equation (28).</li> 85 * <li><a href="http://mathworld.wolfram.com/LanczosApproximation.html"> 86 * Lanczos Approximation</a>, equations (1) through (5).</li> 87 * <li><a href="http://my.fit.edu/~gabdo/gamma.txt">Paul Godfrey, A note on 88 * the computation of the convergent Lanczos complex Gamma approximation 89 * </a></li> 90 * </ul> 91 * 92 * @param x the value. 93 * @return log(Γ(x)) 94 */ 95 public static double logGamma(double x) { 96 double ret; 97 98 if (Double.isNaN(x) || (x <= 0.0)) { 99 ret = Double.NaN; 100 } else { 101 double g = 607.0 / 128.0; 102 103 double sum = 0.0; 104 for (int i = LANCZOS.length - 1; i > 0; --i) { 105 sum = sum + (LANCZOS[i] / (x + i)); 106 } 107 sum = sum + LANCZOS[0]; 108 109 double tmp = x + g + .5; 110 ret = ((x + .5) * FastMath.log(tmp)) - tmp + 111 HALF_LOG_2_PI + FastMath.log(sum / x); 112 } 113 114 return ret; 115 } 116 117 /** 118 * Returns the regularized gamma function P(a, x). 119 * 120 * @param a the a parameter. 121 * @param x the value. 122 * @return the regularized gamma function P(a, x) 123 * @throws MathException if the algorithm fails to converge. 124 */ 125 public static double regularizedGammaP(double a, double x) 126 throws MathException 127 { 128 return regularizedGammaP(a, x, DEFAULT_EPSILON, Integer.MAX_VALUE); 129 } 130 131 132 /** 133 * Returns the regularized gamma function P(a, x). 134 * 135 * The implementation of this method is based on: 136 * <ul> 137 * <li> 138 * <a href="http://mathworld.wolfram.com/RegularizedGammaFunction.html"> 139 * Regularized Gamma Function</a>, equation (1).</li> 140 * <li> 141 * <a href="http://mathworld.wolfram.com/IncompleteGammaFunction.html"> 142 * Incomplete Gamma Function</a>, equation (4).</li> 143 * <li> 144 * <a href="http://mathworld.wolfram.com/ConfluentHypergeometricFunctionoftheFirstKind.html"> 145 * Confluent Hypergeometric Function of the First Kind</a>, equation (1). 146 * </li> 147 * </ul> 148 * 149 * @param a the a parameter. 150 * @param x the value. 151 * @param epsilon When the absolute value of the nth item in the 152 * series is less than epsilon the approximation ceases 153 * to calculate further elements in the series. 154 * @param maxIterations Maximum number of "iterations" to complete. 155 * @return the regularized gamma function P(a, x) 156 * @throws MathException if the algorithm fails to converge. 157 */ 158 public static double regularizedGammaP(double a, 159 double x, 160 double epsilon, 161 int maxIterations) 162 throws MathException 163 { 164 double ret; 165 166 if (Double.isNaN(a) || Double.isNaN(x) || (a <= 0.0) || (x < 0.0)) { 167 ret = Double.NaN; 168 } else if (x == 0.0) { 169 ret = 0.0; 170 } else if (x >= a + 1) { 171 // use regularizedGammaQ because it should converge faster in this 172 // case. 173 ret = 1.0 - regularizedGammaQ(a, x, epsilon, maxIterations); 174 } else { 175 // calculate series 176 double n = 0.0; // current element index 177 double an = 1.0 / a; // n-th element in the series 178 double sum = an; // partial sum 179 while (FastMath.abs(an/sum) > epsilon && n < maxIterations && sum < Double.POSITIVE_INFINITY) { 180 // compute next element in the series 181 n = n + 1.0; 182 an = an * (x / (a + n)); 183 184 // update partial sum 185 sum = sum + an; 186 } 187 if (n >= maxIterations) { 188 throw new MaxIterationsExceededException(maxIterations); 189 } else if (Double.isInfinite(sum)) { 190 ret = 1.0; 191 } else { 192 ret = FastMath.exp(-x + (a * FastMath.log(x)) - logGamma(a)) * sum; 193 } 194 } 195 196 return ret; 197 } 198 199 /** 200 * Returns the regularized gamma function Q(a, x) = 1 - P(a, x). 201 * 202 * @param a the a parameter. 203 * @param x the value. 204 * @return the regularized gamma function Q(a, x) 205 * @throws MathException if the algorithm fails to converge. 206 */ 207 public static double regularizedGammaQ(double a, double x) 208 throws MathException 209 { 210 return regularizedGammaQ(a, x, DEFAULT_EPSILON, Integer.MAX_VALUE); 211 } 212 213 /** 214 * Returns the regularized gamma function Q(a, x) = 1 - P(a, x). 215 * 216 * The implementation of this method is based on: 217 * <ul> 218 * <li> 219 * <a href="http://mathworld.wolfram.com/RegularizedGammaFunction.html"> 220 * Regularized Gamma Function</a>, equation (1).</li> 221 * <li> 222 * <a href="http://functions.wolfram.com/GammaBetaErf/GammaRegularized/10/0003/"> 223 * Regularized incomplete gamma function: Continued fraction representations (formula 06.08.10.0003)</a></li> 224 * </ul> 225 * 226 * @param a the a parameter. 227 * @param x the value. 228 * @param epsilon When the absolute value of the nth item in the 229 * series is less than epsilon the approximation ceases 230 * to calculate further elements in the series. 231 * @param maxIterations Maximum number of "iterations" to complete. 232 * @return the regularized gamma function P(a, x) 233 * @throws MathException if the algorithm fails to converge. 234 */ 235 public static double regularizedGammaQ(final double a, 236 double x, 237 double epsilon, 238 int maxIterations) 239 throws MathException 240 { 241 double ret; 242 243 if (Double.isNaN(a) || Double.isNaN(x) || (a <= 0.0) || (x < 0.0)) { 244 ret = Double.NaN; 245 } else if (x == 0.0) { 246 ret = 1.0; 247 } else if (x < a + 1.0) { 248 // use regularizedGammaP because it should converge faster in this 249 // case. 250 ret = 1.0 - regularizedGammaP(a, x, epsilon, maxIterations); 251 } else { 252 // create continued fraction 253 ContinuedFraction cf = new ContinuedFraction() { 254 255 @Override 256 protected double getA(int n, double x) { 257 return ((2.0 * n) + 1.0) - a + x; 258 } 259 260 @Override 261 protected double getB(int n, double x) { 262 return n * (a - n); 263 } 264 }; 265 266 ret = 1.0 / cf.evaluate(x, epsilon, maxIterations); 267 ret = FastMath.exp(-x + (a * FastMath.log(x)) - logGamma(a)) * ret; 268 } 269 270 return ret; 271 } 272 273 274 /** 275 * <p>Computes the digamma function of x.</p> 276 * 277 * <p>This is an independently written implementation of the algorithm described in 278 * Jose Bernardo, Algorithm AS 103: Psi (Digamma) Function, Applied Statistics, 1976.</p> 279 * 280 * <p>Some of the constants have been changed to increase accuracy at the moderate expense 281 * of run-time. The result should be accurate to within 10^-8 absolute tolerance for 282 * x >= 10^-5 and within 10^-8 relative tolerance for x > 0.</p> 283 * 284 * <p>Performance for large negative values of x will be quite expensive (proportional to 285 * |x|). Accuracy for negative values of x should be about 10^-8 absolute for results 286 * less than 10^5 and 10^-8 relative for results larger than that.</p> 287 * 288 * @param x the argument 289 * @return digamma(x) to within 10-8 relative or absolute error whichever is smaller 290 * @see <a href="http://en.wikipedia.org/wiki/Digamma_function"> Digamma at wikipedia </a> 291 * @see <a href="http://www.uv.es/~bernardo/1976AppStatist.pdf"> Bernardo's original article </a> 292 * @since 2.0 293 */ 294 public static double digamma(double x) { 295 if (x > 0 && x <= S_LIMIT) { 296 // use method 5 from Bernardo AS103 297 // accurate to O(x) 298 return -GAMMA - 1 / x; 299 } 300 301 if (x >= C_LIMIT) { 302 // use method 4 (accurate to O(1/x^8) 303 double inv = 1 / (x * x); 304 // 1 1 1 1 305 // log(x) - --- - ------ + ------- - ------- 306 // 2 x 12 x^2 120 x^4 252 x^6 307 return FastMath.log(x) - 0.5 / x - inv * ((1.0 / 12) + inv * (1.0 / 120 - inv / 252)); 308 } 309 310 return digamma(x + 1) - 1 / x; 311 } 312 313 /** 314 * <p>Computes the trigamma function of x. This function is derived by taking the derivative of 315 * the implementation of digamma.</p> 316 * 317 * @param x the argument 318 * @return trigamma(x) to within 10-8 relative or absolute error whichever is smaller 319 * @see <a href="http://en.wikipedia.org/wiki/Trigamma_function"> Trigamma at wikipedia </a> 320 * @see Gamma#digamma(double) 321 * @since 2.0 322 */ 323 public static double trigamma(double x) { 324 if (x > 0 && x <= S_LIMIT) { 325 return 1 / (x * x); 326 } 327 328 if (x >= C_LIMIT) { 329 double inv = 1 / (x * x); 330 // 1 1 1 1 1 331 // - + ---- + ---- - ----- + ----- 332 // x 2 3 5 7 333 // 2 x 6 x 30 x 42 x 334 return 1 / x + inv / 2 + inv / x * (1.0 / 6 - inv * (1.0 / 30 + inv / 42)); 335 } 336 337 return trigamma(x + 1) + 1 / (x * x); 338 } 339 } 340
