1

There is a very concise algorithm for computing incomplete gamma function:

https://people.sc.fsu.edu/~jburkardt/f_src/asa147/asa147.html

However, when I look extensively over the Internet or trying to find any references I do not find a rationale for this algorithm. For example why we use log gamma function into it? Would be grateful for any comments which explain the algorithm construction and the mathematical basis for it.

rk85
  • 127

1 Answers1

2

There isn't any magic in it, and the source code does contain references: one to a standard book on special functions (a rough Internet equivalent would be this), another one to an article in "Applied Statistics", but the algorithm uses a simple and well-known series expansion (essentially the one here). It computes $$\frac{\gamma(p,x)}{\Gamma(p)}=x^p e^{-x}\sum^\infty_{k=0}\frac{x^k}{\Gamma(k+p+1)},$$ starting from the first term $$\frac{x^p e^{-x}}{\Gamma(p+1)}=\exp(p\ln x-\ln\Gamma(p+1)-x),$$ that's where the function LGAMMA is used.

  • Many thanks. I was unable to see the simple solution. It is, however, unclear why such a substitution is used. Perhaps using LGAMMA function is numerically better than using GAMMA function itself. Again, many thanks. – rk85 Sep 02 '17 at 13:56
  • Just to add. The reference to Applied Statistics does not contain any description, but the original old Fotran code: http://lib.stat.cmu.edu/apstat/147. – rk85 Sep 02 '17 at 14:08