[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: [Bug-gsl] Bug in gsl_cdf_chisq_Pinv for very large degrees of freedo
From: |
John Halley Gotway |
Subject: |
Re: [Bug-gsl] Bug in gsl_cdf_chisq_Pinv for very large degrees of freedom |
Date: |
Thu, 30 Oct 2008 13:20:50 -0600 |
User-agent: |
Thunderbird 2.0.0.17 (X11/20080914) |
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.
>