Bug in incomplete gamma/gamma distribution CDF/Poisson distribution CDF.

[Zerbe, Brandon Scott]

To GSL bug report;


A colleague and I have found what I believe is a bug in numerous related 
functions related to the incomplete gamma function.  Specifically, there are 
return values of a number of related functions that are impossible --- the 
return values should fall within the range [0, 1] but don't.


Attached is a short C program and corresponding Makefile that can be used to 
generate a program that demonstrates the issue for 6 pairs of input parameters. 
I see these issues using GSL version 2.6.


In the attached program, I demonstrate this is for the functions 
gsl_cdf_gamma_Q, gsl_cdf_Poisson_P, and gsl_sf_gamma_inc_P; note, there are 
additional related functions like gsl_cdf_gamma_P, etc. that should likewise 
demonstrate this issue, but I neither examined them nor included them in the 
attached program. Specifically, for gsl_cdf_gamma_Q, the return values for the 
6 provided input parameters are:


1.100938

1.092963

1.584198

1.025459

1.118742

1.155800


When I run these input parameters using scipy.stats.gamma.cdf, I get the 
following return values that fall within the expected [0,1] range:


0.831058047952184

0.8401959763339082

0.8379612526244196

0.8384586343440682

0.8387069514229172

0.8392028379910421


So I do not believe this is a garbage-in/garbage-out situation.  I note, 
however, that the above GSL functions provide fine results on most input 
parameters I've used, so the set of problematic input parameters I found must 
go activate a branch in the code that contains the bug.


I am happy to help as much as I can in tracking down this bug --- if I were to 
do so, I would probably start at looking at the implementation of Sterling’s 
approximation within this code as that is the only branch point I can think of 
without looking into the code for these functions.  Perhaps others have a 
better guess in what could be going wrong/where to start?


Thanks!


Brandon Zerbe


#include <gsl/gsl_cdf.h>
#include <gsl/gsl_sf.h>
#include <stdio.h>

int main() {
  unsigned int ks[6] = {909922, 909954, 991071, 991073, 991074, 991076};
  double ts[6] = {909007.8, 909004.6, 990089.29, 990089.27, 990089.26,
                  990089.24};
  for (unsigned int i=0; i < 6; i++) {
    printf("%f\n", gsl_cdf_gamma_Q(ts[i], ks[i], 1.));
    printf("%f\n", gsl_cdf_poisson_P(ks[i]-1, ts[i]));
    printf("%f\n", gsl_sf_gamma_inc_P(ks[i], ts[i]));
  }
}

Makefile
Description: Makefile