MATLAB: Problem with gammainc function

gammainc

I am calculating incomplete gamma function using gammainc gammainc(0.04344,0,'upper') and I get the result as zero
but when I calculate same integral numerically I get result 2.6. I have also verified this results with other software, which gives the same result.
for example (using MuPAD)
double(feval(symengine,'igamma', 0, 0.043444))
gives me the result 2.6
Why I am not able to get same result with gammainc…

Best Answer

If you go by the integrals on the gammainc() documentation page and examine them with respect to approaching a = 0, you find that you have a term which is 1 / gamma(0); and since gamma(0) is infinity, 1 / gamma(0) approaches 0, which could lead you to expect that gammainc(x,0) (using lower tail) should be 0.
But when you continue the analysis, you realize that that is only true provided that the integral itself does not approach infinity.... which it does. Thus for the case a = 0, the function must proceed either by way of deeper analysis, or by way of approximation.
The symbolic math igamma() function proceeds by way of deeper analysis and converts igamma(0,x) to Ei(x), the exponential integral.
Do we know what gammainc() does? Not for sure. But an important clue is the statement in the documentation that,
For small X and A, gammainc(X,A) is approximately equal to X^A, so gammainc(0,0) = 1.
If we substitute in A=0 for that, we see that gammainc() plausibly considers gammainc(X,0) to be approximately equal to X^0 which is 1. As that is the lower tail, the upper tail would then be approximated as 1 minus that 1, giving 0. And that fits the actual returned value for gammainc(X,0,'upper'), so it probably isn't worth continuing with guessing.
So... in this boundary case, gammainc() is likely approximating and getting the wrong approximation, whereas igamma() knows the correct analytic conversion.
If you know ahead of time that your A will be 0, you should probably use the exponential integral, expint(x)
Related Question