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

6 views (last 30 days)

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!

John D'Errico
on 2 Jun 2020

Edited: John D'Errico
on 3 Jun 2020

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.

- Use a tool that can compute in higher precision, that is either my HPF or syms.
- 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.

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.

David Goodmanson
on 2 Jun 2020

Edited: David Goodmanson
on 3 Jun 2020

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

Steven Lord
on 2 Jun 2020

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)

Opportunities for recent engineering grads.

Apply Today
## 1 Comment

## Direct link to this comment

https://ww2.mathworks.cn/matlabcentral/answers/540815-issue-with-gammainc-x-a-for-small-x-and-larger-a#comment_880934

⋮## Direct link to this comment

https://ww2.mathworks.cn/matlabcentral/answers/540815-issue-with-gammainc-x-a-for-small-x-and-larger-a#comment_880934

Sign in to comment.