[Bug-gsl] GSL 1.15 gsl_chisq_cdf_Pinv fails on small degrees of freedom

[Yan Zhou]

The following code fails and generate errors following it,

====================Begin of code=====================
#include <stdio.h>
#include <gsl/gsl_cdf.h>

int main( int avgc, char **argv ) {

    double nu = 0.01, P = 0.5, v;

    v = gsl_cdf_chisq_Pinv(P, nu );
    printf( "P(%f) = %f\n", P, v );

    return 0;
}
=====================end of code======================

====================Begin of error======================
gsl: gammainv.c:111: ERROR: inverse failed to converge
Default GSL error handler invoked.
Abort trap: 6
=====================end of error======================

The problem is that gsl fails for not too small degree of freedoms.  I think 
this is reported and fixed before 
(http://lists.gnu.org/archive/html/bug-gsl/2008-04/msg00036.html). The GSL 
version I was using is 1.15, compiled with Clang 3.0 (Apple) and the `make 
check` test passed successfully when installing GSL. It appears that the same 
problem remains.

And this time I don't see a simple fix for this problem as in the original 
post. Maybe this is a a flow in the algorithm itself rather the implementation 
and I would suggest taking a look into R's pgamma implementation, which seems 
to be more robust. 

Best,

Yan Zhou