This is the mail archive of the
gsl-discuss@sourceware.org
mailing list for the GSL project.
Problems with gsl_cdf_gamma_Pinv(0.01,a,b)
- From: Benjamin Redelings <benjamin_redelings at ncsu dot edu>
- To: gsl-discuss at sourceware dot org
- Date: Mon, 22 Oct 2007 16:16:22 -0400
- Subject: Problems with gsl_cdf_gamma_Pinv(0.01,a,b)
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