Issue with gammainc(x,a) for small x and larger a

Hi,

I've been having some trouble with the gammainc(x,a) function.

When x = 0.01 and a = 100, MATLAB says gammainc(x,a) = 0.

However, it should equal approx. 1e-358 (according to wolfram alpha).

How can I keep the precision in the value that close to zero?

Thank you!

[David Hill]

You could try programming the gammainc function yourself using symbolic variables; otherwise, the function cannot produce a result better than double precision.

Answer 1

[John D'Errico]

Sadly, you are limited in that no matter how hard you try, you will not be able to compute it in double precision when it is that small. It will underflow. That leaves you with several choices.

1. Use a tool that can compute in higher precision, that is either my HPF or syms.
2. Use logs, and an approximation that will be valid for small x and large a.

If you cannot store the number in double precision, you really have no choice beyond those two cases. The number 1e-358 simply does not exist in the world of double precision. It is just zero. But let me give a few examples of things I might try.

https://en.wikipedia.org/wiki/Incomplete_gamma_function

Also, we need to worry about how MATLAB defines the incomplete gamma. That is, MATLAB normalizes by gamma(a). From the help for gammainc, we see:

The incomplete gamma function is defined as:

gammainc(x,a) = 1 ./ gamma(a) .*

integral from 0 to x of t^(a-1) exp(-t) dt

So, if you want to compute the same result as say Wolfram Alpha, which probably does not have that same normalization built in, we need to be careful.

Looking at the wiki site, I see a recurrence relation for the lower incomplete gamme. I'll call it g(a,x). Thus...

g(a+1,x) = a*g(a,x) - x^a*exp(-x)

How does this help? Say we want a 16 digit approximation to g(a,x). I only want 17 digits, because I'll use gammaln at the very end. I'll do it first using HPF. The trick is to run the recurrence BACKWARDS. My computations suggest we will need to take roughly 4 terms for 16 digits of precision, but I'll take back up the reverse recurrence by 10 terms to be safe.

DefaultNumberOfDigits 17
a = 100;
x = hpf('0.01');
a0 = a + 10; % backing up by 10 terms in the recurrence
gax = hpf(0);
for ai = a0:-1:a
  gax = (gax + x^ai*exp(-x))/ai;
end
% don't forget to put in the normalization, to make
% the result consistent with gammainc in MATLAB.
gax = gax/exp(hpf(gammaln(100)));
gax

gax =
1.0609536274307748e-358

Of course, I could have done the same using the symbolic toolbox. The trick is always to start out at zero, and then run backwards. This is a common trick for many problems. Personally, I've always loved the idea of running a recurrence backwards. It often works remarkably well, at least if you do so on the right problem.

Next, I could have used a simple series approximation. In fact, this should converge pretty well, as long as x is small. We would have:

For S as the sum from 0 to infinity of the terms

x^k/gamma(a + k + 1)

g(a,x) = x^a*gamma(a)*exp(-x)*S

And of course, if we use the MATLAB incomplete gamma form, where we divide by gamma(a), it gets even simpler. I'll do this one using syms and vpa, just for fun, and to show you how you might approach it with syms.

a = sym(100);
x = sym('0.01');
S = 1/gamma(sym(a)+1);
for k = 1:50 % 50 terms should suffice
  S = S + x^k/gamma(a + k + 1);
end
gax = S*x^a*exp(-x);

vpa(gax,17)
ans =
1.0609536274307734e-358

I should be able to do this another way, if I want to compute a result in a log form, so entirely in double precision.

Here it is easiest to use gammainc with the 'scaledlower' option, which takes all the nasty stuff out.

x = 0.01;
a = 100;
gaxln = log(gammainc(x,a,'scaledlower'));
gaxln = gaxln - gammaln(a + 1) - x + a*log(x);

Done entirely in double precision, the natural log of your result is:

gaxln
gaxln =
         -824.266295139666

You cannot exponentiate it of couse as a double, because it will underflow. But we can use vpa to verify the result.

vpa(exp(sym(gaxln)),16)
ans =
1.060953627430852e-358

Which looks reasonable. If you need to work in double precision, this may be your best choice, as long as the log is sufficient.

For a very simple approximation, we can use the asymptotic behavior of the lower incomplete gamma. Again, looking at the wiki page I cited, I see the limit as:

g(a,x)/x^a --> 1/a

as x --> 0.

Therefore, we can compute an approximation to the natural log of that lower incomplete gamma as

x = 0.01;
a = 100;
gaxlnapprox = a*log(x) - log(a) - gammaln(a)
gaxlnapprox =
         -824.256394154373

which agrees to about 5 decimal digits with the value computed above. So not terribly bad as an approximation.

So you have some alternatives. With some thought, I'd bet I can come up with a few more.

[John Fullerton]

Hi John, thanks for your response! That helps a lot!

Answer 2

[David Goodmanson]

Hi John,

Matlab does not appear to have a symbolics version of the incomplete gamma function, but it's still possible to make progress. The function is

1/gamma(a) * integral {0, x0} x^(a-1) exp(-x)

In double precision, the smallest expressible number is

realmin ans = 2.225073858507201e-308

so a calculation down around 1e-358 is going to underflow to zero. However, you can do the integral separately from the gamma function and do some borrowing of the exponent to get a result.

a = 100
fun = @(x) x.^(a-1).*exp(-x)
I = integral(fun,0,.01)
   I = 9.901478680963859e-203

This appears to be fairly accurate, although a version using variable precision arithmetic might do better.

G = gamma(a)
  G = 9.332621544394404e+155

I think we can assume this one to be accurate.

You can't divide I by G because that drops below realmin, but you can do

I200 = 1e200*I;
Gm100 = 1e-100*G;
gam_inc = I200/Gm100
   gam_inc = 1.060953627430776e-58

and putting back the 1e-300 that was borrowed results in a value near 1e-358

[John Fullerton]

Thanks for your help David!

Answer 3

[Steven Lord]

Depending on what you plan to do with this tiny number, using the 'scaledlower' option for gammainc may help. I computed the overall scale factor symbolically as otherwise it would overflow.

x = 0.01;
a = 100;
L = gammainc(x, a, 'scaledlower')
scalefactor = gamma(sym(a)+1)*exp(sym(x))/sym(x)^a;
vpaScalefactor = vpa(scalefactor, 20)
result = vpa(L/scalefactor, 20)

[John Fullerton]

Thanks Steven!