[BUG] wrong output from incomplete upper gamma function

[soerenhm]

Describe the bug The (regularized) incomplete upper gamma function 'Gamma(x, a)' gives wrong values for x > 100 and for certain other special cases.

To Reproduce Here are some examples of wrong behavior:

open FSharp.Stats.SpecialFunctions
Gamma.upperIncomplete 1 infinity = 1. // Correct value: 0
Gamma.upperIncomplete 100 0 = 0.  // Correct value: 1.0
Gamma.upperIncomplete 100 200 = 1. // Correct value: ~0

OS and framework information (please complete the following information):

• OS: Windows 11
• OS Version [22000]
• .Net core SDK version 6.0.100

Additional context When Gamma.upperIncomplete fails to give the correct answer (e.g. when x > 100), the output, call it y, is identical to that produced by Gamma.lowerIncomplete; the correct should've been is 1 - y.

[bvenn]

I'm not quite into the gamma function in the moment but I tried to visualize the function and highlight the lower incomplete and upper incomplete function. In my understanding your first statement (upperIncomplete 1 infinity = 1. // Correct value: 0) is correct but I'm puzzled with latter two... The function always is nonnegative and in the case of upper(100,0) the integral is a tremendous high value or am I wrong? The same applies to upper(100,200) while here lower>upper.

[bvenn]

Two R functions I tested result in high values as well for upper(100,0)

• library expint gammainc(a,x):
  • gammainc(100,0) -> 9.332622e+155
• library pracma gammainc(x,a):
  • gammainc(0,100) -> 9.332622e+155

[soerenhm]

No you are absolutely right, upper(100, 0) is tremendously high. What can be a bit confusing is that the F# implementation is of the regularized incomplete gamma function. This means that upperIncomplete(a,x) is scaled down by gamma(a). The same is true of the lower incomplete version.

For example upper(100, 0) = 9.3...e+155 and is the same as gamma(100), which makes the regularized version identical to 1. In the case of upper(100,200), the value is also extremely high, but gamma(100) is much higher, so the regularized version is close to 0.

[bvenn]

Sorry, I absolutely missed that information :sweat_smile: Thanks for the clarification and the fix :rocket: