Problems with gsl_cdf_gamma_Pinv(0.01,a,b)

Problems with gsl_cdf_gamma_Pinv(0.01,a,b)

[Benjamin Redelings]

Hi,

I have been using GSL 1.10 on a 64-bit linux system (Debian/AMD64), and 
I've been trying to figure out how to handle problems like the following:

    // works fine
    gsl_cdf_gamma_Pinv(0.01, 11.888, 1.0);

    // fails to converge
    gsl_cdf_gamma_Pinv(0.01,11.887411491530846,1.0);

Since version 1.10, this failure is a hard failure - e.g. it returns 
GSL_NAN or it terminates. 

My explorations indicate that the range on 'a' in this function is not 
very large, in the grand scheme of things.  It seems like it can be 
expected to work if 'a' is in [0.25, 10] and also for some values outside.

Q1: I'm very much interested in advice on how to handle this kind of 
situation.  Should I expect this to work in GSL, or is there some reason 
why GSL does not want this to work (e.g. speed issues)?  Should I write 
some code myself to revert to a gaussian when a > 10?

Q2: I notice that the R code for this function is more intricate and 
seems to handle values of alpha that are very high (e.g. 500) and very 
low (0.01).  At least, qgamma(0.0001,500,1) and qgamma(0.0001,0.01,1) 
give me answers inside the R program.

The R algorithm seems to use the Chi^2 distribution (which is, I 
suppose, like the Normal distribution, effectively scale-independent).  
Is there any interest in expanding gsl_cdf_gamma_Pinv( ) to handle 
similar ranges to the code in R?

https://svn.r-project.org/R/trunk/src/nmath/qgamma.c

Thanks!

-BenRI

Re: Problems with gsl_cdf_gamma_Pinv(0.01,a,b)

[Brian Gough]

At Mon, 22 Oct 2007 16:16:22 -0400,
Benjamin Redelings wrote:
> Q1: I'm very much interested in advice on how to handle this kind of 
> situation.  Should I expect this to work in GSL, or is there some reason 
> why GSL does not want this to work (e.g. speed issues)?  Should I write 
> some code myself to revert to a gaussian when a > 10?

Hello,

Thanks for the bug report, I have added a test case for that value.
You should expect this sort of thing to work in GSL -- definitely.

I added the hard failures in 1.10 as the problem was reported when I
was making the release and I didn't have time to fix it fully.  Some
of the other inverse functions were also affected so I thought it best
to report an error.

I would welcome any patches or other suggestions on the best way to
calculate this.  High precision test values to add to test.c are also
useful.

-- 
Brian Gough

GNU Scientific Library -
http://www.gnu.org/software/gsl/