Re: [Bug-gsl] Bug in gsl_cdf_chisq_Pinv for very large degrees of freedom

[John Halley Gotway]

Brian,

Great thanks!

Below is a snippet of code I've written as a workaround for the time being that 
uses a normal approximation.

I found that for p = 0.05 and deg_freedom = 1263131, the result of the normal 
approximation matches the result of the gsl_cdf_chisq_Pinv routine out to the 
5th decimal place.
v(normal_approx)      = 1.26051777114457241e+06
v(gsl_cdf_chisq_Pinv) = 1.26051777113338886e+06

Thanks,
John

+++++++++++++++++++++++++++++++++++++++++++++++++

double large_n = 1263131.0;

//
// The following is a workaround for handling a GSL bug when the
// degrees of freedom are greater than 1,263,131.  As of GSL
// version 1.11, a bug in the gsl_cdf_chisq_Pinv routine causes the
// program to abort for large degrees of freedom.  Once the GSL bug
// is fixed, this check should be removed.  No warning is necessary
// since the difference in the result is out in the 5th decimal place.
//
if(deg_freedom > large_n) {

   cv_normal = normal_cdf_inv(y, 0.0, 1.0);
   x = deg_freedom*
       pow(1.0 - (2.0/(9.0*deg_freedom)) +
           cv_normal*sqrt(2.0/(9.0*deg_freedom)), 3.0);
}
else {
   x = gsl_cdf_chisq_Pinv(y, deg_freedom);
}

++++++++++++++++++++++++++++++++++++++++++++++++++



Brian Gough wrote:
> At Wed, 29 Oct 2008 11:14:48 -0600,
> John Halley Gotway wrote:
>> I've encountered a potential bug in the routine "gsl_cdf_chisq_Pinv"
>> defined in the file gsl-1.11/cdf/chisqinv.c.
>>
>> When calling this routine with the degrees of freedom > 1263131, I
>> get the following error message:
>>
>> gsl: gamma_inc.c:99: ERROR: error Default GSL error handler invoked.
>> Abort
>>
>> For degrees of freedom <= 1263131, it seems to perform fine.
>>
>> I'm running GSL version 1.11 and calling the routine as follows:
>>    v = gsl_cdf_chisq_Pinv(0.05, nu); where nu = the degrees of
>> freedom.
> 
> Thank you, I am logging this in the bug tracker.  
> 
> The routine needs switch to a different approximation for very large
> degrees of freedom, which it currently doesn't do.
>