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

gammagamma functiongammainc

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!

Best Answer

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.

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.

Related Solutions

MATLAB: Exp(infty1) * erfc(infty2) = NaN

You have several choices.

1. Use the symbolic Toolbox, which can handle high precision computations.

x = sym(200);
vpa(erfc(x))
ans =
4.6893592959836014980356997769923e-17375

2. Use my HPF toolbox, which can also handle high precision computations. I thought I implemented erfc, but its been a while. It should be doable though. I do know there is what looks like a nicely convergent expansion for large x for erfc. So if it turned out I did not implement erfc for large arguments, it would be pretty easy to write (I think.)

3. Use approximations, numerical analysis, etc. For example, I'd not be surprised if there is not a method to compute the log(erfc(x)), for large x. No matter what, you will still probably need a tool that can provide wide ranging computations on lots of digits, because you are adding terms there. Will you need thousands of significant digits? Or can some terms simply drop out when they are clearly insignificant? Again this is a numerical analysis question as much as anything.

Let me expand on my last comment. If the last term is insignifcant compared to the others, AND you are willing to work in double precision, you would merely compute the log of that last term. A poor approximation for the log of that expression would suffice there to know IF that term is indeed insignificant.

In the end, you will need to be careful and use good numerical methods no matter which approach you take here. Remember that the issue is not only the dynamic range of a double, but the number of significant digits it carries in the mantissa, and how many correct digits you really need in your result.

MATLAB: Finding the exponential of a very small number

Multiple ways to solve this problem. Sorry, but quad precision is not an option. There is no way to achieve quad precision, at least without recourse to some other tool, or exporting the results to a language where quad precision is available.

Brute force, without resource to the symbolic toolbox ... Use my HPF toolbox, which is free to download. I'll use 100 decimal digits of precision, even though

DefaultNumberOfDigits 100
s = hpf('-6.471606349773182e-23');
d = hpf('-6.604939683106516e-23');

See that HPF has no problem in evaluating expressions like those that you have.

exp(s)
ans =
0.9999999999999999999999352839365022681800000020940844373212284440930971466223949024568510358194049670
log(1+exp(s))
ans =
0.6931471805599453094171997634264277021655001348837763634509871205166792020337447156058633269050605702

So A is just

A = hpf('1.449056603773585e+19')*(log(1+exp(s))-log(1+exp(d)))
A =
0.000009660377358490614968553143308576941771731705532455911053800738735222612516743226005743174365407550000

As a check with the symbolic toolbox, it yields:

s = sym('-6.471606349773182e-23');
d = sym('-6.604939683106516e-23');
A = sym('1.449056603773585e+19')*(log(1+exp(s))-log(1+exp(d)))
A =
0.0000096603773584906138365242807537049

In fact, the symbolic toolbox loses some precision here, because it used limited precision to do some of those computations. And, since you only provided those numbers to about 20 significant digits, what can you expect? But it does reasonably well. The exact result, reporting it as 32 significant digits is:

0.0000096603773584906149685531433085769

Can we get that result by using only numerical methods in double precision? Brute force will of course fail, as you found.

exp(-6.471606349773182e-23)
ans =
     1
A = 1.449056603773585e+19*(log(1+exp(s))-log(1+exp(d)))
A =
     0

So in order to get a viable result for A using double precision, one needs to resort to numerical analysis. We need to understand that even if we could compute this for very tiny numbers s and d, the results log(1+exp(s)) and log(1+exp(d)) are very close to each other. Using HPF, the exact results out to 40 decimal digits are:

log(1+exp(s))
ans =
0.6931471805599453094171997634264277021655
log(1+exp(d))
ans =
0.6931471805599453094171990967597610354955

So, when we subtract the two results, we will see massive subtractive cancellation. HPF is able to give the exact results here because I did the computations in sufficiently high precision, but then only reported the first 40 digits.

Going back to double precision, even tools in MATLAB like log1p must fail, because exp(s) returns exactly 1. As far as MATLAB is concerned, with doubles for s and d, these two results are identical.

log(1+exp(s)) == log(1+exp(d))
ans =
  logical
   1

Ok, so can we do something? Yes. With some mental effort. Note that is s is TINY, then exp(s) is essentially 1. And then we lose accuracy because of the massive subtractive cancellation, because both logs are returning essentially log(2).

What if we computed the result of

log(1+exp(x)) - log(2)

Of course, we need to do this computation accurately, for very tiny x. First see that

log(1+exp(x)) - log(2) = log((1+exp(x))/2)

When we do the subtraction, those -log(2) terms will exactly cancel each other out. But they enable us to do the computation in double precision. We will need to use a Taylor series. I'll use the symbolic toolbox here, because I'm too lazy to do the algebra on paper.

syms x
taylor(log((1+exp(x))/2))
ans =
-x^4/192 + x^2/8 + x/2

For tiny values of x of magnitude smaller than eps, we can truncate that series to one term. But we can see what happens if we use 2 terms anyway.

logapprox1 = @(x) x/2;
logapprox2 = @(x) x/2 + x.^2/8;
s = -6.471606349773182e-23;
d = -6.604939683106516e-23;
format long g
A = 1.449056603773585e+19*(logapprox1(s)-logapprox1(d))
A =
      9.66037735849058e-06
A = 1.449056603773585e+19*(logapprox2(s)-logapprox2(d))
A =
      9.66037735849058e-06

As you can see, the two results are both reasonably accurate. Remember, the exact result (using scientific notation) was:

      9.6603773584906149685531433085769e-6

So our double precision computation was correct out to almost the full precision available using doubles.

The trick was to use finesse. Brute force only works if you are willing and able to do the computation using the proper tools that will support brute force. As I showed, even the symbolic toolbox missed a few digits out there at the end.

Best Answer

Related Solutions

MATLAB: Exp(infty1) * erfc(infty2) = NaN

MATLAB: Finding the exponential of a very small number

Related Question