Implementing gamma incomplete in java

Question

I'm having trouble implementing the gammainc of MATLAB in Java.

I tried doing it part by part. I used these function that uses Lanczos approximation which I got on the web for solving Γ(a):

private static double logGamma(double x) {
          double tmp = (x - 0.5) * Math.log(x + 4.5) - (x + 4.5);
          double ser = 1.0 + 76.18009173    / (x + 0)   - 86.50532033    / (x + 1)
                           + 24.01409822    / (x + 2)   -  1.231739516   / (x + 3)
                           +  0.00120858003 / (x + 4)   -  0.00000536382 / (x + 5);
          return tmp + Math.log(ser * Math.sqrt(2 * Math.PI));
       }

private static  double gamma(double x) { return Math.exp(logGamma(x)); }

Then I used Simpson's Rule on solving the integral part then I combined them to do the gammainc but the output I get is not reasonable.

The integral part can also be seen as a lower incomplete gamma function.

I'm asking for advice for a better solution.

Answer 1

Try using Apache Commons Math, which includes logGamma() and regularizedGammaP() and regularizedGammaQ().

I'm not quite sure exactly which quantity you are looking for (could you be more specific?), but some algebraic manipulation of one or two of these should get you what you want.

If you can't deal with the .jar import of Apache Commons, just include the source file for Gamma.java from it in your project (but make sure the licensing issues don't cause you any problems).